Theoretical model estimation of guest diffusion in metal–organic frameworks (MOFs)

Abstract

Characterizing molecule diffusion in nanoporous matrices is critical to understanding the novel chemical and physical properties of metal–organic frameworks (MOFs). In this paper, we developed a theoretical model to rapidly and accurately compute the diffusion rate of guest molecules in a zeolitic imidazolate framework-8 (ZIF-8). The ideal gas or equilibrium solution diffusion model was modified to contain the effect of host framework via introducing the possibility of guests passing through the framework gate. The only input in our model is the energy barrier of guests passing through the MOF's gate. Molecular dynamics (MD) methods were employed to gather the guest density profile, which then was used to deduce the energy barrier values. This produced reliable results that require a simulation time of 5 picoseconds, which is much shorter when using pure MD methods (in the nanosecond scale). Also, we used density functional theory (DFT) methods to obtain the energy profile of guests passing through gates, as this does not require specification of a force field for the MOF degrees of freedom. In the DFT calculation, we only considered one gate of MOFs each time; as this greatly reduced the computational cost. Based on the obtained energy barrier values we computed the diffusion rate of alkane and alcohol in ZIF-8 using our model, which was in good agreement with experimental test results and the calculation values from a standard MD model. Our model shows the advantage of obtaining accurate diffusion rates for guests in MOFs for a lower computational cost and shorter calculation time. Thus, our analytic model calculation is especially attractive for high-throughput computational screening of the dynamic performance of guests in a framework.

---

1. Introduction

Owing to their flexible and tunable pore geometry, metal–organic frameworks (MOFs) are fascinating materials that have applications in mixture separation, gas storage, sensing and catalytic conversion. These applications involve a process whereby the guest molecules diffuse in a spatial confinement on a length-scale below 1 nm and in a full flexible framework with breathing and gate-opening effects.¹ The diffusion rate influences the exchange phenomena between the host materials and the surroundings and influences the performance of MOFs' in their various applications.^2–4 Therefore, investigating guest diffusion within a MOF framework is an in demand research focus.

Guest diffusivities can be experimentally determined by recording overall uptake and release stimulated by a step change of external pressure. Microscopically, pulsed-field-gradient NMR and related microimaging techniques were developed to capture diffusion paths and transient concentration profiles, followed by estimating mass transfer in nanoporous materials.^3,5,6 Compared with expensive and time-consuming experimental testing, these theoretical methods that involve computing the guest diffusion rate in MOFs are more economical and timely. Molecular dynamics (MD), based on classical mechanics, is the most widely applied tool for the production of trajectories of guests in the phase space and obtainment of measurable macroscopic properties like diffusion.¹ Although many accurate MD simulations have been reported the results mainly depend on the assumption of inter- and intra-molecular potentials.^7–9 A reasonable semi-empirical force field for MOF systems should represent both the stability and the flexibility of the framework, which in principle increases the difficulty of building a force field.^7,9–13 Density functional theory (DFT) methods based on quantum mechanics can be used to fundamentally understand the structure, properties and reactivity of MOFs.¹⁴ In particular, time-dependent DFT (TD-DFT) methods have the potential to investigate guest diffusion in MOFs. However, the use of a full quantum mechanical treatment of guest diffusion in the complete flexible MOFs system, comprising nuclei and electrons, is far from feasible. This is because the transport phenomena in MOFs span a geometrical range from nanometres (intra-crystalline) to hundreds of micrometres (inter-crystalline), using time scales ranging from femtoseconds (chemical reactions) to microseconds (long range diffusion), which is well beyond the ability of current clusters performing DFT calculations. In order to overcome the limitations imposed by both the time and geometrical scales in the typical theoretical method, an approach combining multiple calculation methods has been developed.^15,16 Ab initio molecular dynamics (AIMD) were used to gather the flexible framework structures, followed by employing these structures to describe guest diffusion with a classical molecular dynamics model. Although this model has advantages in that it accounts for framework flexibility without requiring a force field, the AIMD calculation is time-consuming making it challenging when dealing with larger models (thousands of atoms).

