Molecular Simulation Studies of the Diffusion of Methane, Ethane

Article pubs.acs.org/JPCC

Molecular Simulation Studies of the Diffusion of Methane, Ethane, Propane, and Propylene in ZIF‑8 Panagiotis Krokidas,† Marcelo Castier,† Salvador Moncho,‡ Edward Brothers,‡ and Ioannis G. Economou*,† †

Chemical Engineering Program and ‡Science Program, Texas A&M University at Qatar, 23874 Doha, Qatar ABSTRACT: ZIF-8 is a strong candidate for propane/ propylene separation, which is regarded as one of the most industrially demanding. Molecular simulation of this separation must account for the flexibility of the structure, which enables the adsorption and diffusion of molecules with kinetic diameter larger than the apertures of the pores. Moreover, this simulation requires modeling subtle changes since the strong sieving effect upon the mixture depends on the very small differences between propane and propylene molecular sizes (∼0.2 Å). In this work, a new force-field for the ZIF-8 structure has been developed from DFT calculations in simplified structures. The new parameter set reproduces structural properties in very good agreement with the experimental measurements reported in literature. Molecular dynamics simulations and the Widom test particle insertion method were then employed for the calculation of diffusivities, activation energies and adsorption properties of propane and propylene. The results are in agreement with experiments and demonstrate that the sieving of such a mixture is a kinetic driven separation process.

■

INTRODUCTION

Zeolitic imidazolate frameworks (ZIFs) constitute a relatively new family of materials suitable for separation processes. ZIF-8 is one of the most investigated among them because of its excellent thermal and chemical stability.1,2 Experimental measurements of adsorption and diffusion of species with kinetic diameter higher than the aperture size of the framework of ZIF-8 (estimated at 3.4 Å)1 indicate elastic behavior of the framework, as the pores adsorb compounds with sizes from methane up to propane and benzene.3−7 This behavior is attributed to the flipping motion of the ligands that form the aperture. Recent studies have focused on the performance of ZIF-8 for propane/propylene separation, which is regarded as one of the most difficult yet critical processes for the petroleum and chemical industry.8−11 Experiments using membranes have shown that the two gases exhibit a similar sorption affinity in ZIF-8, but their diffusion rates differ by 2 orders of magnitude, revealing a very promising, diffusion-driven selectivity.12,13 The two molecules are almost of the same size, which lies in the critical aperture size range of ZIF-8 (Figure 1), where the sieving effect is dominant, as has been stated in recent works.9 Propylene’s slightly smaller kinetic diameter compared to propane, as shown in Figure 2, allows it to pass more easily through the aperture and diffuse much faster than propane. These detailed molecular characteristics of ZIFs that control macroscopic properties can be elucidated using molecular simulation, which can guide the design of materials for gas separation, thus saving time and money. To the best of our knowledge, this is the first attempt to computationally study the © XXXX American Chemical Society

Figure 1. Corrected diffusivities of various gas molecules versus their molecular diameters. The black dashed line depicts the experimentally measured1 pore aperture diameter of ZIF-8. The stripe defines the effective aperture size range of ZIF-8. The plot is a rendition of the original figure found in ref 11.

propylene adsorption and the diffusion of propane and propylene in ZIF-8. Received: September 2, 2015 Revised: November 9, 2015

A

DOI: 10.1021/acs.jpcc.5b08554 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

agreement needed to accurately model the propane/propylene separation. In the present work, a new force-field is proposed and implemented in molecular dynamics (MD) and sorption properties calculations. The force-field parameters arise from fitting ab initio data for certain parameters in combination with the TraPPE force-field, which is a very accurate model for a broad range of compounds. The model is initially validated by reproducing the structure of ZIF-8 and comparing major structural properties of its framework, such as bond lengths and angles, with literature experimental values. Subsequently, the model is applied to predict the diffusivity of methane and ethane through ZIF-8. The choice of these gases is made because their dimensions lie above the aperture size of ZIF-8 (see Figure 2). Finally, we use various computational approaches for the calculation of self- and corrected-diffusivity, activation energy of diffusion, and isosteric heat of adsorption of propane and propylene in ZIF-8. We compare our results with literature data and conclude upon the validity and performance of our proposed model.

Figure 2. Kinetic diameters for (a) propane and (b) propylene. The “gray cloud” around molecules depicts the van der Waals volume of atoms (values taken from ref 11).

A key element in the simulation of subtly balanced separations such as this is the force-field used to describe interatomic potentials in the ZIF structure. This must be able to reproduce not only the structural properties of the framework itself, but also the response of the framework in the presence of diffusing species of different shape and size. This issue is most clearly demonstrated by comparing simulations where large differences in the diffusion coefficient are obtained for simulations that adopt a rigid lattice versus those that use a flexible lattice.14−16 Simulations reported recently for ZIFs used a variety of force fields. The widely used DREIDING17 and UFF18 force-fields have been used for the sorption and diffusion of Ar, CH4, H2, and CO2 molecules in ZIF-8 by Pantatosaki et al.19,20 Calculations confirmed that DREIDING is accurate for modeling the diffusion of guest molecules in a fully flexible host ZIF.19 Pantatosaki et al. employed DREIDING in additional calculations in which they investigated the effect of the rotation of linkages on the accurate prediction of gas diffusivity.20 However, Hertag et al.15 employed the DREIDING force field in their work and concluded that it provides inaccurate results regarding the flexibility of the framework. Suffriti et al.21 proposed a new set of parameters, based on the AMBER force-field22 and new ab initio calculations for the partial charges, which yielded satisfactory results for predicting the structure and, subsequently, for the simulation of diffusion of noble gases diffusion in ZIF-8.23 In this work, we implemented this force-field for our calculations and found that it underestimates the aperture slightly so that gas molecules do not diffuse fast enough compared to literature experimental values. A similar conclusion concerning the aforementioned force-field was reported by Chokbunpiam et al.24 Another variation of AMBER force-field, proposed by Hertag et al.,15 has been implemented successfully for the calculation of ethane diffusion. It was also used here and predicted a wider aperture, resulting in overestimated diffusion rates for different species compared to experimental data. The flexible framework model of Battisti et al.16 reproduces the structure of the ZIFs unit cells, but the lack of framework charges affects the prediction of the diffusion of polar molecules. Yet another recently developed force-field has reproduced the experimental structure of ZIF-8 accurately.25 Nevertheless, comparisons of reported computational predictions of CO 2 and CH 4 diffusivities with the experimental results resulted in poor agreement. This failure was attributed to the values of the framework charges used.26 Zheng et al.26 suggested a different approach to the swinging motion of the linkers, by employing an additional nonbonding interaction of harmonic form between the carbon atoms of different linkages. Although they reported an accurate description of ethane diffusion, their comparison of propane self-diffusivities with experiments did not yield the level of

