J. Phys. Chem. C 2008, 112, 2805-2811
2805
ARTICLES High-Accuracy Estimation of ‘Slow’ Molecular Diffusion Rates in Zeolite Nanopores, Based on Free Energy Calculations at an Ultrahigh Temperature Ryo Nagumo,*,† Hiromitsu Takaba,‡ and Shin-ichi Nakao† Department of Chemical System Engineering, The UniVersity of Tokyo, 7-3-1 Hongo Bunkyo-ku, Tokyo 113-8656, Japan, and Department of Applied Chemistry, Tohoku UniVersity, 6-6-11-1302 Aoba Aramaki Aoba-ku, Sendai 980-8579, Japan ReceiVed: April 27, 2007; In Final Form: August 13, 2007
We present a methodology to estimate ‘slow’ molecular diffusivities through nanopores, which have been difficult to calculate from conventional molecular dynamics (MD), combining statistical-mechanical theory with MD simulations. In this method, the hopping rate of a guest molecule is calculated by transition state theory (TST) and used to estimate self-diffusivities. To obtain hopping rates from TST, it is essential to obtain free energy profiles along the reaction coordinate. However, the free energy calculation can involve considerable error because of insufficient statistical samplings, especially simulating ‘slow’ molecular diffusion. To solve this problem, Helmholtz free energy profiles at much lower temperatures are predicted from ultrahigh temperature MD, based on classical statistical-mechanical theory. Using these predicted profiles, selfdiffusivities at lower temperatures were estimated. This methodology was applied to estimate molecular diffusivities of CH4 through an LTA-type zeolite. Our predicted diffusivities showed excellent agreement with those obtained from conventional MD. Our methodology is a promising approach to the systematic prediction of ‘slow’ molecular diffusivity in various types of nanoporous material.
1. Introduction Zeolites are among the most significant nanoporous materials, whose unique characteristics, such as molecular sieving, preferential adsorption, and shape-selective catalysis, have given them a remarkable importance in industrial separation and reaction processes.1,2 To promote more applications in various industrial fields, an essential task is the estimation of the performance and properties of zeolites, such as catalytic activity or life, and permeate flux or separation factor. In particular, molecular self-diffusivity of guest molecules is one of the most important factors for such properties. Molecular dynamics (MD) simulation has been a powerful tool to estimate selfdiffusivities.3-10 In our group, molecular diffusivities of inorganic guest molecules11 and linear alkanes12,13 within zeolites have been calculated. Most of these previous studies have dealt with the guest-host systems, whose molecular diffusivities are in the order of more than 10-10 m2/s. However, industrially important guest-host systems, like the separation process of aromatic isomer mixtures through an MFItype zeolite, generally show a ‘slow’ diffusion regime, where molecular diffusivity is much less than 10-12 m2/s. In this regime, classical MD is not practical, because the simulations would take too long (a few or more years) on present computers. Therefore, it is essential to calculate with high accuracy the ‘slow’ diffusion rate of guest molecules. Several methodologies have been developed to simulate ‘slow’ diffusion phenomena, * Corresponding author. E-mail:
[email protected]. † The University of Tokyo. ‡ Tohoku University.
as shown in a previous review by Voter et al.14 One of the representative methods is temperature-accelerated dynamics (TAD), proposed by Sorensen et al.15 TAD has already been applied to some dynamic phenomena, such as diffusion16 or crystal growth17 on a metal surface. Further applications of TAD will certainly be expected in the near future, but TAD involves some assumptions, such as the Arrhenius temperature dependence of the hopping rates. Therefore, some improvements are still necessary for the highly accurate simulation of diffusion phenomena. On the other hand, in the zeolitic host-guest systems, one of the promising methodologies to calculate the ‘slow’ diffusion rates is transition state theory (TST).18-20 For example, Mosell et al. 21-24 calculated the self-diffusivities of guest molecules inside the zeolite NaY by TST. Smit et al.25 also used TST to calculate the self-diffusivities of branched alkanes in MFI-type zeolite. Beerdsen and co-workers26-28 recently used TST to calculate the self-diffusivities of n-alkanes within an LTA-type zeolite, incorporating the dynamical correction by TST. Their method is remarkable in that the concentration dependence of self-diffusivities can be estimated, not only the value at infinite dilution. As shown in these studies, TST is certainly powerful in calculating molecular diffusivity within nanoporous materials, but it still includes some problems in achieving systematic estimations. In TST, a free energy calculation is essential for accurate estimation of diffusivities. However, especially for ‘slow’ diffusion phenomena, it is difficult to sample enough data for the calculation. Therefore, an alternative approach or
10.1021/jp073250b CCC: $40.75 © 2008 American Chemical Society Published on Web 02/05/2008
2806 J. Phys. Chem. C, Vol. 112, No. 8, 2008
Nagumo et al.
Figure 2. Schematic of the algorithm to calculate the transmission coefficient. This is computed as the probability that a guest molecule, originating in the reactant region, R, starting at the transition state (TS) at time zero with a positive velocity, Vx, reaches the product region, P, directly without recrossing the TS. Figure 1. LTA-type zeolite structure as a network of three-dimensional cubic lattices of R-cages. Lattice distance between the adjacent cage centers is 1.2305 nm.
some improvements are needed to realize higher accuracy and a more systematic estimation. To perform a high-accuracy calculation, we propose a methodology to estimate ‘slow’ diffusivity inside nanoporous materials using classical MD, based on statistical-mechanical theory. Our approach is applicable to molecular diffusion inside a nanopore that shows an ‘activated’ temperature dependence, where the self-diffusivity increases drastically with an increase in temperature. Ultrahigh temperature heating of the guesthost system leads to the applicability of conventional MD for free energy calculations. That is, an MD calculation at a higher temperature increases the accuracy of the data sampling for the free energy calculation. In our method, from MD simulations at an extremely high temperature, free energy profiles at lower temperatures can be predicted. Then the ‘slow’ diffusivity at a lower temperature is estimated by TST. In this study, we present the details of our methodology and the data required to estimate self-diffusivities of CH4 through an LTA-type zeolite, such as transmission coefficients and free energy profiles, at a concentration of two molecules per cage. We also show thermodynamic quantities, such as potential energy, internal energy, and entropy. These data enable us to investigate the molecular diffusion mechanism from a thermodynamic viewpoint. 2. Methodology 2.1. Self-Diffusivity. The LTA-type structure consists of nanopores whose pore diameter is ca. 0.41 nm. It is considered as a network of three-dimensional cubic lattices of R-cages, whose diameter is about 1.12 nm, as shown in Figure 1. In this coarse-grained model, self-diffusivity, Ds, can be calculated from
Ds ) kλ2
(1)
where λ is the lattice distance between the adjacent LTA cage centers (ca. 1.23 nm), and k is the hopping rate from one cage center to a neighboring cage via a window. TST formulates k as
k ) κkTST
(2)
where κ is the transmission coefficient and kTST is the TST rate constant. In the next sections, the detailed method to calculate κ and kTST is described.
2.2. Transmission Coefficient. The transmission coefficient, κ, is defined as the probability that a molecule, originating in the reactant region and starting at the transition state at time zero with a positive velocity toward the product region, reaches directly the product region without recrossing the transition state. The detailed algorithms for calculating κ have been described in a previous study.29 They investigated the accuracy and efficiency of several algorithms, estimating the transmission coefficient for the diffusion of the potassium ion through an ion channel. In our study, referring to their investigation, the Anderson algorithm30 was used to calculate κ, although some improvements to the algorithm were made, as shown in Figure 2. In this algorithm, κ is estimated from independent and unconstrained multiple MD simulations. The initial configurations of molecules for each MD simulation are randomly generated, but the initial position of the probe particle is set at the top of the barrier, i.e.,
x(0) ) x*
(3)
where x(t) is the reaction coordinate of the probe at time t. The initial molecular velocities should be determined randomly on the Maxwell-Boltzmann distribution. After carrying out n unconstrained MD simulations, κ is calculated from
κ)
nTSfproduct n
(4)
where the numerator nTSfproduct is the number of trajectories originating in the reactant region, starting at the transition state, and proceeding directly to the product region. It has to be ensured that the n initial configurations are completely uncorrelated. Furthermore, there are some points to be noted for calculating κ, as described in the following paragraphs. First, repeating multiple MD simulations, some criterion is needed for judging the arrival of the probe molecule in the reactant or product region. In this study, the distance between the initial position of a probe molecule (i.e., the transition state) and the other region, dmax, was introduced. However, the movement distance of the probe can remain less than dmax, even if a long time has elapsed. In such a case, it is natural to consider that, if the elapsed time t becomes larger than the given time interval, relaxation from the initial configuration should be achieved. In this study, when the elapsed time reaches tmax, it is determined whether the probe is either in the reactant side or the product side in each MD simulation. Second, in the calculation of the numerator of eq 4, we need to count only the number of probes starting at the transition
Slow Molecular Diffusion Rates in Zeolite Nanopores
J. Phys. Chem. C, Vol. 112, No. 8, 2008 2807
state and reaching directly the product region, without recrossing the transition state. If the probe molecule recrosses the transition state, we do not need to count it and the MD run can be terminated. Third, it should be ensured that the probe molecules, starting at the transition state at t ) 0, originate in the reactant region. Therefore, following the approach of Frenkel et al.,31 we performed backward simulations, reversing the sign of the initial velocities of all the molecules inside the simulation model but starting from the same initial configuration. According to the previous study by Ruiz-Montero et al.,32 the statistical error in κ, σ, is estimated by
σ 1 ∼ κ κxn
(5)
It is noted that when κ is much less than 1, a larger number of MD runs is needed. For example, if κ is about 0.1, 104 runs have to be taken to obtain 10% accuracy. 2.3. TST Rate Constant. kTST is expressed in TST as
kTST )
x
1 2πβm
exp{-βA(T, x*)}
∫cage
(6)
exp{-βA(T, x)}dx
where β ) 1/kBT, kB is the Boltzmann constant, T is the temperature, m is the particle mass involved in the reaction coordinate, x, and A(T, x) is the Helmholtz free energy as a function of x. Here, the position of the top of the free energy profile (i.e., the transition state) is denoted as x*. From eq 6, kTST is estimated by calculating the free energy profile. We used the window sampling method33 to calculate A(T, x). In this method, A(T, x) can be easily calculated from
exp{-βA(T, x)} ∝ n(T, x)
(7)
where n(T, x) is the particle density distribution, obtained by counting the number of particles, N(x), existing between the two positions x and x + ∆x during the MD simulations. Then n(T, x) is calculated from
n(T, x) )
1 ∆x
N(x)
∑ N(x)
1 x2πβm 1 - n(T, x*)∆x n(T, x*)
(11)
where QN(V, T) is the canonical partition function. QN(V, T) is expressed as
1 h3NN!
QN(V, T) )
∫ ∫ exp[-βHN(rN, pN)]drNdpN
(12)
where h is Planck’s constant, rN is the set of position vectors of N particles, pN is the set of momentum vectors of N particles, and HN(rN, pN) is the classical Hamiltonian. In classical statistical mechanics, the integral term on pN of eq 12 can be treated as Gaussian for integration. Therefore, QN(V, T) is given by
QN(V, T) )
{(x ) } βh2 2πm
3N
-1
N!
∫ exp[-βVN(rN)]drN
(13)
where VN(rN) is the potential energy of an ensemble. From eqs 11 and 13, the temperature dependence of the Helmholtz free energy is formulated as
A(T, x) ) RT ln
T A(Tref, x) Tref
( ) T Tref
3N/2
〈 [(
- RT ln exp -
) ]〉
1 1 V RT RTref N
(14)
T)Tref
where R is the gas constant and 〈...〉T)Tref is the ensemble average obtained from the MD at Tref. Performing MD simulations at one temperature, TST rates at different temperatures can be calculated using eqs 6 and 14. 2.4. Thermodynamic Quantities. Using the results from the window sampling method, the profiles of thermodynamic quantities along the reaction coordinate can be obtained. First, the free energy profile A(T, x) is calculated from eq 7, given by
-βA(T, x) ) ln{n(T, x)} + C
(15)
where C is a constant. The internal energy, U(T, x), defined as the sum of the kinetic energy, K(T, x), and the potential energy, V(T, x), is expressed as
U(T, x) ) K(T, x) + V(T, x)
(9)
From eqs 7 and 8, eq 6 is expressed by
kTST ≈
-βA(T, x) ) ln{QN(V, T)}
(8)
Here it is noted that n(T, x) should be normalized by
∫cage n(T, x)dx ) 1
When an NVT ensemble (number of atoms, volume, and temperature in the simulation model are kept constant) is applied, the Helmholtz free energy is formulated as
(10)
Therefore, kTST can be calculated from eq 10. However, in performing these calculations, it is generally difficult to estimate the particle density distribution, i.e., Helmholtz free energy, with a statistically acceptable accuracy. This problem becomes quite serious, particularly for the ‘slow’ diffusion regime. One of the alternative approaches is to increase the temperature of the whole system. This leads to an acceleration of molecular diffusion rates. Therefore, sufficient sampling can be taken from conventional MD simulations. Consequently, A(T, x) at a high enough temperature is easily calculated. In this study, the prediction of free energy profiles at lower temperatures was performed from MD simulations at a higher temperature, by the methodology described below.
(16)
K(T, x) and V(T, x) can be easily obtained from MD simulations. Thus U(T, x) is estimated from eq 16. In addition, the entropy, S(T, x), can also be calculated from
-TS(T, x) ) A(T, x) - U(T, x)
(17)
In this study, these thermodynamic quantities were calculated using eqs 15-17. 3. Simulation We studied the diffusion phenomena of CH4 in a threedimensional LTA-type zeolite. In all our simulations, an LTAtype zeolite was modeled as 2 × 2 × 2 unit cells, where a periodic boundary condition was applied in each direction. The initial coordination of guest molecules was decided from grand canonical ensemble Monte Carlo (GCMC) simulations. The crystallographic unit cell consists of eight R-cages interconnected by windows, whose diameters are about 0.41 nm. To
2808 J. Phys. Chem. C, Vol. 112, No. 8, 2008
Nagumo et al.
Figure 3. Arrhenius temperature dependence of self-diffusivities of CH4 in an LTA-type zeolite calculated from MD simulations at each temperature.
introduce a set of potential functions and interaction parameters, we referred to the previous studies.28 All the guest molecules are considered spherical. It is assumed that guest-host interactions are represented by dispersive forces derived mainly from zeolitic oxygen atoms. Guest-guest and guest-host nonbonded interactions are described with a 12-6 Lennard-Jones (LJ) potential. The LJ potential is truncated at 12.0 Å. The zeolite framework is considered rigid, as it was in our previous reports.11-13 The initial velocity of each molecule was set corresponding to the Boltzmann distribution at the desired temperature. The adsorbate concentration in our MD calculation was set at two molecules per R-cage. The temperatures in the MD simulations were set at 1800 K, 1600 K, 1400 K, 1200, 1000, and 800 K. In the free energy calculations, MD simulations were performed for 1-10 ns in real time, where the unit cell length was sectioned into 240 bins along the reaction coordinate to set the width, ∆x. From these MD simulations, self-diffusivities of each guest species were also calculated using the mean square displacement from Einstein’s equation, as shown in the next section. An NVT ensemble was used with the temperature controlling of Berendsen,34 the time step for each MD step being set at 0.5 fs. On the other hand, in calculating transmission coefficients, an NVE ensemble (number of atoms, volume, and total energy in the simulation model are kept constant) was applied, and the time step for each MD step was set at 0.2 fs. Initial configurations were generated from GCMC simulations, with 300 000 steps, followed by 300 000 steps of NVT ensemble Monte Carlo simulations. The distance to judge whether a probe molecule reaches the product or reactant region, dmax, was set at 0.3 nm and the time interval, tmax, was 10 ps, referring to the previous studies.35 4. Results 4.1. Self-Diffusivity from Conventional MD. The selfdiffusivity is defined using the Einstein expression
Ds ) lim tf∞
1 2 〈|b(t) r - b(0)| r 〉 6t
(18)
where 〈...〉 denotes an ensemble average and b(t) r is the position of the molecule at time t. Using eq 18, we computed the selfdiffusivities of CH4 in an LTA-type zeolite at several temperatures, as shown in Figure 3. For CH4, the self-diffusivity can be calculated from conventional MD simulations at all the investigated temperatures, although a decrease in temperature leads to the difficulty of sampling sufficiently to compute diffusivities.
Figure 4. Calculated thermodynamic quantities along the reaction coordinate: (a) potential energy profile, (b) internal energy profile, (c) Helmholtz free energy profile, and (d) entropy profile (multiplied by -T), for CH4 for two molecules/cage in an LTA-type zeolite at 800 to 1800 K. For each figure, the middle point of the horizontal axis is located at the cage center and both ends are at the center of the window .
In this figure, the diffusivity at 800 K is ca. 3.3 × 10-9 m2/s. On the other hand, the self-diffusivity calculated by MD at 750 K in a 5A zeolite is ca. 6 × 10-9 m2/s.36 Furthermore, the experimental value obtained by a neutron spin echo technique at 475 K is also on the order of 10-9 m2/s.37 It should be noted that these previous reports study the diffusion in a cationexchange zeolite at the other temperature than 800 K, but the fact that all these data are in the order of 10-9 m2/s supports the validity of our simulated values. In Figure 3, when we assume the Arrhenius temperature dependence, the activation energy was estimated at 14.3 kJ/ mol, although these points show a somewhat downward convex curve. This slope certainly shows activated diffusion, corresponding to an apparent activation energy. In Figures 3 and 16 of ref 28, the self-diffusivities at 300 and 600 K are shown for two molecules/cage. In these figures, the diffusivity at 300 K is 0.11-0.13 × 10-9 m2/s, and that at 600 K is 1.1-1.3 × 10-9 m2/s, respectively. Thus the apparent activation energy is ca. 9-12 kJ/mol, although it was estimated from only two points. It is reasonable that the value is smaller than that estimated from 800 to 1800 K, because the curve obtained from the interpolation in Figure 3 is somewhat downward convex. The reason for this nonlinearity is discussed in the next section. 4.2. Profiles of Thermodynamic Quantities. The thermodynamic quantities, potential energy V(T, x), internal energy U(T, x), free energy A(T, x), and entropy multiplied by the negative of the temperature, -TS(T, x) were calculated at each temperature using eqs 15-17, as shown in Figure 4a-d. The horizontal axis represents the reaction coordinate from one center of the window to a neighboring center. Thus, the middle point of the axis is located at the center of the R-cage. As for A(T, x) and -TS(T, x), to compare the shapes of these profiles at different temperatures, the values were set at 0 kJ/mol at the original point of the reaction coordinate. These quantities show a maximum at the position of the center of the window.
Slow Molecular Diffusion Rates in Zeolite Nanopores
J. Phys. Chem. C, Vol. 112, No. 8, 2008 2809
TABLE 1: Temperature Dependence of the Differencea between the Maximum and the Minimum Helmholtz Free Energy, ∆A, Internal Energy, ∆U, Entropy Multiplied by -T, -T∆S, and Potential Energy, ∆V, for CH4 for Two Molecules/Cage in an LTA-Type Zeolite T (K)
∆A
∆U
-T∆S
∆V
1800 1600 1400 1200 1000 800
68.6 62.7 55.5 49.3 42.8 35.6
11.7 9.70 9.21 8.57 7.76 7.38
58.2 53.5 47.0 41.3 35.6 29.3
12.6 11.2 10.6 9.70 8.85 8.30
a Values are in kJ/mol. b All these data were calculated from MD simulations.
TABLE 2: Transmission Coefficient K for CH4 for Two Molecules/Cage in an LTA-Type Zeolite, with Statistical Error σ, Computed from n MD Runs
a
tempa
k
(σ
n
1800 1600 1400 1200 1000 800
0.8850 0.8861 0.8720 0.8885 0.8818 0.8635
(0.0253 (0.0319 (0.0318 (0.0291 (0.0264 (0.0281
1565 983 992 1184 1438 1267
Values are in kelvin.
To gain more insight into these quantities, the differences between the maximum and the minimum of these quantities along the reaction coordinate, ∆A, ∆U, -T∆S, and ∆V, are listed in Table 1. Generally, these differences become larger with an increase in temperature. ∆U is a little smaller than ∆V, resulting from the fact that K is almost constant along the reaction coordinate but actually a little smaller at the center of the window. Considering that both of V and U show a maximum at the center of the window and a minimum in the cage regions, it is concluded that K is almost constant along the reaction coordinate at a given temperature. Moreover, this invariance of K shows that the guest molecules diffuse at almost constant speed along the reaction coordinate. On the other hand, it is probably natural to consider that the increase in ∆V with temperature results from the increase of the thermal vibration of guest molecules. As a result, the collision frequency between molecules and nanopores is raised and the potential energy increases drastically, particularly around the center of the window, where the free volume is remarkably small. Therefore, ∆V also increases with an increase of temperature. In particular, in Table 1, ∆V at 1800 K is about 1.5 times higher than at 800 K. The temperature dependence of ∆V probably correlates with the downward convex curve for the temperature dependence of the self-diffusivities shown in Figure 3. This suggests that increasing the temperature of the system has a greater effect on the thermal motion of each guest molecule than expected from the Arrhenius behavior. This certainly leads to higher self-diffusivities than estimated from the extrapolation over a range of lower temperatures. 4.3. Transmission Coefficient. The transmission coefficient, κ, was computed at each temperature using eq 4. The results are summarized in Table 2. The number of MD runs, n, and the error, σ, obtained from eq 5 are also shown. Table 2 suggests that κ is independent of the temperature between 800 and 1800 K for CH4 in an LTA-type zeolite. We consider that all the errors listed in Table 2 are small enough for substituting into eq 2, because all the ratios of σ to κ are within (0.04. It is noted that the temperature dependence of κ can be changed
Figure 5. Comparison of the Helmholtz free energy profile calculated by MD at temperature T and the profiles predicted from the different temperatures, for CH4 for two molecules/cage in an LTA-type zeolite from 800 to 1800 K. (a) T ) 1800 K, (b) T ) 1600 K, (c) T ) 1400 K, (d) T ) 1200 K, (e) T ) 1000 K, and (f) T ) 800 K.
drastically according to the concentration conditions and the types of guest species or host zeolite. 4.4. Prediction of the Free Energy Profile. Performing MD simulations at a given temperature, we predicted Helmholtz free energy profiles at different temperatures, using eq 14. The predicted profiles are shown in Figures 5a-f, where we added the profiles obtained from MD simulations at the corresponding temperatures, which are shown in Figure 4c. In these figures, the horizontal axis represents the reaction coordinate from one center of the window to the adjacent center, and the values of free energy were set at 0 kJ/mol at the original point of the reaction coordinate. In Figure 5, all the predicted profiles show excellent agreement with those obtained from MD at the corresponding temperatures. In particular, it should be noted that we succeeded in predicting the free energy profiles at higher temperatures from MD calculations performed at lower temperatures, as well as the profiles at lower temperatures from MD simulations carried out at higher temperatures. This reversibility on the free energy profile at different temperatures confirms that the temperature dependence of Helmholtz free energy profiles can be predicted over a wide range of temperature by performing MD at only one temperature. In the future, further investigation into the applicable temperature range of our methodology will be required. In eq 14, it is important how accurately the ensemble average can be calculated from MD simulations. Sufficient data sampling of the potential energy VN is needed to achieve the accurate
2810 J. Phys. Chem. C, Vol. 112, No. 8, 2008
Figure 6. Comparison of the self-diffusivities calculated by MD at each temperature and the predicted values using the free energy profiles estimated from the different temperatures, for CH4 in an LTAtype zeolite. The concentration of guest molecules is two molecules/ cage.
estimation of the average. However, in Figure 5, we confirmed that the predicted free energy profiles show excellent agreement with those obtained from MD at the corresponding temperature. This agreement shows that the ensemble average in eq 14 can be estimated accurately. 4.5. Self-Diffusivity. Using the results shown in Figure 5, we can calculate TST rate constants from eq 6. From these TST rates and the transmission coefficients listed in Table 2, hopping rates k were calculated from eq 2. Consequently, the selfdiffusivities were estimated from eq 1, as shown in Figure 6. In this figure, the calculated diffusivities in Figure 3 are also shown for comparison with the diffusivities predicted from our methodology. The results from our method show excellent agreement with the calculated diffusivities over a wide range of temperatures. For example, the self-diffusivity of CH4 at 800 K predicted from MD at 1800 K is 3.1 × 10-9 m2/s, whereas that simulated from MD at 800 K is 3.3 × 10-9 m2/s. This excellent agreement certainly confirms our methodology. In Figure 3, it is shown that the Arrhenius dependence cannot be always assumed over a wide range of temperatures. Therefore, our methodology will be applied to various diffusion phenomena, regardless of the Arrhenius dependence. Moreover, it strongly suggests that by performing the free energy calculation at a given temperature and computing transmission coefficients, the temperature dependence of molecular diffusivity is comprehensively predicted, including the case of the completely ‘slow’ diffusion regime. 5. Conclusion In this study, we estimated molecular diffusivities for CH4 through an LTA-type zeolite at a concentration of two molecules per cage, combining statistical-mechanical theory with conventional MD calculations. In our method, self-diffusivity is calculated by TST, where it is essential to estimate the free energy profile along the reaction coordinate. However, it is quite difficult to take sufficient statistical samplings for the free energy calculation, especially in the case of the completely ‘slow’ molecular diffusion regime. Therefore, in our methodology, by heating the system the molecular diffusion rate is accelerated in MD simulations. By carrying out MD at higher temperatures, free energy profiles at lower temperatures can be predicted. Using these profiles, self-diffusivities can be calculated by TST over a wide range of temperatures. Here, we have given details of our methodology and the data necessary to estimate self-diffusivities, i.e., transmission coef-
Nagumo et al. ficients and free energy profiles, as well as thermodynamic quantities, such as potential energy, internal energy, and entropy. These data enable us to investigate the molecular diffusion mechanism from a thermodynamic viewpoint. The reason that self-diffusivities at higher temperatures become larger than expected from a linear extrapolation at lower temperatures is certainly due to the fact that the difference between the maximum and minimum of the potential energy profiles becomes larger with an increase in temperature. Acturally, cation-exchange zeolites are of more practical importance. A set of potential functions and interaction parameters has been developed for cation-exchange zeolite in a previous study.38 Using such a set of parameters, the efficient prediction of molecular diffusivities inside various types of zeolite will be achieved. Our approach will open the way to achieving the systematic prediction of ‘slow’ molecular diffusion rates within various types of nanoporous material, performing conventional MD simulations at much lower computational cost. In future work, we will apply our methodology to various diffusion phenomena of other guest-host systems. References and Notes (1) Ka¨rger, J.; Ruthven, D. M. Diffusion in Zeolite and Microporous Solids; Wiley-Interscience: New York, 1992. (2) Davis, M. E. Nature 2002, 417, 813-821. (3) Takaba, H.; Koyama, A.; Nakao, S. J. Phys. Chem. B 2000, 104, 6353-6359. (4) Takaba, H.; Nakao, S. J. Chem. Eng. Jpn. 2003, 36, 1364-1369. (5) Takaba, H.; Mizukami, K.; Oumi, Y.; Kubo, M.; Chatterjee, A.; Fahmi, A.; Miyamoto, A. Catal. Today 1999, 50, 651-660. (6) Takaba, H.; Koshita, R.; Mizukami, K.; Oumi, Y.; Ito, N.; Kubo, M.; Fahmi, A.; Miyamoto, A. J. Membr. Sci. 1997, 134, 127-139. (7) Takaba, H.; Mizukami, K.; Kubo, M.; Stirling, A.; Miyamoto, A. J. Membr. Sci. 1996, 121, 251-259. (8) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. A 2003, 107, 1013210141. (9) Skoulidas, A. I.; Sholl, D. S. J. Phys. Chem. B 2002, 106, 50585067. (10) June, R. L.; Bell, A. T.; Theodorou, D. N. J. Phys. Chem. 1992, 96, 1051-1060. (11) Nagumo, R.; Takaba, H.; Nakao, S. Microporous Mesoporous Mater. 2001, 48, 247-254. (12) Suzuki, S.; Takaba, H.; Yamaguchi, T.; Nakao, S. J. Phys. Chem. B 2000, 104, 1971-1976. (13) Nagumo, R.; Takaba, H.; Nakao, S. J. Phys. Chem. B 2003, 107, 14422-14428. (14) Voter, A. F.; Montalenti, F.; Germann, T. C. Annu. ReV. Mater. Res. 2002, 32, 321-346. (15) Sorensen, M. R.; Voter, A. F. J. Chem. Phys. 2000, 112, 95999606. (16) Montalenti, F.; Voter, A. F. J. Chem. Phys. 2002, 116, 48194828. (17) Montalenti, F.; Sorensen, M. R.; Voter, A. F. Phys. ReV. Lett. 2001, 87, 126101. (18) Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. J. J. Phys. Chem. 1996, 100, 12771-12800. (19) Truhlar, D. G.; Hase, W. L.; Hynes, J. T. J. Phys. Chem. 1983, 87, 2664-2682. (20) Garcia-Viloca, M.; Gao, J.; Karplus, M.; Truhlar, D. G. Science 2004, 303, 186-195. (21) Mosell, T.; Schrimpf, G.; Brickmann, J. J. Phys. Chem. 1996, 100, 4582-4590. (22) Mosell, T.; Schrimpf, G.; Hahn, C.; Brickmann, J. J. Phys. Chem. 1996, 100, 4571-4581. (23) Mosell, T.; Schrimpf, G.; Brickmann, J. J. Phys. Chem. B 1997, 101, 9476-9484. (24) Mosell, T.; Schrimpf, G.; Brickmann, J. J. Phys. Chem. B 1997, 101, 9485-9494. (25) Smit, B.; Loyens, L. D. J. C.; Verbist, G. L. M. M. Faraday Discuss. 1997, 106, 93-104.
Slow Molecular Diffusion Rates in Zeolite Nanopores (26) Beerdsen, E.; Smit, B.; Dubbeldam, D. Phys. ReV. Lett. 2004, 93, 248301. (27) Dubbeldam, D.; Calero, S.; Maesen, T. L. M.; Smit, B. Phys. ReV. Lett. 2003, 90, 245901. (28) Dubbeldam, D.; Beerdsen, E.; Vlugt, T. J. H.; Smit, B. J. Chem. Phys. 2005, 122, 224712. (29) White, G. W. N.; Goldman, S.; Gray, C. G. Mol. Phys. 2000, 98, 1871-1885. (30) Anderson, J. B. AdV. Chem. Phys. 1995, 91, 381-431. (31) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic Press: London, 2002. (32) Ruiz-Montero, M. J.; Frenkel, D.; Brey, J. J. Mol. Phys. 1997, 90, 925-941.
J. Phys. Chem. C, Vol. 112, No. 8, 2008 2811 (33) Schu¨ring, A.; Auerbach, S. M.; Fritzsche, S.; Haberlandt, R. J. Chem. Phys. 2002, 116, 10890-10894. (34) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684-3690. (35) Hu, J.; Goldman, S.; Gray, C. G.; Guy, H. R. Mol. Phys. 2000, 98, 535-547. (36) Krishna, R.; van Baten, J. M. J. Phys. Chem. B 2005, 109, 63866396. (37) Jobic, H; Me´thivier, A.; Ehlers, G.; Farago, B.; Haeussler, W. Angew. Chem., Int. Ed. 2004, 43, 364-366. (38) Garcı´a-Pe´rez, E.; Dubbeldam, D.; Maesen, T. L. M.; Calero, S. J. Phys. Chem. B 2006, 110, 23968-23976.