In this work, the conventional diffusion models in ideal gas or equilibrium fluid are modified to suit computing of diffusion rates of guests in periodic media, i.e. MOFs. We show that the diffusivity is mainly influenced by the guest energy barriers passing through the host gate. Short time MD calculations were used to create the guest density profile, followed by deduction of the energy barrier and diffusion rate. Another strategy to compute the energy barrier in this study is to record the energy profile along the preferential path of guests basing DFT calculation results. The model was validated via computing the diffusion rate of alkane and alcohol guests inside a zeolitic imidazolate framework-8 (ZIF-8), a standard type of MOF with exceptional stability. The computed diffusion values agree well with the experimental results.

2. Model development

2.1. Diffusion model

Molecule movement can be described by two conventional models:¹⁷ the ideal gas model (eqn (1)) and the solution model (eqn (2)).where J is the diffusion flux and λ is the free path as the distance travelled or jumped by a particle before its direction of motion is randomized through collisions; v is the average velocity during one free path; N is the concentration in dimensions of diffusion; and x is the position.where P is the probability of the net distance walked being nλ; τ is the time for one whole jump; and t is the total diffusion time.

When comparing the above two models to Fick's laws of diffusion, the diffusion rate (D) can be described as below:

In eqn (3), τ is the travelling or the jumping time, which ignores the residence time at the terminal point. This assumption is correct for ideal gas and free fluid. However, our model involves periodic medium, i.e. the MOF's gates locate in space uniformly, as shown in Fig. 1. In the MOF framework, λ is the distance between two adjacent gates and τ should include both the travelling time inside the cavity (τ₀) and gate hopping (τ₁). The τ₀ is mainly influenced by the adsorption ability of the framework, while the τ₁ involves the complexity of the host–guest system, including the guest–gate interactions, the gate flexibility (gate-opening) and guest geometry. Although quantification of τ₁ is difficult, its effect on slowing the diffusion rate of guest molecules is clear. This dragging effect can be taken into account via introducing the factor, p, in eqn (4), which refers to the possibility of guests passing through the gate. As such, when p = 0 there is zero diffusion no matter how fast the guest is initially. When no gate barrier is in the matrix (p = 1), the guest will diffuse according to classical collision theory.

where v₀ = λ/τ₀.

Fig. 1 One dimensional model to describe the guest molecule diffusion in porous matrix. This framework consists of periodic gates where the thickness of the gate can be ignored. The diffusion path of the guest molecules is described by arrows.

When the adsorption effect is not very strong, v₀ can be obtained from the temperature, according to the ideal gas equation (eqn (5)).

where k is the Boltzmann constant.

In order to define the possibility factor (p) in eqn (4), the microscopic process for guest diffusion was analysed. In the porous model, the guest diffusion process consists of two parts: intra- and inter-cage movement. The former process can be simplified to the ideal gas model (eqn (5)). For the latter, guest molecules will first approach the gate of the matrix before passing through it. Usually, the long-range attractive interactions between the guest and the host can drive the approaching process.⁸ Once distances between the guest and host are smaller than their equilibrium distance, atomic repulsive interactions become dominant. The kinetic energy of the guests will contribute to compensate for the energy consumption due to increasingly repulsive atomic interactions, which corresponds to an energy barrier. This guarantees the total energy conservation in the host–guest system.

The hopping process of guest molecules in a porous matrix is more likely to be a chemical reaction process, in which molecules are supposed to react if they collide with a relative kinetic energy that exceeds the activation energy. This assumption is also used in other work,¹⁵ in which the hopping rate of CH₄ in ZIF-8 was calculated through using transition state theory (TST). Here, we will employ the Arrhenius equation, a formula for the temperature dependence of reaction rates, to define the possibility (p) in eqn (4):

k = Ae^(−E/RT)

p₀ = k/A = e^(−E/RT) (6)

where k is the number of collisions that result in a successful crossing per second; A is the total number of collisions (leading to a crossing or not) per second, p₀ is the probability that any given collision will result in a crossing; E is the energy barrier; and R is the universal gas constant. It is noted that the boundary in eqn (6) is shown as: E = 0, p₀ = 1; E = +∞, p₀ = 0.

2.2. Energy barrier calculation

Energy barrier values were obtained through using two calculation strategies. One employed the Boltzmann-weighted free energy equation,¹⁸ shown as below:

E_B = −k_BT ln(ρ) (7)

where E_B and ρ are the Boltzmann free energy profile and sorbate probability respectively.

In eqn (7), the only input is the density profile of the guest molecules, which can be obtained by classical MD calculations, whose details can be found in our earlier work.⁸ Compared with the diffusion rate calculation in MD, the density profile requires much shorter calculation times. Our test indicates that the density distribution is reliable even for a 5 ps calculation time, while reliable diffusion rates require calculation times of at least 10 ns. However, the accuracy of the density profile calculation mainly depends on the semi-empirical force field. The different force fields for the atomic framework can lead to significantly different calculated energy barrier values.

Another strategy is to use the DFT calculation for the energy barrier, which takes into account the system comprising of nuclei and electrons. The position of the guest molecules along the preferential path was manually adjusted to where the block from the geometrical contact between host and guests was weak. The system's energy was recorded along the path to obtain the energy barrier for the guest molecules passing through the gate. DFT calculations were performed using the Dmol3 code.^19,20 The exchange correlation energy was described by the PW91 functional within the generalized gradient approximation (GGA-PW91).²¹ DFT semi-core pseudo potentials were used for all atoms, as were the double-numeric basis with polarization functions (DNP).

3. Results and discussion

Fig. 2(a) shows a homogeneous density profile of the guest molecules, which are denoted as the centre of mass. In order to overcome relatively poor statistics problem, we rescaled all the guest molecule coordinates to the primitive unit cell to form the sorbate trajectory. Then, discretizing the simulation box in cubic bins (each of side 0.34 Å), we computed the probability (P₁) of finding a sorbate molecule in the given bin, satisfying the condition where V is the volume of the cubic unit cell. The density profile of sorbate described a strongly structured guest distribution function with the presence of well-defined adsorption sites. The adsorption sites for guest molecules were close to the centre of the six-membered ring window and were rarely found inside the window of ZIF-8 (Fig. 2a). Furthermore, the density distribution was projected as a diagonal line (red line in Fig. 2a) in the cubic cell, which perpendicularly passes through the centre of the gate. The density profile of the guest molecules was used in eqn (7) to obtain the free energy profile with acceptable statistics, shown in Fig. 2b. From this figure, it is easy to obtain the energy barrier for guests passing through the gate via subtracting the lowest energy (adsorption energy) from the highest energy.

Fig. 2 (a) Isoprobability contours for guest molecules inside one unit cell of ZIF-8, projected onto the x–y-plane, with darker areas denoting higher probability. (b) The free energy profile of host–guest system calculated by eqn (7) along the path shown as in (a).

The DFT calculation was also employed to compute the energy barrier of the guest passing through the gate. The calculation model, which is a triclinic primitive cell of ZIF-8, consisted of one six-member ring gate (Fig. 3a), with calculation details described in our previous report.⁸ The red line is the initial path, which is the same as for guest diffusion (Fig. 2a). This path is divided into many points. At every point, one skeleton atom of the guest molecule is fixed and then geometry optimization and single point calculations for the whole system were performed to collect the energy value. For example, propane consisting of three carbons as the skeleton atom, will result in three curves of energy distribution (Fig. 2b).

Fig. 3 The model for DFT calculation of the energy barrier of guests in ZIF-8. (a) The preferential path for the guest molecules passing through one gate of ZIF-8. (b) The calculated energy of the host–guest system along the path as shown in (a), with different fixed skeleton atoms of guest molecules.

For every energy curve (Fig. 2b), two low energy points for the adsorption sites located each side of the gate and one high energy point for the gate-opening can be found. All of the energy points were used to find the lowest energy point at the side and centre respectively, corresponding to real adsorption and gate-opening geometry. The adsorption sites were refined through another time of geometry optimization without fixing any atoms. For the gate-opening point, adjustments were made to the position of the guest molecules to confirm if it is the lowest energy for gate-opening. Once we obtained reliable energy values, the energy barrier was described as the difference between the adsorption energy and the gate-opening energy.

The energy barrier values for alkane guest molecules passing through ZIF-8 gate were calculated using the MD and DFT method (Table 1), where the MD method generates lower energy values than the DFT calculation. This can be attributed to the temperature effect in MD (300 K), compared to 0 K for the DFT calculation.