■

METHODOLOGY Structure Force-Field. A new force-field was developed to adequately describe the structure of the ZIFs. The starting point was AMBER, but new parameters were added; these were derived from electronic structure calculations. Density functional theory (DFT) calculations to obtain the parameters were performed with the Gaussian 09 suite of programs.27 Specifically, the hybrid B3LYP density functional28,29 with the 6-311g++(2d,2p) basis set was used throughout this work. This particular level of theory was selected as it has been used before in related studies on zinc bearing metallo-organic frameworks.30 A very dense integration grid has been chosen, namely, the “ultrafine” pruned grid with 99 radial shells of 590 angular points. A dianionic model complex with one central metal atom, [Zn(N-mIM)4]2− was used for the ab initio calculation of the parameters (Figure 3).

Figure 3. DFT optimized geometry of the [Zn(N-mIM)4]2− model.

The [Zn(N-mIM)4]2− model geometry was fully optimized using DFT and the proposed equilibrium values for the geometrical parameters correspond to the most stable of the considered conformers. Force constants for bonds and angles were computed from the relevant terms of the Hessian matrix in Cartesian coordinates, using the formulas proposed by Seminario31 and are shown in Table 1. For the description of the torsional potential, a Fourier-type AMBER expression was implemented based on the rotation barrier, and the relevant parameters are shown in Table 2. Most of the rotation parameters were obtained from the default AMBER parameter set. The rotation barrier for the metal−ligand bond was B


Article


Table 1. New Structure Force-Field Parameters for Bond Stretching and Bond Angle Bending for ZIF-8 Using the AMBER Potential bond

l0 (Å)

kl (kJ/mol/nm2)

angle

θ0

kθ (kJ/mol/rad2)

Zn−N N−C2 N−C1 C1−C1 C1−H1 C2−C3 C3−H2

2.048 1.360 1.376 1.375 1.077 1.498 1.091

52802.1 257818.1 253048.3 339991.8 327690.9 203760.8 286855.0

N−Zn−N Zn−N−C2 Zn−N−C1 C1−N−C2 C1−C1−N C1−C1−H1 C2−C3−H2 H2−C3−H2 N−C2−N N−C2−C3 N−C1−H1

109.5 130.3 125.1 104.5 107.9 130.6 110.8 108.1 113.8 123.1 121.5

296.23 462.74 475.30 1077.80 909.61 552.29 565.68 317.98 955.63 958.97 549.78

Table 2. New Structure Force-Field Parameters for the Torsional Potential of ZIF-8 Using the AMBER Potential

Table 3. New Structure Force-Field Parameters (Partial Charges) for ZIFs Using the AMBER Potential

dihedral

φ0 (degrees)

m

kφ (kJ/mol)

source

atom type

partial charge (e)

N−Co−N−C2 N−Co−N−C1 N−C1−C1−N N−C1−C1−H1 C1−C1−N−Co C1−C1−N−C2 C3−C2−N−Co C3−C2−N−C1

180 180 180 180 180 180 180 180

6 6 2 2 2 2 2 2

10.46 10.46 90.0 90.0 25.1 25.1 41.8 41.8

this work this work AMBER AMBER AMBER AMBER AMBER AMBER

Zn N C1 C2 H1 C3 H2

1.3429 −0.6822 −0.0622 0.7551 0.0912 −0.2697 0.0499

estimated (using the same DFT method) by freezing different values for the dihedral angle. The bond stretching is described by the following equation: kl (l − l0)2 2 while the angle bending is given by U stretch(l) =

for all atoms. Partial charges in the force-field were obtained by averaging the values for different fully optimized geometries of the complex. They correspond to different conformers with comparable energy (relative energies within 8 kJ/mol) obtained through the rotation of the metal−ligand bonds. We used the average value because it provides a more flexible description of how charge redistributes in different conformations around the metal center. We employed the scaled 1−4 approach, which is a common choice for handling interactions between neighboring atoms in the same chain;21,22 interactions between atoms separated by one or two bonds are excluded. For the interactions between atoms separated by three bonds, the Coulombic potential parameter values are divided by 2, while the van der Waals potential parameters are divided by 1.2. Guest−Guest and Guest−Framework Interactions. The transferable potential for phase equilibria (TraPPE) force-field,36 which is accurate for thermodynamic and transport properties of hydrocarbons in mixtures37,38 and ZIFs,24 was employed for the interactions of all diffusing species in the system. The united atom (UA) approach was adopted, so that hydrogen atoms are lumped into the carbon atoms. In TraPPE, the length of C−C single bonds is 1.54 Å, while the length of CC double bonds is 1.33 Å. Harmonic potentials describe the change in angles between three interconnected atoms according to eq 2, while Lennard-Jones interactions account for the nonbonded inter/intramolecular interactions, according to eq 4. The values for all the parameters of the TraPPE force-field for the guest molecules used in this work can be found in Tables 4 and 5. The Lorentz−Berthelot mixing rules were used for the nonbonded interactions between different hydrocarbon atoms, as well as for the cross-interactions between guest and host atoms.

