Modeling Diffusion of Linear Hydrocarbons in Silica Zeolite LTA Using

Article pubs.acs.org/JPCC

Modeling Diffusion of Linear Hydrocarbons in Silica Zeolite LTA Using Transition Path Sampling Salah Eddine Boulfelfel,† Peter I. Ravikovitch,*,‡ and David S. Sholl*,† †

School of Chemical and Biomolecular Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332-0100, United States Corporate Strategic Research, ExxonMobil Research and Engineering, 1545 Route 22 East, Annandale, New Jersey 08801, United States

‡

S Supporting Information *

ABSTRACT: The diffusivities of linear hydrocarbons (CH4, C2H6, C2H4, C3H8, C3H6, and C4H10) in pure silica zeolite LTA (ITQ-29) are computed at 300 K and infinite dilution. To overcome the time scale problem arising from the slow diffusion process at room temperature, we used transition path sampling (TPS). The influence of framework flexibility on diffusion is investigated by combining TPS simulations with fully flexible molecular dynamics performed in the NpT ensemble. The ensemble of the collected reactive trajectories was used to characterize sets of transition states, and the corresponding configurations were analyzed to construct window size distributions during the molecular hopping events. The diffusion process is affected by framework flexibility, and the influence of framework flexibility on diffusion of propane and butane is much larger than for methane and ethane.

1. INTRODUCTION Microporous crystalline materials such as zeolites are widely used in a variety of petrochemical applications.1 They are employed as adsorbents, catalysts, molecular sieves, or ion exchangers in many chemical processes.2−4 In many of these processes the pore size of the material is comparable to the size of the adsorbed molecules. This implies very different transport and thermodynamic properties of the molecules inside the porous material as compared to the bulk.5 Numerous efforts have been made to synthesize zeolites with different pore dimensions6 to match the needs of various applications. Structures with cages and channels connected through smallpore eight-membered rings (8MRs) are potential candidates for processes such as hydrocarbon separations.7−12 The performance of such processes is influenced by the diffusion rate of guest molecules through the zeolite pores.8,9,13 Therefore, an accurate description of diffusion at the molecular level is imperative for a fundamental understanding of transport processes inside zeolites and for making quantitative predictions for materials screening. Computationally, obtaining accurate diffusion coefficients can be a challenging task.14−16 When diffusion is sufficiently fast, diffusivities can be determined from straightforward molecular dynamics (MD) simulations. However, in many situations of interest, particularly for guest molecules in porous materials with small pores, MD simulations cannot achieve sufficiently long time scales to make molecular diffusion observable. To speed up MD simulations, and sometimes due to the lack of a reliable force field, the vibrations of the material framework are neglected, and molecules in MD simulations propagate inside a rigid structure.17,18 However, when the pore © 2015 American Chemical Society

dimension is comparable to or smaller than the size of the diffusing molecule, framework flexibility can play an important role in diffusion and cannot be disregarded. 19,20 To approximate the effect of framework flexibility, the use of a time-averaged structure as a rigid configuration for MD simulation has been proposed.21,22 Recently, Awati et al.23 introduced a more accurate method to include framework flexibility by performing MD simulations in a set of precalculated rigid framework configurations, replacing configurations after a predetermined number of MD steps. A class of methods suitable for predicting slow diffusion are based on transition-state theory.24 They can be dynamically corrected using the Bennett−Chandler approach.25 Typically, transition-state theory (TST)-based methods when applied to nanoporous crystals are used with a rigid framework, and the calculations rely on the definition of a dividing surface.26−28 Most applications of this class of methods for zeolites have been focused on molecules with a shape that can be approximated as spherical, e.g., CH4 and H2. New modifications of TST-based methods have been developed by Haldoupis et al.29 and Awati et al.23 to assess the role of framework flexibility in molecular diffusion of CH4 and CO2 in ZIF-8 and CH4 in 8MR zeolites, respectively. Despite recent advances in the simulation of molecular diffusion, and the development of new computational tools capable of assessing framework flexibility, studies that investigate the diffusion of larger molecules are still limited.26,30 Received: February 17, 2015 Revised: May 28, 2015 Published: May 28, 2015 15643

DOI: 10.1021/acs.jpcc.5b01633 J. Phys. Chem. C 2015, 119, 15643−15653

Article

The Journal of Physical Chemistry C

trajectory length), and the other is time dependent. The time derivative of C(t) gives the reaction rate:

To tackle molecular diffusion of large adsorbates in a fully flexible framework, it is possible to use path-based approaches such as transition path sampling31−33 (TPS) and transition interface sampling.34,35 In these methods, hopping events are generated using a Monte Carlo algorithm in a dynamical path space. MD simulations are used to produce short dynamical trajectories, and the net hopping rate constant is computed from the correlation function of the sampled path ensemble.31−35 Vlugt et al.36 demonstrated how to use TPS to compute the diffusivity of isobutane in silicalite. Although in their study the porous structure was considered rigid, it is possible to perform the sampling with a fully flexible structure. In this work, we investigate diffusion of methane, ethane, ethylene, propane, propylene, and butane in pure silica LTA, a prototypical 8MR zeolite, using atomistic modeling. Our main objectives are to compute diffusion coefficients for each species and to understand the relationship between framework structure flexibility and molecular diffusion. To this end, we use molecular-dynamics-enabled TPS simulations.37 Zeolite LTA (also known as ITQ-29 in its pure silica form38) has been extensively studied computationally.19,21−23,39−44 However, only limited experimental data on diffusion of different molecules are available.11−13 LTA has a symmetrical cubic structure featuring nearly spherical cages connected via 8MR windows.38 The choice of LTA is also motivated by its structural anharmonicity; i.e., the time-averaged and energyminimized minimum window dimensions are unequal,23 which limits the use of a rigid framework approximation.23,29

k=

⟨hA (x0) hB(xt )⟩ ⟨hA (x0)⟩

(4)

⟨hB(t)⟩FAB(x0,T) and ⟨hB(t′)⟩FAB(x0,T) are calculated from separate simulations and FAB(x0 , T ) ≡ ρ(x0) hA (x0) HB(x0 , T )

(5)

Here HB = 1 if the system starts at t = 0 in region A and visits region B at least once prior to t = T. Since transitions from region A to region B are rare, C(t′) is expected to be small in region B, and can be expressed as the probability of finding the system in region B after time t = t′: P(λ , t ′) =

∫ dx0 f (x0 , t ′) δ[λ − λ(xt ′(x0))] ∫ dx0 f (x0 , t ′)

(6)

where λ is a reaction coordinate describing the progress of the transition. Sampling techniques such as umbrella sampling46−48 can be used to compute this probability. 2.2. Transition Path Sampling. To perform transition path sampling simulations, we use two main trial moves to sample the trajectory phase space. 2.2.1. Shooting Moves. TPS is started from an initial trajectory connecting stable states A and B (see Figure 1). A time slice is randomly chosen, and a perturbation δp drawn from a Gaussian distribution is added to its atomic momenta. This move is accepted or rejected according to the probability