Table 1 The energy barrier values obtained from MD and DFT calculations for guest molecules passing through one six-member ring gate in ZIF-8 and diffusion rate comparison against standard models and experimental tests

			CH4	C2H6	C3H8
Barrier (kcal mol^−1)	MD		4.2	5.5	7.8
	DFT		5.0	5.9	9.2
Diffusion rate (m^2 s^−1)	Current model	MD	3.5 × 10^−10	2.5 × 10^−11	4.2 × 10^−13
		DFT	8.4 × 10^−11	1.3 × 10^−11	4.2 × 10^−14
	Standard MD		3.4 × 10^−10,^8 1.9 × 10^−10(ref. 10)	3.4 × 10^−11,^8 6.0 × 10^−11(ref. 13)	3.0 × 10^−13(ref. 8)
	Experiment		6.8 × 10^−11,^22 (0.9–1.4) × 10^−10(ref. 6)	2.2 × 10^−11,^22 (0.6–2.1) × 10^−11(ref. 11)	1.3 × 10^−13,^22 2.0 × 10^−14(ref. 23)

---

Furthermore, energy barrier values were used to compute the diffusion rate of guest molecules in ZIF-8 via employing eqn (4). The values computed using our model are shown in Table 1. The diffusion rates computed from our model show good agreement with those calculated using the standard MD procedure, where a simulation was required to calculate the rate by mean squared displacement using the Einstein relation (eqn (8)). The standard MD model simulates the 3-dimensional diffusion of the guest molecules with different loading, which is close to a realistic system but requires a long simulation time (50 ns). Comparatively, our analytic model uses a much shorter simulation time (5 ps) yet obtains diffusion rates that agree well with experimental diffusion rates. The energy barrier values from DFT calculations were also used to obtain comparable diffusion rates, except for C₃H₈. The results indicate a clear geometry effect of the zigzag propane molecule when the temperature was set to zero.²⁴

where d₀ is the dimensionality of the system and r_i(t) is the position vector of the sorbate molecule i at time t. D_S is averaged over all N sorbed molecules over multiple time origins t₀ (as symbolized by the brackets).

The diffusion rates (Table 1) were used to compute the separation factors (SF) of alkane molecules by ZIF-8 (Table 2). The experimental SF value for C₂H₆/C₃H₈ was about ∼50 times larger than that for CH₄/C₂H₆. Only the model using the energy barrier from the DFT calculation is in good agreement with experimental results,²² with SF values of 318.8 and 6.3 for C₂H₆/C₃H₈ and CH₄/C₂H₆ respectively. It is worthy of note that temperature was not taken into account in the DFT calculation. Whereas the SF values computed using the energy barrier from MD that consider temperature were 59.5 and 14.2 for C₂H₆/C₃H₈ and CH₄/C₂H₆, respectively. Therefore, the sharp separation between C₂H₆/C₃H₈ and CH₄/C₂H₆, was reproducible using our model, although the absolute values are different between the experiment test and theoretical simulation. This difference can be attributed to the assumption of our model, where a one dimensional framework is considered and the intra-cavity energy barrier is ignored.

Table 2 Comparison of the diffusion SFs of alkane by ZIF-8 from two different calculation models and experimental tests. Diffusion values from the same model were used to calculate the separation factors of alkane

	Separation factor
	Current model		Standard MD^8	Experiment^22
	MD	DFT
CH4/C2H6	14.2	6.3	10.6	3.1
C2H6/C3H8	59.5	318.8	106.7	169.2

---

Our model was further validated through computing the diffusion rate of alcohol molecules within ZIF-8 (Table 3). Compared with the alkane guests, alcohol molecules contain a nonpolar hydroxyl group where the electrostatic host–guest interactions become very important. However, it is difficult to find an appropriate semi-empirical force field to perform MD calculations for the alcohol–MOF system. The DFT method, without requiring a semi-empirical force field, was employed to calculate the energy barrier of alcohol passing through the ZIF-8 gate.

Table 3 The calculated energy barrier values and diffusion rates for alcohol molecules in ZIF-8

	Methanol CH3OH	Ethanol C2H5OH	Propanol C3H7OH