(1)

kθ (θ − θ0)2 (2) 2 The torsional potential in the ZIF-8 framework is given by U bend(θ ) =

U torsion(φ) = kφ[1 + cos(mφ − φ0)2 ]]

(3)

Finally, nonbonded interactions are calculated based on the Lennard-Jones potential, according to the expression: ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij U (rij) = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠ LJ

(4)

All Lennard-Jones parameters for the atoms in the ZIF-8 framework were taken from ref 15. Partial charges shown in Table 3 were calculated using the RESP procedure within the Antechamber program.32,33 The electrostatic potential was calculated and sampled with Gaussian, using the Merz−Singh−Kollman procedure34,35 with UFF radii. During the RESP procedure, all charges in the model compound were optimized, restraining equivalent atoms in the ligands to have the same energy. Dianionic model ML−2 4 implies formal charges of +1 for the metal and −1 for the ligands, but the computed partial charges are different, causing the ML2 fragment to be charged (non-neutral). This was modified by reducing the charge in the ligand’s atoms to a neutral ML2 fragment; the excess charge was reduced equally C


Article


For the validation of the force-field for ZIF-8, 5 ns long simulations were used to obtain equilibrated structural properties that were compared to the experimental structure from the literature. The system was then loaded with the desired number of gas molecules (CH4, C2H6, C3H8, or C3H6) and the diffusion of the species was examined based on MD runs that varied from 10 to 100 ns. The mean square displacement of species was used to calculate the self-diffusion coefficient based on the Einstein relation for a threedimensional system:42,43

Table 4. Lennard-Jones Parameters for the Various UAs in the TraPPE Potential

1 Ds = lim t →∞ 6t

1 N

N

∑ | ri (⃗ t ) − ri (⃗ t0)|2

(5)

i=1

where ri⃗ (t) denotes the position vector of the diffusing molecule in the ZIF-8 supercell at time t, Ds is an average over all N diffusing molecules whose positions are given by the vector ri⃗ (t) and for various time origins t0. In order to have a direct comparison with the experimental results of propane and propylene diffusion in ZIF-8, we calculated the corrected diffusivity, D0, which connects the transport diffusivity with the so-called thermodynamic correction factor in the Darken equation44 through the expression: Dt(c) = D0(c)((∂ ln f)/(∂ ln c))T, where Dt(c) is the transport diffusivity and f and c are the fugacity and concentration of adsorbed species, respectively. We computed D0 for propane and propylene using a methodology developed by Theodorou et al.,45 which provided satisfactory results in simulations of gas diffusion in zeolites.46 The method makes use of an equation similar to eq 5:

Table 5. Parameters of the TraPPE Potential for the Harmonic Angle Potential of Propane and Propylene

■

SIMULATION DETAILS Construction of ZIF-8 in the Atomistic Scale. The experimental data from Park et al.1 were used to construct the unit cell of ZIF-8. It consists of a cubic unit cell of length 16.991 Å that belongs to the I4̅3m space group, which results in a sodalite topology with the following composition: C96H120N48Zn12, where the Zn atom is tetrahedrally connected with four N atoms. The N atoms are bridged by 2-methylimidazolate (mIM) atom units shown in Figure 4, which results in the Zn(N-mIM)4 unit.

N

D0 =

1 1 lim ⟨|∑ ( ri (⃗ t ) − ri (⃗ t0))|2 ⟩ 6N t →∞ t i = 1

(6)

Widom Test Particle Insertion Method. The isosteric heat of adsorption of propane and propylene was calculated by applying the Widom test particle insertion method,47 which has been employed successfully in the past for the diffusion of gas molecules in ZIF-8.20 The technique incorporates the insertion of a single sorbate molecule in the system (ZIF framework), at fixed volume, V, and temperature, T. Different positions of the guest molecule are tried for a large number of conformations of the ZIF structure, generated by the MD simulations. In the case of the isochoric−isothermal (NVT) ensemble, this method provides the excess chemical potential, μex i , of component i, according to the expression:48,49

Figure 4. Basic Zn(N-mIM)4 unit of ZIF-8 framework.