2. COMPUTATIONAL METHODS In the following section, the theoretical background of the transition path sampling method31−33,37 is briefly reviewed with the focus on sampling trajectories with deterministic dynamics.32 The basic concepts for the calculation of hopping rates,31,33,45 necessary for diffusion constant computation, and the identification of transition states36,37 are also presented. 2.1. Reaction Rate. For a dynamical system at equilibrium with two stable states, A and B, kinetic information about the transition between these states can be derived from the correlation function of the state population fluctuations: C(t ) ≡

d[⟨hB(t )⟩FAB(x0 , T )] dC(t ) C(t ′) = dt ⟨hB(t ′)⟩FAB(x0 , T ) dt

(1)

where hA and hB are state functions with hA,B(x) = 1 if x ∈ A, B and hA,B(x) = 0 otherwise, ⟨...⟩ denotes an average over the equilibrium ensemble, and xt = qt, pt is a set of atomic coordinates and atomic momenta specifying the system at time t. According to the fluctuation−dissipation theorem37,46 C(t ) ≈ ⟨hB⟩(1 − e−t / τrxn) ≈ kA → Bt

(2)

where τrxn is the characteristic reaction time for the system A → B and kA→B is the corresponding rate constant. Because dynamical trajectories are fully determined by the initial conditions, x0, the ensemble average in eq 1 can be rewritten in terms of the equilibrium phase space distribution ρ(x0): C(t ) =

∫ dx0 ρ(x0) hA (x0) hB(xt ) ∫ dx0 ρ(x0) hA (x0)

Figure 1. Schematic representation of the simulation system. The red/ yellow tubes represent the pure silica LTA zeolite framework. The large blue ball represents a CH4 molecule as a united atom. The trace of the molecule along the dynamical hopping trajectory between cage A (solid circle) and cage B (dashed circle) is represented as a chain of small white spheres.

(3)