Energy barrier (kcal mol^−1)	5.1	6.3	10.8
Diffusion rate (m^2 s^−1)	7.1 × 10^−11	9.3 × 10^−12	4.7 × 10^−15

---

The calculated energy of the host–guest system during the propanol molecule moving along the preferential path indicates that the adsorption energy on each side of the gate is different (Fig. 4). In fact, the left and right adsorption sites correspond to the methyl terminal and the hydroxyl terminal of the alcohol molecule respectively. The hydroxyl group is more electron withdrawing than the methyl group, therefore resulting in stronger adsorption when directionality is fixed in favour of the hydroxyl group entering the gate first.

Fig. 4 The calculated energy profile for the host–guest system of a propanol molecule passing through the ZIF-8 gate.

From Fig. 4 it is clear that the energy barrier for a guest crossing the gate actually depends on which terminal end of the alcohol molecule approaches the gate first. The energy difference between the two cases is about 2–3 kcal mol⁻¹, which is so small that the priority of the hydroxyl terminal can be ignored and the probability for the hydroxyl and methyl terminal approaching the gate first is considered to be equal. The two adsorption energy values were employed to compute the energy barrier values and their arithmetical average was considered as the final energy barrier (Table 3). By comparing the energy barriers of alkane (Table 1) and alcohol (Table 3), the same rising trend between barrier energy and increasing the length of the carbon chain was obtained.

Furthermore, the energy barrier of alcohol guests was employed to compute the diffusion rate of alcohol in ZIF-8 using our model (eqn (4)). The computed diffusion rate (Table 3) agrees well with the experimental values and the diffusion of the alcohols sharply decreases between ethanol to propanol, as was similarly found for alkane.

A sharp decrease in diffusion rate from C2 to C3 was observed with a larger inter-cavity energy barrier difference (3–4 kcal mol⁻¹), compared to the ∼1 kcal mol⁻¹ difference between C1 and C2 (Tables 1 and 3). Although sharp separation was obtained, the 3–4 kcal mol⁻¹ energy difference is still relatively small. In order to isolate contributions from the barrier difference, eqn (9) was used to compute the SF (Fig. 5). Three points, 3, 4 and 5 kcal mol⁻¹, correspond to the SFs of 147, 779 and 4117 respectively. This describes an exponential increase in SF with increasing energy differences, when considering guest geometry, framework flexibility, guest–host interactions and temperature effects. The current model provides a simple route to analysing the guest's separation by MOFs, which can guide both experimental testing and practical applications.

Fig. 5 The separation factors vs. inter-cavity energy barrier of guests in ZIF-8.

Finally, we present a rough estimate of the time efficiency attainable using our analytical model to compute the guest diffusion in ZIF-8. The most time-consuming aspect of the calculation is in computing the energy barrier. Using the DFT method, 3 energy points (two adsorption points and one transition state point) are required to obtain the barrier. Every energy point of CH₄–ZIF system in Fig. 2b takes about 5 hours on 8 core clusters (Inter Xeon X5570). Using the MD method, a 5 ps MD simulation of a 2 × 2 × 2 ZIF-8 to obtain reliable distribution of CH₄ guest takes about 0.1 hours. On the other hand, frequently-used MD methods that obtain the diffusion values from the slope of the calculated mean square displacement, require a 50 ns simulation of the same system and takes about 125 hours. The same trend holds also for the remaining guest–ZIF systems.

Conclusions

In summary, a simple model to compute the diffusion rate of guest molecules in MOFs was developed. In this model, the diffusion rate mainly depends on the energy barrier of guests passing through the MOFs' gate. Two strategies, using either a density profile from MD or energy profile from DFT, were used to obtain the energy barrier values. The proposed model, which incorporates the energy profile from DFT and MD, was applied to compute the diffusion rate of alkane and alcohol guests in ZIF-8, agreeing well with the experimental test results. Using our model, a sharp C2/C3 separation in both alkane and alcohol guests by ZIF-8 was observed. This simple model can generate accurate results using short simulation times and small framework sizes when compared to pure MD or DFT models. Future research will aim to develop three-dimensional models for computing the diffusion rates of guests in other MOF systems to understand and guide experimental work.