1 inter )⟩Widom, N , V , T ] μiex = − ln[⟨exp( −βUtest β

(7)

where β = 1/kBT and Uinter test is the intermolecular energy of the test particle when it is inserted in the system. The brackets indicate averaging over all ZIF framework configurations and spatial averaging over all the degrees of freedom of the inserted guest molecule (translational and rotational). In this work, 1000 frames were generated by MD simulations on the empty ZIF-8 supercell. In each frame 50,000 insertions of the test particle (C3H8 or C3H6 molecule) were tried in order to get adequate statistics on the conformations of the guest−host system. The energy of the system was calculated in each insertion and an average over these different conformations provided the excess chemical potential, μex i , of the inserted species according to eq 7.

A super cell of 2 × 2 × 2 unit cells, shown in Figure 5, was used in the simulation box. The supercell and the ZnN4 tetrahedra have been visualized with the use of the Discovery Studio Visualizer.39 Molecular Dynamics Simulations. Cubic boundary conditions were applied in all three directions of the box and a cutoff of 13 Å was imposed for the van der Waals interactions. MD simulations were carried out in the NVT ensemble using the Nosé−Hoover thermostat,40,41 with a step of 1.0 fs. For the electrostatic interactions, the particle mesh Ewald method (PME) was used. Interactions between bonded atoms were treated with the 1−4 convention mentioned above, which is a standard approach in the literature for these systems.21,22· D


Article


Figure 5. Snapshot of the reconstructed supercell of ZIF-8 as used in our simulations. Gray balls depict hydrogen atoms. The other atoms follow the color sequence of Figure 4.

temperature was set at 258 K, to allow a direct comparison with available experimental data.1 In Table 6, experimental data

In the limit of very low pressure, the sorption capacity exhibits a linear dependency on fugacity, which can be expressed according to

KH = lim

xi → 0

fi xi

Table 6. Representative Bond Lengths and Bond Angles of the ZIF-8 framework from XRD Measurements1 and MD Simulations

(8)

where KH is the Henry constant, and f i and xi are the fugacity and the mole fraction of sorbed species i, respectively. The Henry constant of a guest molecule in the ZIF framework can be evaluated by the following expression:48 ⎛ρ ⎞ KH = lim ⎜ ZIF exp(βμiex )⎟ xi → 0⎝ β ⎠

(9)

where ρZIF is the molecular density of the empty framework, which is the number of ZIF-8 molecules in the simulation cell over the cell volume. The isosteric heat of adsorption, Q, is given by ⎛ Q ⎞ ⎟ KH = K 0 exp⎜ − ⎝ RT ⎠

(10)

where K0 is a constant and R is the gas constant. The isosteric heat was obtained by performing Widom insertion calculations at various temperatures and plotting ln KH as a function of the inverse temperature (1/T). The slope of the Arrhenius plot provides Q/R. MD simulations and Widom test particle insertions were carried out using the GROMACS open-source molecular simulation platform (version 4.6.5).50 The assignment of the force-field types and parameters to each atom of the system requires the construction of a file that, in the case of a complex chemical system such as the ZIF-8 supercell, requires an automated procedure. For this reason, we used the g_x2top routine of GROMACS, for the production of topologies for our calculations.

and MD predictions for various structural properties are shown. The agreement is very good in all cases. Deviation between experimental data and the MD simulation is lower than 0.8% for bond lengths and 0.5% for bond angles. An important structure property is the aperture, which gives access to the large cavities of the system. In this work, we performed MD simulations using two additional force-fields, which have been successfully employed in the literature to simulate gas diffusion in ZIF-8.16,21 In both cases, we used the scaling factor for the 1−4 interactions mentioned before. In all cases, we calculated the diameter of the aperture following the procedure proposed by Chokbunpiam et al.24 Two hydrogen atoms chosen are separated by a distance r1. A third hydrogen atom is picked and the diameter of the circle where these three points lie on is calculated from the formula:

■

RESULTS AND DISCUSSION Validation of the Force-Field. In order to validate the accuracy of the force-field to represent the ZIF-8 structure, MD simulations in the NVT ensemble were performed. The

d=2

E

(r1r2r3) (r1 + r2 + r3)(r2 + r3 − r1)(r3 + r2 − r2)(r1 + r2 − r3) (11) DOI: 10.1021/acs.jpcc.5b08554 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article


general, if an atomistic reconstruction of the material on the computer is feasible (from experimental data, etc.), then the NVT/NVE ensemble is the common choice for investigating the various questions from experiment about diffusion in such materials. Normally, a model capable of predicting structural changes upon pressure variation is needed for reproducing the structure in conditions outside available experimental data. This is a task to be treated in a future work. The force field was used subsequently for the calculation of the diffusion coefficient of light alkanes and propylene. Initially, methane and ethane were examined because experimental data are readily available,5,4 and MD simulations predict selfdiffusivities that agree well with these data.15,24 The loading used matches the experimental literature, with five methane molecules per unit cell and eight ethane molecules per unit cell. These occupancies resulted in 40 CH4 and 64 C2H6 molecules in the supercell. Four to five runs of 20 and 30 ns for methane and ethane, respectively, were adequate to get reliable statistics for the calculations of diffusion coefficients. The diffusivity was calculated from the mean square displacement curves. Representative mean square displacements as a function of time for both species are shown in Figure 7.

The same procedure was repeated for all combinations of hydrogen atoms, retaining the smallest circle found, as shown in Figure 6(a), which is regarded as the geometrical circum-

Figure 6. Aperture diameter of ZIF-8: (a) the calculated geometrical circumference (green dashed circle) leads to (b) the effective aperture diameter (red dashed circle) by subtraction of the hydrogen van der Waals radius.

ference. The effective diameter is the geometrical diameter minus the van der Waals radius of hydrogen atoms, which here is estimated to be 1.2 Å (see Figure 6b). This is a common approach to depict the effective diameter of the aperture.1 In Table 7, experimental data and MD simulation results are shown for the aperture diameter. Table 7. Aperture Diameter from Experimental Data and MD Simulations Using Various Force-Fields MD simulation experimental d(Å)

∼3.4

1

new force-field

from ref 15

from ref 22

3.44

3.59

3.21

The value of 3.44 Å agrees very well with the experiments of Park et al.1 The force-field from the work of Hertag et al., when implemented in our calculations, gave a value of 3.59 Å, which is bigger than the experimental value and resulted in an overestimation of diffusion of species in our calculations. A possible explanation for this discrepancy can be attributed to the 1−4 interactions, for which the authors do not state explicitly their approach. However, implementation of the force-field proposed by Zheng et al.21 in our MD simulations gave an underestimated aperture size (3.20 Å). The sieving effect of ZIF-8 on propane/propylene, which depends on the slight difference on the two molecule sizes, can only be reproduced accurately if the aperture size is predicted reliably by simulations. All MD simulations in this work were performed in the NVT statistical ensemble. The majority of flexible metal−organic framework (MOF) models are very sensitive to the application of pressure, having a tendency to shrink or even collapse the framework in the course of NPT simulations, despite their satisfactory operation in NVT calculations. Pantatosaki et al.20 reported reasonable results by employing DREIDING forcefield17 in their NVT simulations. However, they had to modify the angle, stretch, and torsion parameters in the ZnN4 tetrahedra when they switched to NPT calculations, in order to prevent the shrinkage of the framework. The force-field proposed by Hertag et al., which gives good results in the diffusion of gas molecules at fixed NVE,15,24 causes an unreasonable shrinkage of the unit cell when employed in the NPT ensemble. Few proposed force-fields in literature demonstrate a realistic response to temperature and pressure changes (thermal expansion coefficient and bulk modulus). In

Figure 7. Representative mean square displacement of methane () and ethane (···) as function of time in ZIF-8 at 298 K.

Self-diffusion coefficients from MD simulations for both species are in very good agreement with experimental data, as shown in Table 8. The experimental values for CH4 are transport diffusivities (DT) measured with infrared microscopy at very low loading4 (where they coincide with self-diffusivities) and self-diffusivities through pulsed-field gradient nuclear magnetic resonance technique in a loading that corresponds to 5 CH4 mol/u.c.20 The experimental value of ethane Table 8. Self-Diffusivities of CH4 and C2H6 from Our Simulations and Comparison with Experiments and Simulations from Literature DCH4 (m2/s) experimental data

4,5,20

MD simulations (this work) MD simulations (literature)14,24 F

−10

∼1 × 10 ∼1.5 × 10−10 (1.3 ± 0.1) × 10−10 (1.6 ± 0.2) × 10−10

DC2H6 (m2/s) ∼2 × 10−11 (1.0 ± 0.1) × 10−11 6 × 10−11


Article


systematically higher than their self-diffusivities. Ds depends on the loading of the diffusing molecules in ZIF, while D0 is regarded as concentration independent (although this is a matter of discussion and simulations of single component diffusion in zeolites have shown that this approximation does not hold always51). D0 and Ds become equal in the limit of very low adsorbate concentration.45 Regarding the above, it would be interesting to study the dependency of these quantities on the loading of propane and propylene. Work is underway toward this. Our calculations allow the estimation of the activation energy of diffusion. The values of activation energy for both species can be found in Table 12, together with experimental values reported by Liu et al.12 and Hara et al.13 The extracted values for activation energy of diffusion (Eα) from the Arrhenius analysis (Figure 9) agree well with the experimental values reported in literature (Table 12). Propane’s higher values for activation energy are a sign of its hindered motion. The diffusing propane molecules have to overcome a high energy barrier in order to diffuse in ZIF-8. However, propylene’s lower value shows that the transitions along each trajectory as it passes through the structural “bottlenecks” are much smoother and, therefore, propylene diffuses faster. The comparison is a confirmation of our model’s efficiency in describing the interactions that govern the system as a whole. A weak point in the majority of molecular mechanics models can be the calculation of the “guest−host” interactions, which is based on Lorentz−Berthelot or other combining rules. In this approach, parameters from different force fields are combined (in this work TraPPE with the new force field from DFT calculations). The TraPPE force-field was developed mainly for the calculation of equilibrium properties of gases. The question arises as to whether a more detailed, full atom description of the guest molecules is needed in the polar environment of a ZIF framework. In such an approach, each hydrogen atom is described explicitly, and partial charges are assigned to the atoms of the hydrocarbon molecules. In the context of isolating the guest−host interactions, the isosteric heat of adsorption for both species in ZIF-8 was calculated. This property is directly connected with the sorption affinity of guest molecules in a porous medium at the Henry regime, where the guest−host interactions govern the sorption. Therefore, a successful calculation of isosteric heat of adsorption is a proof of valid guest−host cross-interactions. The calculation of Henry constant is carried out by employing the Widom insertion method as described in the methodology section. The calculations extended over temperatures ranging from room temperature up to 450 K, and then the Arrhenius plot of lnK vs 1/T was constructed as shown in Figure 10, where the curve demonstrates an almost linear behavior. From the slope, we extract the isosteric heat of adsorption according to eq 10. A comparison of our results with experimental data is shown in Table 13. In all cases, the two species exhibit the same sorption affinity with ZIF-8. Consequently, separation of the mixture is a diffusion-driven process.

corresponds to the self-diffusivity at a loading of 8 mol/u.c, as measured by PFG NMR experiments.5 It should be mentioned that the molecular diameters of both methane and ethane (3.8 and 4.0 Å, respectively) are larger than the aperture of 3.44 Å. Since our model predicts accurately the diameter for the pore aperture, diffusion of those species means that our model describes the framework properly. The main aim of this study is to model accurately propane and propylene diffusion in ZIF-8. The reconstructed super cell was loaded with 110 C3H8 and 110 C3H6 molecules, for the calculation of the self- and corrected diffusivities. These loadings correspond to the maximum uptake of C3H8 and C3H6 in ZIF-8 at 310 K, as measured by Zhang et al. (∼5.0 mmol/gr at ∼4 bar).11 We chose the highest available loading at 310 K in order to get better statistics since propane diffuses very slowly in ZIF-8. Liu et al.12 measured pure gas diffusivities of propane and propylene in ZIF-8 with steady state membrane permeation. Experiments by Zhang et al. reported corrected diffusivities (D0) obtained from kinetic uptake rate measurements at 310 K. In the same work, they also reported steady state measurements11 at the same temperature. Self-diffusivities and corrected diffusivities for both species from molecular simulations are in satisfactory agreement with experimental values reported in literature as shown in Table 9. Table 9. Self- and Corrected Diffusivities of C3H8 and C3H6 from Molecular Simulations and Comparison with Experiments from Literature at 310 K DC3H8 (m2/s) gas permeance measurements12 kinetic uptake measurements11 steady state measurements11 MD simulations, selfdiffusivities MD simulations, corrected diffusivities a

DC3H6 (m2/s)

3.99 × 10−14

1.25 × 10−12

2.0 × 10−14

2.9 × 10−12

(1.7 ± 0.8) × 10−14

(1.6 ± 0.3) × 10−12

(5.1 ± 0.9) × 10−14a

(1.6 ± 0.1) × 10−12

(8 ± 1) × 10−14a

(2.0 ± 0.7) × 10−12

extrapolated from higher temperature values using an Arrhenius plot.

This is the first validation of the delicate sieving effect of ZIF-8 from MD simulations. Although propane and propylene differ marginally in size, their diffusion at low temperatures differs by 2 orders of magnitude. Representative results for the mean square displacement for the two diffusivities (Ds and D0) at various temperatures are shown in Figure 8. In all cases, D0 is higher than Ds. These properties make the material unique with respect to its separation performance. Propane diffuses very slowly at low temperatures. As a result, MD simulations cannot provide a reliable, direct prediction of the diffusion coefficient in room temperature. For this purpose, MD simulations at higher temperatures were performed and an Arrhenius representation was used to estimate both the selfdiffusivity and corrected diffusivity at 310 K. The values for all the calculated self-diffusivities and corrected diffusivities at various temperatures are provided in Tables 10 and 11, respectively. Corrected diffusivities are measured experimentally using various methods, such as through kinetic uptake measurements.11 Although the calculated corrected and self-diffusivities are of the same order of magnitude over the entire temperature range, the corrected diffusivities of propane and propylene are

■

CONCLUSIONS A new force-field for the ZIF-8 structure has been developed by combining AMBER force field for the Lennard-Jones interactions and ab initio calculations for the bonded interactions and the partial charges of the framework. The force field was used for the detailed modeling and simulation of G


Article


Figure 8. Representative mean square displacement for Ds () and Do (···) of propane and propylene as a function of time in ZIF-8 at three different temperatures.

Table 10. Self-Diffusivities of C3H8 and C3H6 As-Calculated by Simulations of This Work at Various Temperatures T (K) 310 330 350 370 400 430 470 530

DC3H8 (m2/s)

(1.6 ± 0.5) × 10−13 (2.7 ± 0.8) × 10−13 (6 ± 0.8) × 10−13 (1.3 ± 0.1) × 10−12 (2.1 ± 0.2) × 10−12 (5.09 ± 0.09) × 10−12

Table 12. Activation Energy of Diffusion As-Calculated by Simulation from This Work, along with Experimentally Reported Values

DC3H6 (m2/s) (1.6 ± 0.1) × (3.1 ± 0.2) × (5.1 ± 0.4) × (7.2 ± 0.9) × (1.20 ± 0.05) (1.6 ± 0.1) × (2.5 ± 0.2) × (4.5 ± 0.3) ×

10−12 10−12 10−12 10−12 × 10−11 10−11 10−11 10−11

experiments12 experiments13 MD simulations

EC3H8 (kJ/mol)

EC3H6 (kJ/mol)

38.8 26.6 30.0 ± 1.0

12.7 15.1 18.0 ± 1.1

Table 11. Corrected Diffusivities of C3H8 and C3H6 AsCalculated by Simulations of This Work at Various Temperatures T (K) 310 330 350 370 400 430 470 530

DC3H8 (m2/s)

DC3H6 (m2/s)

(2.9 ± 0.9) × 10−13 (3.7 ± 0.8) × 10−13 (8 ± 1) × 10−13 (1.4 ± 0.8) × 10−12 (2.5 ± 0.8) × 10−12 (5.6 ± 0.5) × 10−12

(2.0 ± 0.7) × 10−12 (3.5 ± 0.8) × 10−12 (5.4 ± 0.7) × 10−12 (8 ± 2) × 10−12 (1.5 ± 0.5) × 10−11 (1.9 ± 0.6) × 10−11 (3.2 ± 1.0) × 10−11 (5.0 ± 0.8) × 10−11

Figure 9. Arrhenius plot of self-diffusivities for propane (■) and propylene (●).

propane/propylene separation. Important structural parameters, such as bond lengths and angles of the framework, were predicted with high accuracy. The critical aperture size, which governs the sieving efficiency of the material, compared well with the experimentally measured diameter. The flexibility of the framework with the new force-field was modeled in a very

detailed manner, allowing methane and ethane molecules to pass through the apertures. The agreement of the estimated self- and corrected-diffusivities with experiments in the literature enhanced the efficacy of the proposed model in achieving proper flexible behavior. H


Article


Stability of Zeolitic Imidazolate Frameworks. Coord. Chem. Rev. 2011, 255, 1791−1823. (3) Bux, H.; Feldhoff, A.; Cravillon, J.; Wiebcke, M.; Li, Y. S.; Caro, J. Oriented Zeolitic Imidazolate Framework-8 Membrane with Sharp H2/C3H8Molecular Sieve Separation. Chem. Mater. 2011, 23, 2262− 2269. (4) Bux, H.; Chmelik, C.; Krishna, R.; Caro, J. Ethene/Ethane Separation by the MOF Membrane ZIF-8: Molecular Correlation of Permeation, Adsorption, Diffusion. J. Membr. Sci. 2011, 369, 284−289. (5) Bux, H.; Chmelik, C.; Van Baten, J. M.; Krishna, R.; Caro, J. Novel MOF-Membrane for Molecular Sieving Predicted by IRDiffusion Studies and Molecular Modeling. Adv. Mater. 2010, 22, 4741−4743. (6) Gücüyener, C.; Van Den Bergh, J.; Gascon, J.; Kapteijn, F. Ethane/Ethene Separation Turned on Its Head: Selective Ethane Adsorption on the Metal−Organic Framework ZIF-7 through a GateOpening Mechanism. J. Am. Chem. Soc. 2010, 132, 17704−17706. (7) Diestel, L.; Bux, H.; Wachsmuth, D.; Caro, J. Pervaporation Atudies of n-hexane, Benzene, Mesitylene and their Mixtures on Zeoolitic Imidazolate Framework-8 Membranes. Microporous Mesoporous Mater. 2012, 164, 288−293. (8) M. Robert, A. Handbook of Petroleum Refining Processes; McGrawHill: New York, 1986. (9) Ma, X. L.; Lin, B. K.; Wei, X. T.; Kniep, J.; Lin, Y. S. GammaAlumina Supported Carbon Molecular Sieve Membrane for Propylene/Propane Separation. Ind. Eng. Chem. Res. 2013, 52, 4297−4305. (10) Bernardo, P.; Drioli, E.; Golemme, G. Membrane Gas Separation: A Review/State of the Art. Ind. Eng. Chem. Res. 2009, 48, 4638−4663. (11) Zhang, C.; Lively, R. P.; Zhang, K.; Johnson, J. R.; Karvan, O.; Koros, W. Unexpected Molecular Sieving Properties of Zeolitic Imidazolate Framework-8. J. Phys. Chem. Lett. 2012, 3, 2130−2134. (12) Liu, D.; Xiaoli, M.; Hongxia, X.; Lin, Y. S. Gas Transport Properties and Propylene/Propane Separation Characteristics of ZIF-8 Membranes. J. Membr. Sci. 2014, 451, 85−93. (13) Hara, N.; Yoshimune, M.; Hideyki, N.; Haraya, K.; Hara, S.; Yamaguchi, T. Diffusive Separation of Propylene/Propane with ZIF-8 Membranes. J. Membr. Sci. 2014, 450, 215−223. (14) Seehamart, K.; Nanok, T.; Krishna, R.; van Baten, J. M.; Remsungnen, T.; Fritzsche, S. A Molecular Dynamics Investigation of the Influence of Framework Flexibility on Self-diffusivity of Ethane in Zn(tbip) Frameworks. Microporous Mesoporous Mater. 2009, 125, 97− 100. (15) Hertag, L.; Bux, H.; Caro, J.; Chmelik, C.; Remsungen, T.; Knauth, M.; Fritzsche, S. Diffusion of CH4 and H2 in ZIF-8. J. Membr. Sci. 2011, 377, 36−41. (16) Battisti, A.; Taioli, S.; Garberoglio, G. Zeolitic Imidazolate Frameworks for Separation of Binary Mixtures of CO2, CH4, N2 and H2: A Computer Simulation Investigation. Microporous Mesoporous Mater. 2011, 143, 46−53. (17) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: a Generic Force Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (18) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (19) Pantatosaki, E.; Pazzona, F. G.; Megariotis, G.; F; Papadopoulos, G. K. Atomistic Simulation Studies on the Dynamics and Thermodynamics of Nonpolar Molecules within the Zeolite Imidazolate Framework-8. J. Phys. Chem. B 2010, 114, 2493−2503. (20) Pantatosaki, E.; Megariotis, G.; Pusch, A. K.; Chmelik, C.; Stallmach, F.; Papadopoulos, G. K. On the Impact of Sorbent Mobility on the Sorbed Phase Equilibria and Dynamics: A Study of Methane and Carbon Dioxide within the Zeolite Imidazolate Framework-8. J. Phys. Chem. C 2012, 116, 201−207.

Figure 10. Arrhenius plot of Henry constant for C3H8 (■) and C3H6 (□).

Table 13. Comparison of Isosteric Heat of Adsorption from Simulations of This Work with Experimental Values from Literature ΔHC3H8 (kJ/mol)

ΔHC3H6 (kJ/mol)

18.9 30 25.9

18.4 34 26.3

13

experiment experiment52 MD simulation

The properties for the propane/propylene separation were simulated and yielded results that agree with literature experiments. The difference in the activation energy of diffusion of the two species characterizes this separation as being kinetic driven. This is corroborated by the calculation of the isosteric heats of adsorption, which is almost the same for both species, as a consequence of the similar affinity of molecules with the surrounding walls of the pores. Modification of the framework, by replacement of specific atoms or units (like the linkages), is expected to enhance the separation efficiency of ZIF-8. Such problems can be addressed with proper adjustment and modification of the force field, which will allow the treatment of systems with significant structural changes.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This publication was made possible by NPRP grant number 7042-2-021 from the Qatar National Research Fund (a member of the Qatar Foundation). The statements made herein are solely the responsibility of the authors. We are grateful to the High Performance Computing Center of Texas A&M University at Qatar for generous resource allocation.

■

REFERENCES

(1) Park, K. S.; Ni, Z.; Cote, A. P.; Choi, J. Y.; Huang, R.; UribeRomo, F. J.; Chae, H. K.; O’Keeffe, M.; Yaghi, O. M. Exceptional Chemical and Thermal Stability of Zeolitic Imidazolate Frameworks. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 10186−10911. (2) Li, J. R.; Ma, Y.; McCarthy, M. C.; Sculley, J.; Yu, J.; Jeong, H. K.; Balbuena, P. B.; Zhou, H. C. Exceptional Chemical and Thermal I


Article


(42) Allen, M.; Tildesley, D. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (43) Einstein, A. The Motion of Elements Suspended in Static Liquids as Claimed in the Molecular Kinetic Theory of Heat. Ann. Phys. 1905, 17, 549−560. (44) Darken, L. S. Diffusion, Mobility and Their Interrelation through Free Energy in Binary Metallic Systems. Trans. AIME 1948, 175, 184−201. (45) Theodorou, D. N.; Snurr, R. Q.; Bell, A. T. In Comprehensive Supramolecular Chemistry; Alberti, G., Bein, T., Eds.; Pergamon Press: New York, 1996; Vol. 7, pp 507−548. (46) Skoulidas, A. I.; Sholl, D. S. Transport Diffusivities of CH4, CF4, He, Ne, Ar, Xe, and SF6 in Silicalite from Atomistic Simulations. J. Phys. Chem. B 2002, 106, 5058−5067. (47) Widom, B. Some Topics in the Theory of Fluids. J. Chem. Phys. 1963, 39, 2808−2812. (48) Theodorou, D. N. In Materials Science of Membranes for Gas and Vapor Separation; Yampolskii, Y., Pinnau, I., Freeman, B. D., Eds.; John Wiley & Sons, Inc.: New York, 2006; pp 49−94. (49) Raptis, N. E.; Economou, I. G.; Theodorou, D. N.; Petrou, J.; Petropoulos, J. H. Molecular Dynamics Simulation of Structure and Thermodynamic Properties of Poly(dimethylsilamethylene) and Hydrocarbon Solubility Therein: Toward the Development of Novel Membrane Materials for Hydrocarbon Separation. Macromolecules 2004, 37, 1102−1112. (50) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435−447. (51) Skoulidas, A. I.; Sholl, D. S. Transport Diffusivities of CH4, CF4, He, Ne, Ar, Xe, and SF6 in Silicalite from Atomistic Simulations. J. Phys. Chem. B 2002, 106, 5058−5067. (52) Li, K.; Olson, H. D.; Seidel, J.; Emge, J. T.; Gong, H.; Zeng, H.; Li, J. Zeolitic Imidazolate Frameworks for Kinetic Separation of Propane and Propene. J. Am. Chem. Soc. 2009, 131, 10368−10369.

(21) Zheng, B.; Sant, M.; Demontis, P.; Suffritti, G. B. Force Field for Molecular Dynamics Computations in Flexible ZIF-8 Framework. J. Phys. Chem. C 2012, 116, 933−938. (22) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (23) Parkes, M. V.; Hakan, D.; Teich-McGoldrick, S. L.; Sholl, D. S.; Greathouse, J. A.; Allendorf, M. D. Molecular Dynamics Simulation of Framework Flexibility Effects on Noble Gas Diffusion in HKUST-1 and ZIF-8. Microporous Mesoporous Mater. 2014, 194, 190−199. (24) Chokbunpiam, T.; Chanajaree, R.; Saengsawang, O.; Reimann, S.; Chmelik, C.; Fritzsche, S.; Caro, J.; Remsungnen, T.; Hannongbua, S. The Importance of Lattice Flexibility for the Migration of Ethane in ZIF-8: Molecular Dynamics Simulations. Microporous Mesoporous Mater. 2013, 174, 126−134. (25) Wu, X.; Huang, J.; Cai, W.; Jaroniec, M. Force Field for ZIF-8 Flexible Frameworks: Atomistic Simulation of Adsorption, Diffusion of Pure Gases as CH4, H2, CO2 and N2. RSC Adv. 2014, 4, 16503− 16511. (26) Zheng, B.; Pan, Y.; Lai, Z.; Huang, K. W. Molecular Dynamics Simulations on Gate Opening in ZIF-8: Identification of Factors for Ethane and Propane Separation. Langmuir 2013, 29, 8865−8872. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A. Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision D.01; Gaussian Inc.: Wallingford, CT, 2013. (28) Becke, A. D. Density-functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (29) Lee, C.; Yang, W.; Parr, R. G. Development of the Colle-Salvetti Correlation-energy Formula into a Functional of the Electron-density. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (30) Lin, F.; Wang, R. Systematic Derivation of AMBER Force Field Parameters Applicable to Zinc-Containing Systems. J. Chem. Theory Comput. 2010, 6, 1852−1870. (31) Seminario, J. M. I. Calculation of Intramolecular Force Fields from Second-derivative Tensors. Int. J. Quantum Chem. 1996, 60, 1271−1277. (32) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (33) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (34) Singh, U. C.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129− 145. (35) Besler, B. H.; Merz, K. M., Jr.; Kollman, P. A. Atomic Charges Derived from Semiempirical Methods. J. Comput. Chem. 1990, 11, 431−439. (36) Martin, M.; Siepmann, J. Transferable Potentials for Phase Equilibria. 1. United-atom Description of n-alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (37) Michalis, V. K.; Moultos, O. A.; Tsimpanogiannis, I. N.; Economou, I. G. Molecular Dynamics Simulations of the Diffusion Coefficients of Light n-alkanes in Water Over a Wide Range of Temperature and Pressure. Fluid Phase Equilib. 2016, 407, 236−242. (38) Makrodimitri, Z. A.; Unruh, D. J. M.; Economou, I. G. Molecular Simulation and Macroscopic Modeling of the Diffusion of Hydrogen, Carbon Monoxide and Water in Heavy n-alkane Mixtures. Phys. Chem. Chem. Phys. 2012, 14, 4133−4141. (39) Discovery Studio Modeling Environment, Release 4.5; Dassault Systèmes: San Diego, CA, 2015. (40) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (41) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697. J


Molecular Simulation Studies of the Diffusion of Methane, Ethane

Recommend Documents