The expression of the correlation function C(t) can be factorized into two terms: one is time independent and is computed for a chosen time t′ (t′ < T; T is the maximal 15644


Article

The Journal of Physical Chemistry C n→o ⎞ ⎛ F(x0n , T )Pgen o→n ⎟ Pacc = min⎜⎜1, o o→n ⎟ ⎝ F(x0 , T )Pgen ⎠

pB(N ) (r ) ≈ (7)

⎛ F(x0n , T ) ⎞ o→n Pacc = min⎜1, ⎟ o ⎝ F (x 0 , T ) ⎠

(8)

(11)

(12)

3. SIMULATION DETAILS In this work, the TPS method was implemented as a standalone program that can in principle be combined with any molecular dynamics code. In our case, all MD calculations were performed with LAMMPS. 3.1. Defining the Stable States. TPS relies on an order parameter to distinguish stable states A and B and, therefore, state functions hA and hB. In this work, we use λ, the ratio of the difference to the sum of the distances (dA,dB) between the center of mass of the molecule and the geometrical centers of the cages (considered as spheres) associated with states A and B, as illustrated in Figure 1:

(9)

2.3. Transition State Ensemble. To shed light on the effect of framework flexibility on molecular diffusion in LTA, we need to examine configurations from dynamical trajectories at times when molecules cross the 8MR. This implies identification of transition states and requires a special procedure to identify them from a large number of trajectories harvested during TPS simulations. In a simple system, the transition state is defined as a saddle point in the potential energy surface (PES). However, the saddle points of a more complex system cannot be easily located and enumerated. Transition-state identification in systems of this type is affected by the choice of a reaction coordinate describing the progress of the transition A → B. Unlike in TST-based methods, TPS does not require a priori knowledge of the hopping mechanism, and the dividing surface and transition states are the output of the simulation rather than inputs. To determine transition states, configurations from which the system relaxes into one or the other stable state with equal probability are identified. Since we do not rely on a predefined reaction coordinate, this probability is quantified from an ensemble of dynamical trajectories using the committor probability49

∫ dx(ts) P[x(ts)] δ(r0 − r ) hB(xts) ∫ dx(ts) P[x(ts)] δ(r0 − r )

i=1

NB N

Here, σA ≈ (pA(1 − pA)/N)1/2 and σB ≈ (pB(1 − pB)/N)1/2 are errors on the estimated probabilities pA and pB, respectively. For an accurate identification of a transition state, σA and σB should be chosen as small as possible. Therefore, the only way to satisfy eq 12 is to have pA = pB = 0.5, which leads to one state (or configuration) along a single trajectory.

This implies that if a trajectory is initiated in region A at time t = 0, and visits region B at least once during time T, the move is accepted. After application of a shooting move to a time slice, the initial conditions x0 are set by integrating the equations of motion back to time t = 0. 2.2.2. Shifting Moves. These moves consist of translating the initial conditions by an amount Δt that can be drawn from a Gaussian distribution. The corresponding acceptance probability can be written as

pB (r , ts) =

s

|pA − pB | < σA + σB

⎛ ρ(x0n) ⎞ n n = min⎜1, o hA (x 0 ) HB(x 0 , T )⎟ ⎝ ρ (x 0 ) ⎠

⎛ ρ(x0n) ⎞ o→n n n = min⎜1, Pacc o hA (x 0 ) HB(x 0 , T )⎟ ⎝ ρ (x 0 ) ⎠

N

∑ hB(xt(i)) ≈

where N represents the number of trajectories initiated from a configuration r selected from a reactive path collected during TPS simulations. NB is the number of trajectories that propagate into stable state B. A point is considered as a transition state if it satisfies

n→o where Po→n gen and Pgen are the probabilities of generation of a new (n) configuration from the old (o) one and vice versa, respectively. For deterministic dynamics, the acceptance probability can be written as

= hA (x0n) HB(x0n , T )

1 N

λ=

dA − dB dA + dB

(13)

By definition, the order parameter λ varies between −1 and +1. Negative values of this parameter correspond to the configuration belonging to stable state A, while positive values are associated with region B. We performed a set of long MD runs (>1000 ps) with different initial conditions starting in both regions to determine the range of λ corresponding to each region. As a result, the stable states A and B are defined by intervals −1.0 ≤ λA ≤ −0.1 and 0.1 ≤ λB ≤ 1.0, respectively (Figure 2). 3.2. Force Fields. The united atom (UA) TraPPE (transferable potential for phase equilibria) model50,51 is used to describe internal degrees of freedom in the hydrocarbon molecules studied (CH4,50 C2H6,50 C2H4,51 C3H8,50 C3H6,51 and C4H1050). The UAs are charge neutral, and intramolecular interactions include bond stretching, angle bending, and torsion. Nonbonded interactions are modeled using a 12−6 Lennard-Jones pair potential with a cutoff of 12 Å. This set of force field parameters was demonstrated to accurately reproduce the bulk equilibrium properties of a series of linear and branched hydrocarbons.50,51 The flexible LTA zeolite was described by the Hill−Sauer forcefield,52,53 which includes bond stretching (Si−O), angle bending (Si−O−Si and O−Si−O), and torsion (O−Si−O−Si) and coupling terms between these degrees of freedom. The charges assigned to Si and O ions are q (e) = 0.5236 and −0.2618, respectively, and a 9−6 Lennard-Jones potential with a cutoff of 12 Å is used for short-range nonbonded interactions. Ewald summation was used for computing long-range Coulomb interactions.47

(10)

where r is the configuration space vector and ts is the time at the transition state. The committor probability is especially useful for the determination of transition states in flexible framework simulations and also for molecules with irregular or nonspherical shape. For TPS combined with deterministic dynamics, we can define 15645


Article


temperature can be used to increase the frequency of the transition event and obtain a first pathway.49,70 However, lowering the temperature to the desired value during the sampling can lead to frequent unsuccessful transition paths, reducing the overall efficiency of the calculations. The symmetry of the stable states A and B can be used to generate geometric intermediate configurations to derive a first trajectory.71−74 In this work, the molecule was placed at the center of one of the cubic cell faces at crystallographic special position 3c (0, 1/2, 1/2), which corresponds to the center of the 8MR. Then the system was propagated forward (+t) and backward (−t) in time. If the corresponding MD trajectories do not each end in one of the stable states A and B, the previous procedure is repeated with a different velocity distribution until an initial successful trajectory is obtained. 3.5. Collecting the Trajectory Ensemble. In this work, each TPS run involved at least 104 cycles. For each cycle, a shooting move (90% probability) or a shifting move (10% probability) was applied. The acceptance of the combined moves was at least 30% for each TPS run. No other type of move (e.g., displacement of the molecule or trajectory reversal33,37) was applied during the sampling. The average length of trajectories was 5 ps for methane and ethane and 10 ps for propane and butane. The length of the trajectories was chosen to ensure that both ends of each path, from forward and backward integration in time, are well localized in states A and B. The computational effort required to generate these ensembles of trajectories should be compared to the effort required to simulate molecular diffusion using standard MD. Although the sampling of the trajectory ensemble of methane was more time-consuming than using straightforward MD, the efficiency of TPS exceeds that of MD as we move to larger molecules such as propane and butane. It was demonstrated that the typical cage-to-cage hopping time for ethane in LTA at 300 K is around 500−1000 ps.30 To accurately determine this rate using standard MD, at least five hopping events per molecule are desirable. Although this is reasonable for methane, the number of cage-to-cage hopping events recorded for propane and butane are almost equal to zero using standard MD for the same computing time as TPS at 300 K. 3.6. Hopping Rate and Diffusion Constant. The rate constants are obtained from the time derivative of the correlation function of the trajectory ensemble as explained in section 2.2. Diffusion constants at infinite dilution are computed using Dx = (1/2)kxa2, Dy = (1/2)kyb2, and Dz = (1/ 2)kzc2. Here, a, b, and c are the unit vectors in the diffusion lattice of pure silica zeolite LTA (a = b = c = 12.04 Å). The net diffusivity D is given by D = (Dx + Dy + Dz)/3.36 Diffusion is investigated under the assumption that two successive hopping events are uncorrelated. Molecular dynamics simulations performed on the diffusion of small hydrocarbons showed that the time needed to observe two successive hopping events is much larger than the total length of each of the trajectories collected during the sampling.30 The waiting times between these events increases dramatically as the size of the molecule gets larger. On the other hand, the analysis of the path ensemble of methane showed than only a few trajectories (less than 10 out of 10 000) displayed 2 successive hopping events that belonged to the same trajectory. For ethane, propane, and butane, no such behavior was recorded during the sampling. This strongly supports the assumption of uncorrelated hopping.

Figure 2. Variation of the order parameter λ defined as λ = (dA − dB)/ (dA + dB). The data are from a trajectory of one CH4 molecule hopping from cage A to cage B over a period of time equal to 7 ps. Regions A and B are associated with negative and positive values of the order parameter, respectively. For λ = 0.0, the molecule is at the center of the 8MR separating the two cages.

For adsorbate−adsorbent interactions, many sets of interactions were reported in the literature.54−62 We used LennardJones parameters from refs 58 and 62: CH4−O [ε/kB = 115 K, σ = 3.47 Å], CH3−O [ε/kB = 93 K, σ = 3.48 Å], CH2−O [ε/kB = 60.5 K, σ = 3.58 Å]. These parameters were optimized to reproduce inflection points on the adsorption isotherms of linear and branched alkanes.62 A more detailed description of the force field models and the interaction parameters used in this work is presented in the Supporting Information. 3.3. Molecular Dynamics Simulations. Classical molecular dynamics simulations of flexible LTA zeolite were carried out using LAMMPS.63 The equations of motion used are those of Shinoda et al.,64 which combine the hydrostatic equations of Martyna, Tobias, and Klein65 with the strain energy proposed by Parrinello and Rahman.66 The velocity-Verlet algorithm47 is used to integrate the equations of motion. To ensure good time reversibility, an integration time step of 0.2 fs is used. All simulations are performed at temperature T = 300 K and pressure P = 1 bar. The thermostat and barostat are applied to the system using Nose−Hoover algorithms67,68 implemented as described by Melchiona et al.,69 allowing anisotropic shape changes of the simulation box. The simulation box contained 1944 atoms (648 Si and 1296 O atoms) in addition to one linear hydrocarbon. This corresponds to a 3 × 3 × 3 supercell of the LTA zeolite. For the calculation of diffusion constants from standard MD simulations, we performed a set of 30 MD runs for methane, ethane, and ethylene. Each run was propagated for at least 200 ns with an integration time step equal to 0.5 fs using the velocity-Verlet algorithm in the NpT ensemble. For methane, an LTA 2 × 2 × 2 supercell was used with 16 noninteracting molecules. For ethane and ethylene, the same supercell contained 8 molecules. The mean-squared displacement plots are reported in Figures S1 and S2 in the Supporting Information. 3.4. Generating a First Trajectory. The transition path sampling procedure develops from an initial trajectory connecting stable regions A and B. The slow diffusion rate at the target temperature (T = 300 K) can make the task of obtaining such a trajectory difficult. Molecular dynamics at high 15646


Article


results of all molecules are summarized in Table 1. The first step in the calculation of a rate constant from TPS consists of performing a Monte Carlo sampling of the distribution FAB(x0,T) from eq 5. A subset of the collected trajectories is illustrated in Figure 4. The function ⟨hB(t)⟩FAB(x0,T) associated

3.7. Identification of Transition States. Every 100 Monte Carlo moves during the sampling of the distribution FAB(x0,T) from eq 5 in section 2.1, a successful trajectory is analyzed according to the procedure described in section 2.3. For each configuration along the chosen trajectory, 100 MD runs (N = 100) are initialized by randomly assigning velocities using a Maxwell−Boltzmann distribution. Each run is propagated for at least 20 000 MD steps with a time step equal to 0.2 fs. In Figure 3, we have plotted the variation of the committor probability used for the determination of transition states for

Figure 4. Variation of the order parameter λ for trajectories collected during the sampling of CH4 hopping in LTA. For clarity, only 200 trajectories are shown.

Figure 3. Variation of the committor probability pB as a function of the order parameter λ along a transition path for CH4 in silica LTA zeolite. The committor probability is computed by the determination of the fraction of fleeting trial trajectories initiated with randomized momenta that reach region B in a time ts = 4 ps. More than 100 fleeting trajectories were generated for each time slice to obtain pB with sufficient accuracy.

methane in LTA. Along the transition path, any configuration that satisfies the condition |pA − pB| < σA + σB is considered as a transition state. For each identified transition state, the window size of the starting configuration is recorded for later analysis. The window size is identified in the eight-membered ring crossed by the molecule. Only two O−O distances are considered, vertical and horizontal (diagonals are not included). The ionic radius of one oxygen ion (2.7 Å) is subtracted from the distance between the centers of the oxygen ions.

Figure 5. Function ⟨hB(t)⟩FAB(x0,T) from TPS simulations of CH4 hopping in LTA.

4. RESULTS The different steps of the procedure of computing the rate constant and diffusivity are explained for methane, and then the

with this ensemble is plotted in Figure 5 for methane with T = 5 ps. Next, the time derivative of the previous function is

Table 1. Diffusion Constants for Linear Hydrocarbons in LTA at Infinite Dilution at 300 K from Transition Path Sampling and Standard MD Simulationsa molecule CH4 C2H6 C2H4 C3H8 C3H6 C4H10

diffusion constant (TPS) (m2·s−1) (1.95 (1.15 (1.04 (2.30 (9.30 (7.1

± ± ± ± ± ±

0.45) 0.54) 0.95) 0.70) 0.98) 1.10)

× × × × × ×

acceptance rate (%)

diffusion constant (standard MD) (m2·s−1)

65 59 55 50 54 61

(2.10 ± 0.86) × 10−11 (1.76 ± 0.73) × 10−12 (3.45 ± 0.79) × 10−12 b b b

10−11 10−12 10−12 10−13 10−13 10−16

a

TPS results are based on the analysis of path ensembles of at least 10 000 trajectories for each molecule. bMD results are not reported due to the very slow diffusion at 300 K. 15647


Article


obtain for methane k = 5.52 × 107 s−1. Using eqs 4 and 5, we obtain D = 1.96 × 10−11 m2·s−1. For validating the results obtained for CH4, we performed standard MD simulations at 300 K for the infinite dilution case (details are given in section 3.3), and the obtained diffusion constant was D = 2.10 × 10−11 ± 0.86 × 10−12 m2·s−1. This shows an adequate agreement between TPS and standard MD results. The TPS results for ethane, propane, and butane are summarized in Table 1 with a comparison to diffusivities obtained from standard MD simulations. TPS simulations are able to access diffusivities for adsorbate−adsorbent systems far beyond the time scales accessible to standard MD. This is best illustrated by our results for n-butane, which is found to diffuse more than 4 orders of magnitude slower than methane at 300 K. It is worth noting that significant error bars are introduced by the calculation of P(λ,t′) in the computation of rate constants. This is mainly due to the sensitivity of umbrella sampling calculations to the choice of the reaction coordinate and the definition of the interval delimiting the stable region B. In this study, the parameter λ defined in section 3.1 better captures the progress of the hopping of methane than that of propane or butane. The diffusivities of propane and butane are on the order of 10−13−10−14 and 10−15−10−16 m2·s−1, respectively.

computed as shown in Figure 6. Then the probability P(λ,t′) is computed for t′ = 3 ps using umbrella sampling and integrated

Figure 6. Time derivative of the function ⟨hB(t)⟩FAB(x0,T) from transition sampling simulations of CH4 hopping in LTA.

over region B (t′ must be chosen shorter than the total trajectory length but long enough to observe a hopping event). By combining the results from the second and third steps, we

Figure 7. Variation of the two minimal O−O distances characterizing the 8MR in LTA zeolite (termed W1 and W2) for representative trajectories of (a) methane, (b) ethane, (c) propane, and (d) butane. Data were recorded for a trajectory corresponding to the hopping of the molecule between two adjacent cages. The lower panel illustrates the variation of the order parameter along the trajectory. 15648


Article

The Journal of Physical Chemistry C In general, the diffusion coefficients of n-alkanes in zeolites decrease as the chain length increases.5 When the kinetic diameter of a hydrocarbon molecule is comparable to the pore size, the diffusion process is strongly controlled by the pore size. In this case, vibrations of framework atoms defining the window that the molecule goes through play a crucial role in the diffusion. Earlier studies have well demonstrated the limitations of the use of a rigid framework for computing diffusion in this regime.23,29 The magnitude of the framework flexibility effect on diffusion is connected to the harmonicity of the vibrations of the atoms forming the window and also to the number of distinct directions used to define its minimum size. Using a statistically representative series of configurations of an empty adsorbent framework, these effects on the diffusion of spherical molecules were accurately captured.23 While this helped to gain valuable insights into molecular diffusion, the relationship between adsorbates with internal degrees of freedom and framework flexibility remains less clear. To understand factors controlling the decrease in selfdiffusivities and also shed light on possible correlations between framework flexibility and intramolecular degrees of freedom, we used a series of four hydrocarbons with increasing kinetic diameter,75 CH4 (3.8 Å), C2H6 (3.9 Å), C3H8 (4.3 Å), and C4H10 (4.3 Å), to investigate the use of TPS simulations. Combined with NpT molecular dynamics, all degrees of freedom of both the adsorbate and adsorbent have been accounted for. To solely focus on the adsorbate−adsorbent interactions, we performed our simulations at infinite dilution. This also allowed us to compare our results for small molecules with those of other methods that approximate framework flexibility. First, we started from the simplest case, methane, which is approximated as a spherical adsorbate. Our computed diffusion coefficient (see Table 1) is in good agreement with both standard MD results and the value reported by Awati et al.,23 who used the same interaction models and force field parameters. Garcia-Sanchez et al.19 compared the results obtained using the Hill and Sauer52,53 and Nicholas et al.76 force fields for methane adsorption and diffusion in LTA zeolites. They reported diffusivities on the order of 10−9−10−10 m2·s−1 for pure silica LTA at 500 K. Ethane has one more degree of freedom than methane (in the united atom representation). This change in adsorbate shape from spherical to linear resulted in a decrease by an order of magnitude for diffusivity (see Table 1). Unlike methane, ethane has a preferred orientation for crossing the 8MR window. Again, our TPS results are in good agreement with standard MD simulations. Another decrease in diffusivity by an order of magnitude is observed in changing from ethane to propane (see Table 1). The loss of adsorbate linearity from ethane to propane is accompanied by an increase in the kinetic diameter. The analysis of the set of successful trajectories collected during the sampling indicated longer residence times in the window when compared to those of CH4. The lower panels in Figure 7 illustrate the variation of the order parameter (defined above) used for TPS simulations. Unlike C2H6 and C3H8, CH4 required less than 1000 fs to cross from one cage to the adjacent one through the 8MR window. Considering that the typical breathing frequency for pure silica LTA is around 150 cm−1,23 we see that the window underwent on the order of five vibrations (size changes) while the methane was crossing over. Longer chains sample more pore size variations during a

crossing event (e.g., Figure 7d). This result highlights the importance of framework flexibility for diffusion simulations of extended molecules passing through narrow pores. We also examined the situation of adsorbates with similar degrees of freedom and shapes but with slightly different kinetic diameters. In Table 1, we compared the diffusivities of ethane and propane to those of ethylene and propylene. Although the computed diffusivities of ethane and ethylene are almost identical, propylene was around 4 times faster than propane. This different behavior cannot be explained by differences in the kinetic diameter alone. Below, an analysis of the window size changes and residence times at the window will shed more light on this point. The diffusivities computed in this work show a trend different from that of the results reported in experiments on diffusion in the cationic zeolites NaCaA12,77,78 and CaA12 and in cation-free AlPO-LTA.79 This difference in trends is unsurprising because of the considerable influences that cations have on adsorbate diffusion in these materials. Self-diffusivities reported by Hedin et al.11,12 from pulsed field gradient (PFG) NMR experiments in Si-LTA are faster than our calculated values by 1 order of magnitude for methane and ethane/ethylene and slower by 2 orders of magnitude for propylene. Using X-ray diffraction, they determined the cell parameter, volume, and window size of the Si-LTA structure to be 11.84 Å, 1660.0 Å3, and 4.00 Å, respectively. The Hill−Sauer force field used in this work predicts a slightly larger cell parameter and volume, 12.037 Å (1.66% larger than that from experiment) and 1744.0 Å3 (5% larger than that from experiment), respectively. The average window size from our NpT simulations is 3.78 Å. This is smaller than the size reported in the experiments of Hedin et al.11 (4.00 Å) and Corma et al.38 (3.995 Å). This difference, although small in absolute magnitude, may have a considerable influence on molecular diffusion. If our calculations underestimated the 8MR window size relative to that of the real material, this would lead to calculated diffusivities that are lower than the experimental values. This explanation appears plausible for the data from Hedin et al.11 for methane, ethane, and ethylene. The same simple explanation cannot be applied to propylene, for which the diffusivity measured by Hedin et al. is much slower (4.7 × 10−15 m2·s−1) than our computed infinite dilution diffusivity (9.3 × 10−13 m2·s−1) Previous studies44,80−83 have shown that the diffusivity of small adsorbates in cage-type zeolites typically increases as the loading is increased for low and moderate loadings before decreasing at high loadings due to steric effects. Hedin et al. conducted their measurements for propylene at a loading of 3.4 molecules/cage,11 which is almost the saturation loading of 3.6 molecules/cage as determined by Corma et al.38 at 298 K. It therefore seems likely that the infinite dilution diffusivity associated with these experimental data is higher than Hedin et al.’s finite loading measurement. This lessens the discrepancy between our predicted diffusivity and experimental observation, but it still appears that our computed value corresponds to considerably faster diffusion of propylene in SiLTA than has been seen experimentally. In the context of the comparison with the experiment above, it is important to note that our calculations used force field parameters for framework−adsorbate interactions58,62 in combination with TraPPE50 and Hill−Sauer52,53 force fields without refitting. Simulation methods to assess the diffusivity of molecules such as ethylene and propylene in flexible but tightfitting 8MR zeolites were not previously available. Our 15649


Article


Figure 8. Histograms of the LTA 8MR window size distribution for (a) methane, (b) ethane, (c) propane, and (d) butane. For comparison, the histogram corresponding to an empty LTA framework is shown in blue. Data were collected from a configuration corresponding to transition states along the transition of a molecule from cage A to cage B. Each histogram is normalized to have the same area under the curve.

kinetic diameter (4.3 and 3.9 Å for C3H8 and C3H6, respectively) has a greater effect on diffusion than in the case of ethane/ethylene because of the magnitude of the window size changes during the hopping, which reflect stronger adsorbate−adsorbent interactions. Another important observation from Figure 7 is the appearance of an order parameter plateau for n-butane (Figure 7d). This indicates a longer crossing time spent at the window (∼2 ps) even when compared to that of propane (∼1 ps). To achieve a statistically better description of the qualitative information obtained above from single trajectories, the set of transition-state ensembles corresponding to each adsorbate has been determined using a systematic approach based on the committor analysis described in section 3.7. This approach consists of monitoring the 8MR window size behavior at times of molecular hopping. This requires identification of possible transition states for each set of trajectory ensembles collected during TPS simulations. Then the fluctuations of the window size are recorded and analyzed. The collected data are represented in the form of histograms (Figure 8) and compared to the window size distribution for an empty LTA framework (Figure 8, blue curve) obtained from a 100 ns NpT MD run. Histograms corresponding to methane and ethane confirm that the flexibility of the framework can be easily approximated using a series of snapshots collected from

calculations demonstrate that relatively direct comparisons between computed and measured diffusivities in these systems can now be made. It will be useful in the future to extend these kinds of calculations to explore the variation among predicted diffusivities as a function of the force fields used to describe adsorbates and zeolite frameworks. Framework−adsorbate interactions such as those we have used here are typically developed to accurately describe adsorption isotherms. While this is a necessary attribute of an accurate force field, this approach focuses attention on the favorable energy states in a material, and as a result, there is no reason to expect that these force fields will automatically describe activated diffusion accurately.82 To better capture the effect of framework flexibility on diffusion while increasing the chain length, representative trajectories for each molecule were selected from the path ensemble after sampling convergence and the corresponding window size fluctuations were recorded along the transition (Figure 7). It is easy to distinguish two different patterns of window size variation from Figure 7. The first difference is the tendency of the longer adsorbates to spend more time to achieve a single hop between adjacent cages. Second, an increase of window size during the hopping of propane and butane is observed. This helps explain the difference between the diffusivities of propane and propylene. The difference in the 15650



■

the empty framework.23,29 No coupling or correlation between adsorbate hops and window size fluctuations is evident during the diffusion process. The decrease in diffusivity of ethane is mainly due to a larger kinetic diameter and the linear geometry that requires the adsorbate to be somewhat perpendicular to the plane defining the 8MR window, unlike the spherically shaped methane molecule. The appearance of a tail in the window size distribution at higher values of the window size is observed for propane and butane in parts c and d, respectively, of Figure 8. This is compelling evidence of the coupling between the adsorbates and the zeolite framework as molecules pass through the 8MR windows. Simplistically, the molecule can be thought of as “holding the pore open” as it passes through. Simulation methods that rely on the assumption that the framework dynamics is decoupled from the dynamics of the diffusing molecule cannot accurately predict hopping rates in this situation. Our results show that TPS is a useful approach for characterizing molecular diffusion, especially in situations where there is a significant influence of the adsorbate on the framework dynamics as molecules pass through transition states, as in the case of propane or butane in LTA.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Phone: +1-404-8942822. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Vermeiren, W.; Gilson, J. P. Impact of Zeolites on the Petroleum and Petrochemical Industry. Top. Catal. 2009, 52, 1131−1161. (2) Song, L.; Sun, Z.; Duan, L.; Gui, H.; McDougall, G. S. Adsorption and Diffusion Properties of Hydrocarbons. Microporous Mesoporous Mater. 2007, 104, 115−128. (3) Li, Z. J.; Johnson, M. C.; Sun, M. W.; Ryan, E. T.; Earl, D. J.; Maichen, W.; Martin, J. I.; Li, S.; Lew, C. M.; Wang, J.; et al. Mechanical and Dielectric Properties of Pure-Silica-Zeolite Low-K Materials. Angew. Chem., Int. Ed. 2006, 45, 6329−6332. (4) Corma, A. State of the Art and Future Challenges of Zeolites as Catalysts. J. Catal. 2003, 216, 298−312. (5) Smit, B.; Maesen, T. L. M. Molecular Simulations of Zeolites: Adsorption, Diffusion, and Shape Selectivity. Chem. Rev. 2008, 108, 4125−4184. (6) Sun, J. L.; Bonneau, C.; Cantin, A.; Corma, A.; Diaz-Cabanas, M. J.; Moliner, M.; Zhang, D. L.; Li, M. R.; Zou, X. D. The ITQ-37 Mesoporous Chiral Zeolite. Nature 2009, 458, 1154−U90. (7) Camblor, M. A.; Corma, A.; Lightfoot, P.; Villaescusa, L. A.; Wright, P. A. Synthesis and Structure of ITQ-3, the First Pure Silica Polymorph with a Two-Dimensional System of Straight Eight-Ring Channels. Angew. Chem., Int. Ed. Engl. 1997, 36, 2659−2661. (8) Barrett, P. A.; Boix, T.; Puche, M.; Olson, D. H.; Jordan, E.; Koller, H.; Camblor, M. A. ITQ-12: A New Microporous Silica Polymorph Potentially Useful for Light Hydrocarbon Separations. Chem. Commun. 2003, 2114−2115. (9) Olson, D. H.; Camblor, M. A.; Villaescusa, L. A.; Kuehl, G. H. Light Hydrocarbon Sorption Properties of Pure Silica Si-CHA and ITQ-3 and High Silica ZSM-58. Microporous Mesoporous Mater. 2004, 67, 27−33. (10) Palomino, M.; Cantin, A.; Corma, A.; Leiva, S.; Rey, F.; Valencia, S. Pure Silica ITQ-32 Zeolite Allows Separation of Linear Olefins from Paraffins. Chem. Commun. 2007, 1233−1235. (11) Hedin, N.; DeMartin, G. J.; Roth, W. J.; Strohmaier, K. G.; Reyes, S. C. PFG NMR Self-Diffusion of Small Hydrocarbons in High Silica DDR, CHA and LTA Structures. Microporous Mesoporous Mater. 2008, 109, 327−334. (12) Hedin, N.; DeMartin, G. J.; Strohmaier, K. G.; Reyes, S. C. PFG NMR Self-Diffusion of Propylene in ITQ-29, CaA and NaCaA: Window Size and Cation Effects. Microporous Mesoporous Mater. 2007, 98, 182−188. (13) Cantin, A.; Corma, A.; Leiva, S.; Rey, F.; Rius, J.; Valencia, S. Synthesis and Structure of the Bidimensional Zeolite ITQ-32 with Small and Large Pores. J. Am. Chem. Soc. 2005, 127, 11560−11561. (14) Krishna, R. Describing the Diffusion of Guest Molecules inside Porous Structures. J. Phys. Chem. C 2009, 113, 19756−19781. (15) Sholl, D. S. Understanding Macroscopic Diffusion of Adsorbed Molecules in Crystalline Nanoporous Materials via Atomistic Simulations. Acc. Chem. Res. 2006, 39, 403−411. (16) Dubbeldam, D.; Snurr, R. Q. Recent Developments in the Molecular Modeling of Diffusion in Nanoporous Materials. Mol. Simul. 2007, 33, 305−325. (17) Smit, B.; Siepmann, J. I. Computer-Simulations of the Energetics and Siting of N-Alkanes in Zeolites. J. Phys. Chem. 1994, 98, 8442−8452. (18) June, R. L.; Bell, A. T.; Theodorou, D. N. Molecular-Dynamics Studies of Butane and Hexane in Silicalite. J. Phys. Chem. 1992, 96, 1051−1060. (19) Garcia-Sanchez, A.; Dubbeldam, D.; Calero, S. Modeling Adsorption and Self-Diffusion of Methane in LTA Zeolites: The

5. SUMMARY In this paper, we have examined the role of the adsorbate and framework flexibility on the diffusion of linear alkanes in pure silica LTA at room temperature. While the decrease in diffusivity of n-alkanes as a function of the chain length is easy to predict qualitatively, accurate computational characterization of this diffusion is hindered by the time-scale problem, especially for nonspherical and nonlinear adsorbates. Using TPS molecular dynamics simulations, we computed the diffusivities for methane, ethane, ethylene, propane, propylene, and butane. The simulation procedure used in this work is based on fully flexible NpT MD trajectories, which enabled us to collect unbiased information about window size changes during the crossing of different adsorbates. The analysis of the corresponding configurations revealed that for small molecules (methane, ethane, and ethylene) the framework flexibility could be accurately approximated using snapshots from a long run of molecular dynamics for an empty LTA structure. Propane, propylene, and n-butane, however, exhibited a longer residence time at the window. These longer residence times cause significant coupling between the adsorbate and framework, making approximate methods that assume that these degrees of freedom are decoupled inapplicable. This is likely to be a generic situation when extended molecules diffuse through highly constricted pores in zeolites and other nanoporous materials. The TPS method we have used is a useful member of the set of simulation tools that are available for computing molecular diffusivities in these materials. TPS is particularly well suited to situations such as butane in LTA, where including full coupling between the adsorbate and framework is required to achieve an accurate description while such a slow diffusion cannot be readily treated with standard MD simulations.

■

Article

ASSOCIATED CONTENT

* Supporting Information S

Detailed description of the force field models and the interaction parameters used in this work. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.5b01633. 15651


Article

The Journal of Physical Chemistry C Influence of Framework Flexibility. J. Phys. Chem. C 2010, 114, 15068−15074. (20) Leroy, F.; Rousseau, B. Self-Diffusion of n-Alkanes in MFI Type Zeolite Using Molecular Dynamics Simulations with an Anisotropic United Atom (AUA) Forcefield. Mol. Simul. 2004, 30, 617−620. (21) Krishna, R.; van Baten, J. M. A Molecular Dynamics Investigation of the Diffusion Characteristics of Cavity-Type Zeolites with 8-Ring Windows. Microporous Mesoporous Mater. 2011, 137, 83− 91. (22) Krishna, R.; van Baten, J. M. Comment on “Modeling Adsorption and Self-Diffusion of Methane in Lta Zeolites: The Influence of Framework Flexibility”. J. Phys. Chem. C 2010, 114, 18017−18021. (23) Awati, R. V.; Ravikovitch, P. I.; Sholl, D. S. Efficient and Accurate Methods for Characterizing Effects of Framework Flexibility on Molecular Diffusion in Zeolites: CH4 Diffusion in Eight Member Ring Zeolites. J. Phys. Chem. C 2013, 117, 13462−13473. (24) Eyring, H. The Activated Complex in Chemical Reactions. J. Chem. Phys. 1935, 3, 107−115. (25) Chandler, D. Statistical-Mechanics of Isomerization Dynamics in Liquids and Transition-State Approximation. J. Chem. Phys. 1978, 68, 2959−2970. (26) Dubbeldam, D.; Smit, B. Computer Simulation of Incommensurate Diffusion in Zeolites: Understanding Window Effects. J. Phys. Chem. B 2003, 107, 12138−12152. (27) Dubbeldam, D.; Beerdsen, E.; Calero, S.; Smit, B. Dynamically Corrected Transition State Theory Calculations of Self-Diffusion in Anisotropic Nanoporous Materials. J. Phys. Chem. B 2006, 110, 3164− 3172. (28) Abouelnasr, M. K. F.; Smit, B. Diffusion in Confinement: Kinetic Simulations of Self- and Collective Diffusion Behavior of Adsorbed Gases. Phys. Chem. Chem. Phys. 2012, 14, 11600−11609. (29) Haldoupis, E.; Watanabe, T.; Nair, S.; Sholl, D. S. Quantifying Large Effects of Framework Flexibility on Diffusion in MOFs: CH4 and CO2 in ZIF-8. ChemPhysChem 2012, 13, 3449−3452. (30) Combariza, A. F.; Sastre, G.; Corma, A. Molecular Dynamics Simulations of the Diffusion of Small Chain Hydrocarbons in 8-Ring Zeolites. J. Phys. Chem. C 2011, 115, 875−884. (31) Dellago, C.; Bolhuis, P. G.; Csajka, F. S.; Chandler, D. Transition Path Sampling and the Calculation of Rate Constants. J. Chem. Phys. 1998, 108, 1964−1977. (32) Bolhuis, P. G.; Dellago, C.; Chandler, D. Sampling Ensembles of Deterministic Transition Pathways. Faraday Discuss. 1998, 110, 421− 436. (33) Dellago, C.; Bolhuis, P. G.; Chandler, D. On the Calculation of Reaction Rate Constants in the Transition Path Ensemble. J. Chem. Phys. 1999, 110, 6617−6625. (34) van Erp, T. S.; Moroni, D.; Bolhuis, P. G. A Novel Path Sampling Method for the Calculation of Rate Constants. J. Chem. Phys. 2003, 118, 7762−7774. (35) van Erp, T. S.; Bolhuis, P. G. Elaborating Transition Interface Sampling Methods. J. Comput. Phys. 2005, 205, 157−181. (36) Vlugt, T. J. H.; Dellago, C.; Smit, B. Diffusion of Isobutane in Silicalite Studied by Transition Path Sampling. J. Chem. Phys. 2000, 113, 8791−8799. (37) Bolhuis, P. G.; Chandler, D.; Dellago, C.; Geissler, P. L. Transition Path Sampling: Throwing Ropes over Rough Mountain Passes, in the Dark. Annu. Rev. Phys. Chem. 2002, 53, 291−318. (38) Corma, A.; Rey, F.; Rius, J.; Sabater, M. J.; Valencia, S. Supramolecular Self-Assembled Molecules as Organic Directing Agent for Synthesis of Zeolites. Nature 2004, 431, 287−290. (39) Nagumo, R.; Takaba, H.; Nakao, S. I. High-Accuracy Estimation of ‘Slow’ Molecular Diffusion Rates in Zeolite Nanopores, Based on Free Energy Calculations at an Ultrahigh Temperature. J. Phys. Chem. C 2008, 112, 2805−2811. (40) Beerdsen, E.; Smit, B.; Dubbeldam, D., Molecular Simulation of Loading Dependent Slow Diffusion in Confined Systems. Phys. Rev. Lett. 2004, 93.

(41) Astala, R.; Auerbach, S. M.; Monson, P. A. Density Functional Theory Study of Silica Zeolite Structures: Stabilities and Mechanical Properties of SOD, LTA, CHA, MOR, and MFI. J. Phys. Chem. B 2004, 108, 9208−9215. (42) Fritzsche, S.; Haberlandt, R.; Hofmann, G.; Karger, J.; Heinzinger, K.; Wolfsberg, M. An MD Study on the Diffusion of Methane in a Cation-Free LTA Zeolite. Illustrations and New Results. Chem. Phys. Lett. 1997, 265, 253−258. (43) Haberlandt, R. Transport Processes in Porous Media: Diffusion in Zeolites. Thin Solid Films 1998, 330, 34−45. (44) Demontis, P.; Suffritti, G. B. Sorbate-Loading Dependence of Diffusion Mechanism in a Cubic Symmetry Zeolite of Type ZK4. A Molecular Dynamics Study. J. Phys. Chem. B 1997, 101, 5789−5793. (45) Dellago, C.; Bolhuis, P. G.; Chandler, D. Efficient Transition Path Sampling: Application to Lennard-Jones Cluster Rearrangements. J. Chem. Phys. 1998, 108, 9236−9245. (46) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press: Oxford, U.K., 1987. (47) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, CA, 2001. (48) Ding, K. J.; Valleau, J. P. Umbrella-Sampling Realization of Widom Chemical-Potential Estimation. J. Chem. Phys. 1993, 98, 3306−3312. (49) Dellago, C.; Bolhuis, P. G.; Geissler, P. L. Transition Path Sampling. Adv. Chem. Phys. 2002, 123, 1−78. (50) Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 1. United-Atom Description of n-Alkanes. J. Phys. Chem. B 1998, 102, 2569−2577. (51) Wick, C. D.; Martin, M. G.; Siepmann, J. I. Transferable Potentials for Phase Equilibria. 4. United-Atom Description of Linear and Branched Alkenes and Alkylbenzenes. J. Phys. Chem. B 2000, 104, 8008−8016. (52) Hill, J. R.; Sauer, J. Molecular Mechanics Potential for Silica and Zeolite Catalysts Based on Ab-Initio Calculations. 1. Dense and Microporous Silica. J. Phys. Chem. 1994, 98, 1238−1244. (53) Hill, J. R.; Sauer, J. Molecular Mechanics Potential for Silica and Zeolite Catalysts Based on Ab-Initio Calculations. 2. Aluminosilicates. J. Phys. Chem. 1995, 99, 9536−9550. (54) Vlugt, T. J. H.; Krishna, R.; Smit, B. Molecular Simulations of Adsorption Isotherms for Linear and Branched Alkanes and Their Mixtures in Silicalite. J. Phys. Chem. B 1999, 103, 1102−1118. (55) June, R. L.; Bell, A. T.; Theodorou, D. N. Molecular Dynamics Studies of Butane and Hexane in Silicalite. J. Phys. Chem. 1992, 96, 1051−1060. (56) Smit, B.; Daniel, J. C.; Loyens, L.; L. M, M.; Verbist, G. Simulation of Adsorption and Diffusion of Hydrocarbons in Zeolites. Faraday Discuss. 1997, 106, 93−104. (57) Macedonia, M. D.; Maginn, E. J. A Biased Grand Canonical Monte Carlo Method for Simulating Adsorption Using All-Atom and Branched United Atom Models. Mol. Phys. 1999, 96, 1375−1390. (58) Calero, S.; Dubbeldam, D.; Krishna, R.; Smit, B.; Vlugt, T. J. H.; Denayer, J. F. M.; Martens, J. A.; Maesen, T. L. M. Understanding the Role of Sodium during Adsorption: A Force Field for Alkanes in Sodium-Exchanged Faujasites. J. Am. Chem. Soc. 2004, 126, 11377− 11386. (59) Dubbeldam, D.; Calero, S.; Vlugt, T. J. H.; Krishna, R.; Maesen, T. L. M.; Beerdsen, E.; Smit, B. Force Field Parametrization through Fitting on Inflection Points in Isotherms. Phys. Rev. Lett. 2004, 93, 088302. (60) Bai, P.; Tsapatsis, M.; Siepmann, J. I. TraPPE-zeo: Transferable Potentials for Phase Equilibria Force Field for All-Silica Zeolites. J. Phys. Chem. C 2013, 117, 24375−24387. (61) Pascual, P.; Ungerer, P.; Tavitian, B.; Pernot, P.; Boutin, A. Development of a Transferable Guest-Host Force Field for Adsorption of Hydrocarbons in Zeolites I. Reinvestigation of Alkane Adsorption in Silicalite by Grand Canonical Monte Carlo Simulation. Phys. Chem. Chem. Phys. 2003, 5, 3684−3693. 15652


Article

The Journal of Physical Chemistry C (62) Dubbeldam, D.; Calero, S.; Vlugt, T. J. H.; Krishna, R.; Maesen, T. L. M.; Smit, B. United Atom Force Field for Alkanes in Nanoporous Materials. J. Phys. Chem. B 2004, 108, 12301−12313. (63) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular-Dynamics. J. Comput. Phys. 1995, 117, 1−19. (64) Shinoda, W.; Shiga, M.; Mikami, M., Rapid Estimation of Elastic Constants by Molecular Dynamics Simulation under Constant Stress. Phys. Rev. B 2004, 69. (65) Martyna, G. J.; Tobias, D. J.; Klein, M. L. Constant-Pressure Molecular-Dynamics Algorithms. J. Chem. Phys. 1994, 101, 4177− 4189. (66) Parrinello, M.; Rahman, A. Polymorphic Transitions in SingleCrystalsA New Molecular-Dynamics Method. J. Appl. Phys. 1981, 52, 7182−7190. (67) Hoover, W. G. Canonical DynamicsEquilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (68) Hoover, W. G. Constant-Pressure Equations of Motion. Phys. Rev. A 1986, 34, 2499−2500. (69) Melchionna, S.; Ciccotti, G.; Holian, B. L. Hoover NPT Dynamics for Systems Varying in Shape and Size. Mol. Phys. 1993, 78, 533−544. (70) Lazaridis, T.; Karplus, M. “New View” of Protein Folding Reconciled with the Old through Multiple Unfolding Simulations. Science 1997, 278, 1928−1931. (71) Boulfelfel, S. E.; Zahn, D.; Hochrein, O.; Grin, Y.; Leoni, S. Low-Dimensional Sublattice Melting by Pressure: Superionic Conduction in the Phase Interfaces of the Fluorite-to-Cotunnite Transition of CaF2. Phys. Rev. B 2006, 74, 094106. (72) Boulfelfel, S. E.; Seifert, G.; Grin, Y.; Leoni, S. Squeezing Lone Pairs: The A17 to A7 Pressure-Induced Phase Transition in Black Phosphorus. Phys. Rev. B 2012, 85, 014110. (73) Boulfelfel, S. E.; Oganov, A. R.; Leoni, S. Understanding the Nature of “Superhard Graphite”. Sci. Rep. 2012, 2, No. 471. (74) Leoni, S.; Boulfelfel, S. E. Pathways of Structural Transformations in Reconstructive Phase Transitions: Insights from Transition Path Sampling Molecular Dynamics. Modern Methods of Crystal Structure Prediction; Wiley-VCH Verlag GmbH & Co. KGaA: Weinheim, Germany, 2010; pp 181−221. (75) Matteucci, S.; Yampolskii, Y.; Freeman, B. D.; Pinnau, I. Transport of Gases and Vapors in Glassy and Rubbery Polymers. Materials Science of Membranes for Gas and Vapor Separation; John Wiley & Sons, Ltd.: Chichester, U.K., 2006; pp 1−47. (76) Nicholas, J. B.; Hopfinger, A. J.; Trouw, F. R.; Iton, L. E. Molecular Modeling of Zeolite Structure. 2. Structure and Dynamics of Silica Sodalite and Silicate Force Field. J. Am. Chem. Soc. 1991, 113, 4792−4800. (77) Heink, W.; Karger, J.; Pfeifer, H.; Salverda, P.; Datema, K. P.; Nowak, A. Self-Diffusion Measurements of n-Alkanes in Zeolite NaCaA by Pulsed-Field Gradient Nuclear Magnetic Resonance. J. Chem. Soc., Faraday Trans. 1992, 88, 515−519. (78) Karger, J. Comment on “PFG NMR Self-Diffusion of Small Hydrocarbons in High Silica DDR, CHA and LTA Structures” [Micropor. Mesopor. Mater. 109 (2008) 327]. Microporous Mesoporous Mater. 2008, 116, 715−717. (79) Hibbe, F.; Caro, J.; Chmelik, C.; Huang, A. S.; Kirchner, T.; Ruthven, D.; Valiullin, R.; Karger, J. Monitoring Molecular Mass Transfer in Cation-Free Nanoporous Host Crystals of Type AlPOLTA. J. Am. Chem. Soc. 2012, 134, 7725−7732. (80) Beerdsen, E.; Dubbeldam, D.; Smit, B. Loading Dependence of the Diffusion Coefficient of Methane in Nanoporous Materials. J. Phys. Chem. B 2006, 110, 22754−22772. (81) van den Bergh, J.; Ban, S.; Vlugt, T. J. H.; Kapteijn, F. Modeling the Loading Dependency of Diffusion in Zeolites: The Relevant Site Model. J. Phys. Chem. C 2009, 113, 17840−17850. (82) Jee, S. E.; Sholl, D. S. Carbon Dioxide and Methane Transport in DDR Zeolite: Insights from Molecular Simulations into Carbon Dioxide Separations in Small Pore Zeolites. J. Am. Chem. Soc. 2009, 131, 7896−7904.

(83) Skoulidas, A. I.; Sholl, D. S. Molecular Dynamics Simulations of Self-Diffusivities, Corrected Diffusivities, and Transport Diffusivities of Light Gases in Four Silica Zeolites To Assess Influences of Pore Shape and Connectivity. J. Phys. Chem. A 2003, 107, 10132−10141.

15653


Modeling Diffusion of Linear Hydrocarbons in Silica Zeolite LTA Using

Recommend Documents