Acknowledgements

This work was supported by the Natural Science Foundation of China under grant 51372197, the Key Innovation Team of Shaanxi Province (2014KCT-04), the Major International Joint Research Program of Shaanxi Province (2012KW-10), and the Doctoral Starting up Foundation of Xi'an University of Science and Technology (2015QDJ023).

Notes and references

1. P. Demontis and G. B. Suffritti, Microporous Mesoporous Mater., 2009, 125, 160–168.
2. M. Tsotsalas, P. Hejcik, K. Sumida, Z. Kalay, S. Furukawa and S. Kitagawa, J. Am. Chem. Soc., 2013, 135, 4608–4611.
3. J. Karger, T. Binder, C. Chmelik, F. Hibbe, H. Krautscheid, R. Krishna and J. Weitkamp, Nat. Mater., 2014, 13, 333–343.
4. T. Titze, C. Chmelik, J. Kullmann, L. Prager, E. Miersemann, R. Gläser, D. Enke, J. Weitkamp and J. Kärger, Angew. Chem., Int. Ed., 2015, 54, 5060–5064.
5. C. Chmelik, D. Freude, H. Bux and J. Haase, Microporous Mesoporous Mater., 2012, 147, 135–141.
6. A.-K. Pusch, T. Splith, L. Moschkowitz, S. Karmakar, R. Biniwale, M. Sant, G. Suffritti, P. Demontis, J. Cravillon, E. Pantatosaki and F. Stallmach, Adsorption, 2012, 18, 359–366.
7. X. Wu, J. Huang, W. Cai and M. Jaroniec, RSC Adv., 2014, 4, 16503–16511.
8. B. Zheng, Y. Pan, Z. Lai and K.-W. Huang, Langmuir, 2013, 29, 8865–8872.
9. B. Zheng, M. Sant, P. Demontis and G. B. Suffritti, J. Phys. Chem. C, 2012, 116, 933–938.
10. L. Hertäg, H. Bux, J. Caro, C. Chmelik, T. Remsungnen, M. Knauth and S. Fritzsche, J. Membr. Sci., 2011, 377, 36–41.
11. M. V. Parkes, H. Demir, S. L. Teich-McGoldrick, D. S. Sholl, J. A. Greathouse and M. D. Allendorf, Microporous Mesoporous Mater., 2014, 194, 190–199.
12. L. Sarkisov, R. L. Martin, M. Haranczyk and B. Smit, J. Am. Chem. Soc., 2014, 136, 2228–2231.
13. T. Chokbunpiam, R. Chanajaree, O. Saengsawang, S. Reimann, C. Chmelik, S. Fritzsche, J. Caro, T. Remsungnen and S. Hannongbua, Microporous Mesoporous Mater., 2013, 174, 126–134.
14. M. K. Rana, F. G. Pazzona, G. B. Suffritti, P. Demontis and M. Masia, J. Chem. Theory Comput., 2011, 7, 1575–1582.
15. E. Haldoupis, T. Watanabe, S. Nair and D. S. Sholl, ChemPhysChem, 2012, 13, 3449–3452.
16. M. Sant, G. K. Papadopoulos and D. N. Theodorou, J. Chem. Phys., 2010, 132, 134108.
17. P. Atkins and J. De Paula, Physical Chemistry, Oxford University Press, Oxford, 2006.
18. M. Sant, G. K. Papadopoulos and D. N. Theodorou, J. Chem. Phys., 2008, 128, 024504.
19. B. Delley, J. Chem. Phys., 1990, 92, 508–517.
20. B. Delley, J. Chem. Phys., 2000, 113, 7756–7764.
21. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter, 1992, 45, 13244–13249.
22. Y. Pan and Z. Lai, Chem. Commun., 2011, 47, 10275–10277.
23. C. Zhang, R. P. Lively, K. Zhang, J. R. Johnson, O. Karvan and W. J. Koros, J. Phys. Chem. Lett., 2012, 3, 2130–2134.
24. Y. Ma and P. B. Balbuena, Chem. Phys. Lett., 2012, 552, 136–140.

---

Footnote

† The authors declare no competing financial interest.

---