Article pubs.acs.org/JPCC
Role of Ion Mobility in Molecular Sieving of CO2 over N2 with Zeolite NaKA Amber Mace, Niklas Hedin, and Aatto Laaksonen* Department of Materials and Environmental Chemistry and Berzelii Centre EXSELENT on Porous Materials, Stockholm University, Stockholm, Sweden S Supporting Information *
ABSTRACT: Classical molecular dynamics and Grand Canonical Monte Carlo simulations are carried out for sorbates (CO2 and N2) in zeolite NaKA using a universal type ab initio force field. By combining the results of these simulations, we reproduce the CO2 uptake as a function of the K+ content in the zeolite NaKA as measured experimentally by Liu et al.1 The experiment yielded an exceptionally high CO 2 -over-N 2 selectivity of >200 at a specific K+/(K+ + Na+) ratio of 17 atom %. This high selectivity could be attributed to the nonlinear uptake dependency of the K+/(K+ + Na+) ratio measured for both CO2 and N2. Additionally, our simulations show a strong coupling between the self-diffusion of CO2 and the site-to-site jumping rate of the extra-framework cations. These results show that this nonlinear uptake dependency of CO2 is the result of molecular sieving. Following this, our simulations conclude that both thermodynamic and kinetic contributions must be taken into account when modeling the uptake of this and similar materials with the same functionalities.
■
INTRODUCTION Adsorption-driven separation of CO2 from N2-rich flue gas could potentially lower the cost for CO2 capture as compared with aqueous solutions of amine-based absorption/reaction processes.2 Most adsorption studies have revolved around equilibrium based separation principles. It appears that an intermediate value of the isosteric heat of adsorption gives the lowest parasitic energy for a generic separation process, which was defined by Lin et al.3 as the minimum electric load imposed on a power plant by a temperature−pressure swing capture process. This study relates to the possiblility that parasitic energy could further be lowered by employing kinetic aspects to CO2-over-N2 selection by adsorption. Zeolite A (framework type code LTA) is a three-dimensional, cage-type microporous aluminosilicate with an Al-to-Si ratio of 1 resulting in a comparatively high cation content acting to neutralize the framework. The zeolite A framework with point group Fm3c has a cubic periodicity built up by so-called α- and β-cages. As illustrated in Figure 1, one unit cell consists of eight (2 × 2 × 2) α-cages, and the smaller β-cage is formed in the cubic center of the conformation, where each unit cell comprises a total of five β-cages. The α-cages intersect with 8ring (8R) pore windows, the β-cages with 4-rings (4R), and the α-cages intersect the β-cages with 6-rings (6R). The number of oxygen atoms within the ring, confining the pore window, determines the window ring number. It is only the α-cages that can host CO2 or N2 molecules, as the 4R and 6R pore windows intersecting the β-cages are too narrow to pass. In the unit cell © XXXX American Chemical Society
Figure 1. LTA framework is built up by the larger α-cages (bottom right) and the smaller β-cages (top right). One LTA unit cell (left) consists of eight α-cages, which in turn construct the five β-cages in their cubic intersections.
of zeolite A, the cations are distributed over three types of sites (I, II, and III).1,4 Site I is centered with the 6-ring, site II is in the plane of the 8R, and site III is centered with the 4R. These sites will be referred to as the 6R, 8R, and 4R site, respectively. In general, each of the 24 8R and 64 6R per unit cell have one Received: October 8, 2012 Revised: October 22, 2013
A
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Several articles have been published on the simulation of equilibrium uptakes of CO2 and N2 gases in cation containing zeolites using GCMC.15−20 Further, there are studies, which predict the diffusion of CO2 within the framework of 8-ring zeolites,21−23 as well as the mobility of the extra-framework cations in zeolites upon gas adsorption,24−26 with MD simulations. These studies on extra-framework ion behavior showed that the mobility of the Na+ ions in the zeolite was, to some extent, dependent on the loading of the respective adsorbate. However, the current work is, to the best of our knowledge, first of its kind in the aspect of showing a connection between the diffusion rate of the adsorbed gas molecules and the mobility of the extra-framework cations and how this can be tuned through ion-exchange. We combine the quantitative results from equilibrium GCMC computations with those from MD diffusion simulations in order to get a deep view of the dynamics involved in the high CO2-over-N2 selectivity of zeolite NaKA observed by Liu et al.1 We will show how the site-to-site jumping rates of the extra-framework cations correlates with the amount of adsorbed gas molecules measured experimentally. It appears that the sieving mechanism in zeolite A has actually much more dynamic character than earlier recognized.
cation coordinated with it, while the eight remaining cations will be coordinated with a 4R distributed throughout the unit cell. As it is only the 8R that is large enough to allow the passage of CO2 or N2, it is the cations occupying site II, which regulates whether the passage of CO2 or N2 is possible. The material with Na+ is referred to as zeolite NaA or 4A and with K+, zeolite KA or 3A, where 4A and 3A denotes the resulting pore window diameters of approximately 4 and 3 Å, respectively. Zeolite NaKA refers to a mixed Na+ and K+ extra-framework distribution. Recently, Liu et al.1 communicated particularly appealing results showing that through K+-exchange in zeolite NaA, an exceptionally high CO2-over-N2 selectivity could be reached, while the total CO2 uptake still remained high. They proposed that this high selectivity was due to a sieving effect. An intermediate K+/(K+ + Na+) ratio of 17 atom % resulted in an effective pore diameter lying right in between the kinetic diameters of CO2 and N2, which are 3.3 Å and 3.6 Å respectively,5 allowing the material to adsorb CO2, while effectively blocking out N2. For CO2-over-N2 selectivity, this phenomenon has been scarcely investigated; however, two important works are those done by Breck et al. in 1956 and Yeh and Yang in 1989. Breck et al. showed experimentally that the CO2 uptake was dependent on the K+/(K+ + Na+) ratio in the ion-exchanged zeolite A.6 Yeh and Yang rationalized the sorbent diffusivity with a theoretical approach on poreblocking.7 If the sorbent shows a sieving effect, which can be tuned through ion exchange, both thermodynamic and kinetic effects are suggested to contribute to the dynamics involved in the separation process. Here, we attempt to simulate these two contributions to the adsorption dynamics using an ab initio based force field to achieve qualitatively robust results. The main methods used for such investigation are Grand Canonical Monte Carlo (GCMC) simulations and Molecular Dynamics (MD). GCMC is commonly used when simulating the equilibrium level of the gas uptake in porous sorbates and is based purely on the descriptors of statistical thermodynamics. We use such GCMC simulations to estimate the thermodynamic contribution to the adsorptive selection of CO2-over-N2 on zeolite NaKA with different K+ levels. To be able to compare with the experiments by Liu et al.,1 we defined selectivity conservatively to the ratio of the uptake of CO2 and N2 as single component gases. We leave mixed gas simulations for future work. The actual selectivity is expected to be much higher than this ratio according to ideal adsorbed solution theory (IAST).8 The GCMC results can be compared with selectivities in similar aluminum rich zeolites without size restricting pores.9−14 Such zeolites exhibit single component selectivities generally in the range of 10 as a result of CO2 possessing a quadrupolar moment approximately three times that of N2. As we also are interested in how the gas molecules move through zeolite NaKA over time, we performed MD simulations to obtain the self-diffusion rate of the gas molecules through the material and, in turn, its sieving functionality. By combining the results from these two methods, we attempt to explain the experimental data and to further reveal new information on the detailed mechanisms of the sieving functionality. This information may be a key in the investigation and development of new adsorbents with a reduced parasitic energy for CO2-over-N2 separation from flue gas.
■
COMPUTATIONAL METHODS All software used throughout this work is a part of the Materials Studio software suite (version 6.0, Accelrys Inc.)27 The all-atom COMPASS-force field28 is used in all simulations performed in this work. It is of universal type, derived by using ab initio data parametrized from a variety of molecular classes and extensively validated using MD and energy minimizing procedures. The COMPASS-force field has a 9−6 functional form of LennardJones-type and utilizes the Lorentz−Berthelot combination rules of sixth order29 for nonbonded interactions. Covalent bonds in the CO2 and N2 molecules are described by valence terms allowing bond-stretching and bending (for CO2).30,31 For the zeolite framework, the semi-ionic model32 is used where the framework interatomic interaction is modeled with a Morsedispersion function, while a simple ionic model33 is used for the extra-framework cations. Hence, a fully flexible system is permitted. More details regarding the parameters and functional forms used by the COMPASS force field can be found in Supporting Information (SI). We chose to use the COMPASS force field despite the fact that numerous empirically derived force fields have been developed specifically for zeolite A, which accurately reproduces experimental adsorption isotherms.15,18−20 As the objective of this work is to predict the extent, to which the uptake is facilitated/hindered by thermodynamics and kinetics, respectively, we are looking to avoid any bias which may be the case for a force field which is empirically parametrized to describe the gas adsorption as solely a thermodynamic effect. The framework coordinates for zeolite NaKA were taken from the Materials Studio structure database as originally presented by Gramlich and Meier.34 In the composition (Al96Si96O394), the framework-type LTA has a cubic unit cell in the crystallographic space group Fm3c with cell dimensions a = b = c = 24.61 Å. The extra-framework cations were added to the framework by simulated annealing using an energy minimizing canonical Monte Carlo algorithm,35 which is part of the Sorption software. The 96 cations needed to neutralize the framework were positioned leaving the composition Na96−xKxAl96Si96O394, 0 < x < 96. A total of 11 structure B
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
resource demanding, so instead effective charges were used. The effective charges were determined by averaging the DFT computed charges for each atom type in each zeolite NaKA structure, which were then plotted as a function of the K+/(K+ + Na+) ratio of respective structure to which a linear regression line was estimated by minimizing the composite sum of the squared deviations. These effective atomic charge values, for each K+/(K+ + Na+) ratio, were then appointed to respective structure to be used for the MD simulations. The averaged atomic charges are shown in Figure 2 together with their
compositions were created with cation ratios ranging from 0 to 1 K+/(K+ + Na+) ratio in steps of 0.1. We used the atomic ́ nchez et al.15 on zeolites charges presented by Garcia-Sá containing Na+ as a starting point for the simulated annealing. We also gave the Na+ and K+ equal charges, initially. This procedure gave very similar positions as the Reitveld refinement of PXRD data measured by Pluth and Smith4 and Liu et al.1 When the ions were positioned in the structure, the atomic charges for use in the GCMC and MD simulations, were determined through Hirshfeld populations analysis.36 The population analysis was carried out by performing a periodic density functional theory (DFT) energy calculation with the generalized gradient approximation (GGA) functional PBE37 using the DMOL3 software. For N2 and CO2 molecules, the bond lengths were optimized with DFT calculations using the B3LYP functional.38 For CO2, the atomic charges were determined with the electrostatic potential (ESP) fitting method,39 while the homonuclear diatomic N2 molecule was left neutral. To compare the equilibrium uptake of our model with that measured by Liu et al.,1 the adsorption of CO2 and N2 at 0.85 bar and 298 K was computed using GCMC. This comparison was done for each composition of zeolite NaKA through the K+/(K+ + Na+) range, with partial atomic charges as determined by DFT. Note that the GCMC method is purely statistical and does not take the gas molecules ability to move through the framework into account when estimating the spatial distribution of molecules in the framework. In each simulation cell, we excluded the forbidden volumes of the small beta cages for insertions, in order to restrict the gas molecules to be unphysically placed in these cages, which have pore openings too small to allow the passage of both CO2 and N2. Each GCMC simulation ran for 1 × 106 equilibrium steps and 1 × 107 production steps, and during the GCMC simulations, the framework was kept rigid. Keeping the framework rigid for the GCMC simulations can be done without affecting the results ́ nchez et al.40 substantially, which has been shown by Garcia-Sá Further, due to limitations of the Sorption software, the extraframework cations were not allowed to move during the GCMC simulations. In 2004, Calero et al.41 showed cases where nonmobile cations cause artifacts to the simulated adsorption isotherm. To estimate potential errors resulting from the fixed cation positions, we executed test cases on zeolite NaA and KA where respective number of CO2 and N2 molecules, as previously estimated by GCMC, were positioned simultaneously with the cations during the simulated annealing procedure. The gas molecules were then removed and the uptake could be measured on the structure with the new ion distribution corrected for the expected uptake. The deviations from the uptake of the noncorrected ion distributions were negligible (>0.15 mmol/g); hence, fixed cations are fully acceptable for the quality of this study. Additional computational details are presented in SI. To investigate how the self-diffusion of the adsorbates effect the uptake in the zeolite NaKA, MD simulations were performed. For these simulations, the simulation cell consisted of 2 × 2 × 2 unit cells as the use of the single unit cell yielded results of poor accuracy due to too few CO2 molecules. An even larger simulation cell would have been desirable but was currently not used due to limited computation time. These simulation cells were created in the same manner as those for the Monte Carlo simulations. Performing a Hirshfeld population analysis on this large simulation cell was too
Figure 2. Averaged Hirshfeld charges. The line is derived by minimizing the composite summed squared deviations for atomic type and fraction of K+. These effective atomic charges are applied to the system for the MD simulations.
respective regression line. The Hirshfeld charging procedure assigns the K+ ions an average charge of around 0.2e higher than the Na+ ions, which can reasonably be explained by the fact that the Na+ ions are smaller and more embedded in the framework, allowing a larger extent of charge transfer than for the larger K+ ions. We see from Figure 2 that the positively charged ions successively become effectively slightly less charged with an increasing K+ content in the framework, while the negative O-atoms strengthen their charges. This regression is distinctly linear, and we can assign the regression line with very small error. The Hirshfeld population analysis shows that atomic charges for the ions may vary in the three different ion positions, which cannot be taken into account as the ions diffuse and change positions. Furthermore, we would expect the partial atomic charges within the CO2 molecules to vary, due to local changes in electrical fields and field-gradients, when they are adsorbed on different sites within the framework. However, we are of the opinion that assigning the effective atomic charge distribution for each K+/(K+ + Na+) ratio is the most reasonable approximation given the classical MD software available today, where all charges specific to a certain atom-type are held constant throughout the length of the simulation. In the MD simulations, 52 gas molecules per unit cell (a total of 416 per simulation cell) were loaded and equilibrated by using the same energy minimization procedure that was used for the positioning of the cations. With an NVT ensemble at a temperature of 298 K and using the Nosé-Hoover thermostat, MD simulations were conducted with a simulation time of 2.0 ns and 1.0 fs time steps for each fully mobile (structure) composition. A mobile composition was needed to allow for zeolite deformations, which is necessary to simulate accurately the diffusion in contrary to adsorption simulations as was ́ shown by Garci a-Sá nchez et al.40 From the resulting trajectories, the self-diffusion coefficient (Ds) was determined by calculating the slope of the mean square displacement C
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
ion, and (3) with the Na+ ion removed. Also, the energies of the single Na+ and K+ ions were computed. The potential energy wells ΔE for each ion/ion-site combination was then computed according to
(displayed in Figure 3) in the linear region (t = 0.2−1.8 ns), according to the Einstein relation.
ΔE = E FW + E ion − E FW + ion
(3)
where EFW is the energy of the geometry optimized system 1, EFW+ion is the geometry optimized system 2 or 3, and Eion is the computed energy of the respective ion at the respective ion site.
■
RESULTS AND DISCUSSION The simulated equilibrium uptakes of CO2 and N2 for the pure zeolite NaA show a particularly strong agreement to those measured by Liu et al.1 (point of 0 atom % K+ content). This similarity indicates that our model should be reasonable. For zeolite NaA, we expect equilibrium adsorption simulations and experiment to agree as both CO2 and N2 translate rather unrestricted throughout the pore windows. When following the computed uptakes for an increasing K+/(K+ + Na+) ratio in the zeolite NaKA, as shown in Figure 4, we observe linearly Figure 3. Mean square displacement curves for CO2, Na+, and K+ estimated from molecular dynamics (MD) simulations.
Ds =
1 d lim 6 t →∞ dt
N
∑ ⟨[ri(t ) − ri(0)]2 ⟩
(1)
i=1
(Dαβ d )
Analogously, the “distinct diffusion coefficient” between the Na+ ions and CO2 molecules was computed, eq 2, in order to investigate the presence of any codependent diffusion between CO2 and Na+. Ddαβ =
1 d lim 6 t →∞ dt
Nα
Nβ
∑ ∑ ⟨[ra (t ) − ra (0)]·[rb (t ) − rb (0)]⟩ j
j
j
j
Figure 4. Uptake of CO2 and N2, experimental values, and Grand Canonical Monte Carlo simulated ones (298.15 K, 0.85 bar) in zeolite NaKA. Lines for respective data set are provided to guide the eye only: (a) CO2 sim, (b) CO2 exp, (c) N2 sim, and (d) N2 exp.
i=1 j=1
(2)
The “distinct diffusion coefficient” is a measure of the crosscorrelation in time between two distinct particles, a and b, of type α and β or, in our case, CO2 and Na+, in an equilibrated system.42 The “distinct diffusion coefficient” should not be confused with the “mutual diffusion coefficient”, which is commonly used to measure the codependence of multicomponent transport diffusion and correlates the nonequilibrium mass transfer of two species over time.43 The “distinct diffusion coefficient”, however, resulted in a negligible value, hence, no sign of codiffusion. Convergence of the diffusion coefficients were ensured by extending the simulation time for four selected configurations to 3.0 ns, two at low and two at high K+ content (0, 20, 80, and 100%). Prolongation of the simulations did not affect the resulting diffusion coefficients; hence, convergence was reached. To compare the potential energy well-depths of the Na+ and + K ions situated in the three different ion sites, DFT computations were performed on a cluster model of an αcage. The system was created by cutting out one α-cage unit from the energy-minimized zeolite NaA model used for the GCMC simulations. The cell was then saturated with H-atoms. One ion site was chosen, and all remaining atoms in the system were spatially fixed, while the ring atoms, as well as the occupying ion of the chosen ion site, were held flexible. Using the PBE functional and the DND basis set, geometry optimizations of the following three systems were performed: (1) the system as is, (2) with the Na+ ion exchanged for a K+
decreasing uptakes for both CO2 and N2. These linear dependencies are, however, not in agreement with the uptake drop shown experimentally around 10−20 atom % and 15−25 atom % K+ content for N2 and CO2, respectively. These deviations indicate that the dynamics involved in the gas separation is more complex than based on thermodynamic effects alone. Especially large is the relative deviation for the uptake of N2, which has a larger kinetic diameter than CO2 within microporous zeolites. The GCMC results show a CO2over-N2 selectivity, as estimated from single-component uptake between the gases, of 5−15, which agree with the selectivity measured experimentally in similar sorbents without size restricting pores, such as MFI-type9−12 and FAU-type13,14 zeolites, hence, allowing both gases to reach equilibrium. The next step in investigating the cation-dependent gas uptakes in zeolite NaKA is the analysis of the NVT-MD results. When analyzing the trajectories for CO2 (SI, Figures 1−3 and 7−9), Na+ (SI, Figures 4−6), and K+ (SI, Figures 10−12) presented in SI, it is apparent that the movements of the CO2 gas molecules and that of the cations are characteristically different. The CO2 molecules show a typical diffusive behavior in a 3-dimensional manner, cage-to-cage, penetrating the pore. The cations, on the other hand, move two-dimensionally on the surface of the framework, jumping from one site to the next closest site, whether located in the same pore window or the D
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
As the site-to-site ion-jumping is a nonlocal movement, seen from the pore windows, the rate of the ion-jumping can be measured with the same quantity as used for the CO2 diffusion, the self-diffusion coefficient. The local movements will only contribute to the MSD curve with an initial increase, while the nonlocal movements will sustain the slope of the MSD curve beyond diameters of the nonlocal movements. This enables a direct comparison between the ion jumping rates and the CO2 molecules possibility to diffuse between pores. The dependency of the self-diffusion coefficients on the K+/ (K+ + Na+) ratio are presented in Figure 6. What is remarkable
next neighboring pore window. When exiting an 8R window plane, the Na+ ion will first visit 6R or 4R sites and then eventually continue to a new 8R. The trajectories in the SI show that the 8R and 6R ion sites are energetically more favorable than the 4R site. The ions stay and vibrate around these sites for longer periods while constantly hopping between the energetically and geometrically equivalent sites. The 4R sites, on the other hand, merely act as stepping stones. This is also confirmed by the DFT calculations, which compare the potential energy wells of the different ion sites, as presented in Table 1. The 6R site has the deepest potential energy wells for Table 1. Potential Energy Well Depths for Na+ and K+ Ions at 4R, 6R, and 8R Sites of Zeolite NaA and KA Computed with Ab Initio Density Functional Theory (DFT)a
a
ion site
Na+
K+
4R 6R 8R
256 366 349
152 371 376
The energy values are presented in kcal/mol.
both cation types, 17 and 5 kcal/mol deeper than the 8R wells for Na+ and K+, respectively. Further, the 4R well depths drop substantially compared to the 6R wells, with decreases of 110 kcal/mol for Na+ and 219 kcal/mol for K+. These trends are also in accordance with the findings of Abrioux et al.26 that the mobility of Na+ ions in zeolite NaA increased in the order 6R, 8R, and then 4R, in respect to their initial ions site. Following this, the cation mobility can be divided into two different types: 1. Local movements: corresponding to jumps between equivalent sites within the same pore windows, as well as vibrations around the sites, as illustrated in Figure 5. 2. Nonlocal movements: corresponding to jumps between sites of different pore windows, which we refer to as siteto-site ion-jumping.
Figure 6. Self-diffusion coefficients of CO2, Na+, K+, and N+ + K+ combined in zeolite NaKA loaded with 52 CO2 molecules/UC, as estimated by NVT MD. Lines for respective data set are provided to guide the eye only.
with the dependencies are the differences in the mobility of the Na+ and K+ extra-framework cations and how the CO2 selfdiffusion, in turn, seemingly depends on them. The K+ ions are comparably static, independent of the K+/(K+ + Na+) ratio of the structure, where the K+ ions tend to sit firmly in their positions, restricted to move out of the 8R plane and, hence, blocking the passage of any CO2 molecule. For Na+, the selfdiffusion coefficient is notably higher for zeolite NaA as compared with NaKA. The difference is reasonably explained by the increasing number of cation sites being occupied by the relatively nonmobile K+, leaving fewer energetically favorable sites for the Na+ ions. This diffusional hindrance for Na+ is schematically shown for different cases in one-dimension for simplification by the cartoon showing three adjacent pore windows in Figure 7. In the cartoon, we assume that the Na+ ion only jumps to neighboring sites according to what is seen in the MD trajectories presented in SI. The results from the DFT cluster computations presented in Table 1 confirm that the K+ ion bonds stronger to the framework than the Na+ ion at the 8R site, with a potential energy well that is 27 kcal/mol deeper for K+ compared to that of Na+, which explains the immobility of the K+ ions. Also, at the 6R site, the K+ well is deeper; however, the difference is smaller, only 5 kcal/mol. The Na+ at the 4R site, on the other hand, bonds stronger to the framework than the K+ at the same site, with a difference of 103 kcal/mol. This can be explained by the larger size of K+ and that it is forced out approximately 0.8 Å farther from the plane as compared to Na+; hence, the interaction of K+ with the framework decreases at these 4R sites. Comparing the CO2 self-diffusion curve in Figure 6 with the weighted average of the apparent self-diffusion of Na+ and K+ cations, we observe that for pure zeolite NaA, CO2 has a selfdiffusion coefficient of more than double that of the cations. However, the shape of the CO2 self-diffusion coefficient curve
Figure 5. Three cation sites (I−III) located in the three different porewindow types (4R, 6R, and 8R) of zeolite A, together with the cross sections of the nonlocal ion movements within the pore windows, as approximated from the Molecular Dynamics trajectories presented in SI. There are four equivalent off-center Na+ ion sites and one centered K+ site in the plane of the 8R pore window (site II), while the 6R has two equivalent sites (site I) for both Na+ and K+, one on each side of the pore window plane as the ions are slightly displaced out from the plane. The 4R site (site III) also coordinates the ions out from the 4R plane and does not have an equivalent site, as the 4R is too narrow for the cations to jump through. E
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
trajectories. The number of visited windows drops successively for Na+, with an increasing K+ content according to the rootmean-square displacements (RMSD) displayed in Table 2, Table 2. Root Mean Square Displacement (RMSD) of CO2, K+, and Na+ at 1.8 ns of Simulation Time in the MD Simulationsa
a
K+/(K+ + Na+) (%)
CO2
Na+
0 10 20 30 40 50 60 70 80 90 100
15.6 14.6 14.0 11.3 10.5 9.5 7.0 6.7 6.6 6.0 6.2
9.8 8.4 7.5 5.4 5.8 4.7 3.7 4.0 3.9 3.8
K+ 2.4 2.2 2.3 2.1 2.3 2.1 2.1 2.0 2.1 2.1
The distance values are presented in Å.
while besides jumping between equivalent sites within the initial plane, the K+ ions in zeolite KA will only visit on an average 1.4 different sites during the simulation time. The RMSD for K+ ranges between 2 and 2.4 Å throughout the range of K+ coverage. From the K+ trajectories in SI, Figures 10−12, the vibrational diameters around the 4R and 8R sites and the distance between equivalent 6R sites were estimated to 0.7, 1.6, and 2.0 Å, respectively, as presented in Figure 5. Hence, the majority of the K+ movements are localized to one pore window, while the remainder of the RMSD is due to the relatively few site-to-site jumps, resulting in a small but constant slope of the MSD curve. To confirm that this slope is contained sufficiently far beyond RMSD distances, which could be attributed to local movements, the MD simulation for CO2 in zeolite KA was extended to 15.0 ns. The slope of the MSD curve remained constant, and the RMSD reached 5.0 Å, which is well beyond the distances of local movements. Hence, the rate of nonlocal K+ jumps was sustained. These differences in mobility patterns between the cations and gas molecules explains why the codiffusion coefficient is negligible. Hence, the migration of the Na+ ions facilitate the pore-to-pore movement of CO2 without involving any codiffusion between the species. We also computed the diffusion coefficients of Na+ and K+ in the pure zeolite NaA and KA, respectively, without the presence of CO2, for a comparison to see if the CO2 molecules had an effect on the mobilities of the ions. From the simulations, we observed that the diffusion of Na+ and K+ increased by 45 and 35%, respectively, in the presence of CO2. This increase is in line with the recent work by Shang et al.44 who showed with plane-wave DFT computations that the transition state of a K+ ion departing from the 8R site is lowered by the presence of CO2 in the pores. Corresponding MD simulations are also performed for N2; however, the simulations are not able to detect any macroscopic diffusion for the gas in the adsorbent even for low K+ content where we would expect to see certain differentiation following the same reasoning as for CO2. This inability is likely to be due to the model used in this work. It appears to not be sufficiently sensitive to measure the very slow diffusion of N2 in this material. A larger simulation cell but, most importantly, longer simulations would improve the statistics to obtain
Figure 7. For cases a, b, the CO2 molecule is free to pass through the 8R pore window as the Na+ ion situated in the center is free to jump to the ion sites both sides (a) or one side (b) as the Na+ ions situating these sites, in turn, are mobile in contrast to the K+ ion. For cases c, d, the CO2 molecule will not pass through the pore window. In case c, the Na+ ion in the center has no sites to jump to and is locked in its position, while in case d, the rigid K+ ion in the center pore-opening will not budge for the CO2 molecule.
resembles that of the cations, both starting with a drop and then to an extent planing out at a similar low value as for the pure zeolite KA. This shape suggests that the mobility of CO2 is restricted by that of the extra-framework ions. The mobile Na+ ions located in the 8R window allow the gas to diffuse relatively unrestricted through the zeolite, while the K+, in contrast, effectively blocks out the 8R window and, in turn, hinders the diffusion of the gas. When analyzing the self-diffusion of CO2, the curve can be divided into two linear regions. The first region between 0 and 60% K+ content has a steep descent which, following our reasoning, mainly results from the decreasing mobility of the pore-blocking cations. In the second linear region, between 60 and 100% K+ content, the slope decreases, as there is very little pore-to-pore diffusion in this region. The remaining negative slope can mainly be explained by the decreasing pore volumes following the increasing content of the larger K+ ions leaving less space for movement of the gas molecules within the pores. On average, the Na+ ions in the pure zeolite NaA visit 5.1 different pore windows (1.3 8R, 2.0 6R, and 1.8 4R) during the simulation time determined from 300 randomly selected F
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
■
CONCLUSIONS We have performed Grand Canonical Monte Carlo and MD simulations for CO2 in zeolite NaKA to estimate how the molecular sieving functionality of the material effects the total uptake for different K+/(K+ + Na+) ratios. It was necessary to combine the results from both thermodynamic (GCMC) and dynamic (MD) simulations in order to rationalize the uptake behavior of CO2, which contributes to the high CO2-over-N2 selectivity measured on zeolite NaKA. This uptake behavior can be ascribed to the difference in mobility of the ions. This difference results in the nonlinear uptake drop seen experimentally for both CO2 and N2 at critical K+-contents, that for N2 coming before that for CO2, leaving an opening of exceptionally high selectivity while the CO2 uptake still remains high. Hence, the sieving ability for CO2 in this zeolite NaKA is a consequence of the differences in mobility of the different cation types and not simply a size effect. We expect similar effects also for other 8-ring zeolites with low Si-to-Al ratios. As we were unsuccessful in simulating the diffusion of N2, no conclusions regarding the total uptake of N2 or the resulting CO2-over-N2 selectivity could be drawn. However, we would expect to see the same effect on N2, which would explain the high selectivity measured experimentally. Macroscopic experiments are sensitive to effects on multiple length and time scales and cannot be always directly correlated to the microscopic level. For this reason, it is essential to ask specific questions regarding certain dynamic patterns, which may be simulated on an atomistic level and, in turn, predict or explain a large scale pattern. Despite that the time scales involved in such experiments are many magnitudes longer compared to those which are possible to achieve with atomistic simulations, the distinct differentiation between the diffusion of the two ion types, which we have predicted in our simulations, should hold true even at larger time scales. However, we cannot exclude additional effects which may influence the adsorption dynamics on a larger scale. Work is in progress by two of us, taking a more quantitative approach to describe this mechanism. We are performing Ab Initio Molecular Dynamics (AIMD) simulations on this system for the purpose of accurately computing the free energy barriers for diffusion for the different sorbates and cations,46 allowing us to determine accurately the transition state hopping rates. This information may be beneficially used as input for Transition State Theory based coarse-graining methods such as the kinetic Monte Carlo, enabling accurate diffusion simulations on up to a microsecond time scale, even for N2. Important work on coarse-graining methods treating gas diffusion in zeolites have been done by Tunca and Ford47 among others.48−50 In addition, we are planning experimental work to investigate the phenomenon presented in this paper. 23Na NMR spectroscopy can be used to estimate the mobility of the Na+ ions as a function of the K+ content in the zeolite A composition. If we can show that the Na+ ions in the zeolite have a decreasing diffusion constant with increasing K+ content, this should be sufficient to support our thesis that the K+ ions hinder the diffusion of the Na+ ions and, in turn, the CO2 diffusion and total uptake in the sorbent.
sufficient data to sense such slow diffusion. Combining the identified thermodynamic and kinetic contributions to the CO2 uptake dependence on the K+ content, measured by two different simulation methods, is likely complex. Nonetheless, we will use a very simple, straightforward, and qualitative approach in an attempt to combine and compare the results from the simulations with the experimental results. The GCMC results provide a measure of the equilibrium uptake. Given infinite time in experimental measurements these values should agree, providing that the potential landscape is accurately described. However, in experiments as well as in any potential applications the adsorption cycles are defined in time and are relatively short and constant. The blocking of the 8R pore windows hinders the free diffusion of the gas molecules within the zeolite framework. This will limit the possibility of the system reaching its equilibrium loading within this time. It is this aspect that potentially allows kinetic effects to be employed to lower the parasitic energy contributions to gas separation as compared to fully equilibrium based processes. Similar reasoning was used by Reyes and Ruthven when analyzing adsorption driven gas separation that is kinetically enhanced.45 Following this argument, we can assume that there is a direct correlation between the diffusion rate of the adsorbate and the rate at which the system will approach the equilibrium loading level. As a correction factor to the equilibrium uptake Ueq, we define the root relative diffusion constant (√Drel s ), according to eq 4, as a measure of the extent of freedom that the gas molecules have to diffuse and reach equilibrium. rel Ds
⎛ D ⎞1/2 = ⎜ s0 ⎟ ⎝ Ds ⎠
(4)
Hence, the relative diffusion for the pure zeolite NaA is set to 1, where D0s refers to the self-diffusion coefficient value of zeolite NaA, as we assume that the diffusion is unrestricted by poreblocking at this composition. The total uptake Utot can then be estimated by combining the results from the GCMC and NVTMD simulations by scaling the equilibrium uptakes with this correction factor according to Utot =
rel
D s Ueq
Article
(5)
We see a substantial improvement to the agreement with the experimental CO2 uptake curve in Figure 8 when scaling the computed adsorption of CO2 in this manner.
■
Figure 8. Experimental (a) gas uptakes and combined NVT-MD and GCMC simulated uptakes according to eq 5 (b) for CO2 (0.85 bar, 298 K). Lines for respective data set are provided to guide the eye only.
ASSOCIATED CONTENT
* Supporting Information S
The first section, “Computational Details”, provides computational details for the MC and MD simulations. The second G
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
section, “Force Field Details” provides the reader with the COMPASS force field functional forms and certain parameters used in the MD and MC simulations. The third section, “MD Trajectories”, provides a graphical view of the CO2 (Figures 1− 3) and Na+ (Figures 4−6) trajectories in zeolite NaA, as well as CO2 (Figures 7−9) and K+ (Figures 10−12) in zeolite KA. Three figures for each species is provided, each figure showing eight randomly selected trajectories of the selected species. This material is available free of charge via the Internet at http:// pubs.acs.org.
■
(15) García-Sánchez, A.; Ania, C. O.; Parra, J. B.; Dubbeldam, D.; Vlugt, T. J. H.; Krishna, R.; Calero, S. Transferable Force Field for Carbon Dioxide Adsorption in Zeolites. J. Phys. Chem. C 2009, 113, 8814−8820. (16) Pillai, R. S.; Peter, S. A.; Jasra, R. V. Correlation of Sorption Behavior of Nitrogen, Oxygen, and Argon with Ca2+ Locations in Zeolite A: A Grand Canonical Monte Carlo Simulation Study. Langmuir 2007, 23, 8899−8909. (17) Maurin, G.; Belmabkhout, Y.; Pirngruber, G.; Gaberova, L.; Llewellyn, P. CO2 Adsorption in LiY and NaY at High Temperature: Molecular Simulations Compared to Experiments. Adsorption 2007, 13, 453−460. (18) Maurin, G.; Llewellyn, P. L.; Bell, R. G. Adsorption Mechanism of Carbon Dioxide in Faujasites: Grand Canonical Monte Carlo Simulations and Microcalorimetry Measurements. J. Phys. Chem. B 2005, 109, 16084−16091. (19) Jaramillo, E.; Chandross, M. Adsorption of Small Molecules in LTA Zeolites. 1. NH3, CO2, and H2O in Zeolite 4 Å. J. Phys. Chem. B 2004, 20155−20159. (20) Akten, E. D.; Siriwardane, R.; Sholl, D. S. Monte Carlo Simulation of Single- and Binary-Component Adsorption of CO2, N2, and H2 in Zeolite Na-4A. Energy Fuels 2003, 17, 977−983. (21) Krishna, R.; van Baten, J. M. A Molecular Dynamics Investigation of a Variety of Influences of Temperature on Diffusion in Zeolites. Microporous Mesoporous Mater. 2009, 125, 126−134. (22) 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. (23) Madison, L.; Heitzer, H.; Russell, C.; Kohen, D. Atomistic Simulations of CO2 and N2 within Cage-Type Silica Zeolites. Langmuir 2011, 27, 1954−1963. (24) Ramsahye, N. A.; Bell, R. G. Cation Mobility and the Sorption of Chloroform in Zeolite NaY: Molecular Dynamics Study. J. Phys. Chem. B 2005, 109, 4738−4747. (25) Maurin, G.; Plant, D. F.; Henn, F.; Bell, R. G. Cation Migration upon Adsorption of Methanol in NaY and NaX Faujasite Systems: A Molecular Dynamics Approach. J. Phys. Chem. B 2006, 110, 18447−54. (26) Abrioux, C.; Coasne, B.; Maurin, G.; Henn, F.; Jeffroy, M.; Boutin, A. Cation Behavior in Faujasite Zeolites upon Water Adsorption: A Combination of Monte Carlo and Molecular Dynamics Simulations. J. Phys. Chem. C 2009, 113, 10696−10705. (27) www.accelrys.com/products/materials-studio. (28) http://accelrys.com/products/datasheets/compass.pdf. (29) Waldman, M.; Harver, A. T. New Combining Rules for Rare Gas van der Waals Parameters. J. Comput. Chem. 1993, 14, 1077− 1084. (30) Sun, H. COMPASS: An Ab Initio Forcefield Optimized for Condensed-Phase Applications: Overview with Details on Alkane and Benzene Compounds. J. Phys. Chem. B 1998, 102, 7338−7364. (31) Yang, J.; Ren, Y.; Tian, A. COMPASS Force Field for 14 Inorganic Molecules, He, Ne, Ar, Kr, Xe, H2, O2, N2, NO, CO, CO2, NO2, CS2, and SO2, in Liquid Phases. J. Phys. Chem. B 2000, 104, 4951−4957. (32) Zhao, L.; Liu, L.; Sun, H. Semi-Ionic Model for Metal Oxides and Their Interfaces with Organic Molecules. J. Phys. Chem. C 2007, 111, 10610−10617. (33) Peng, Z.; Ewig, C. S.; Hwang, M.-J.; Waldman, M.; Hagler, A. T. Derivation of Class II Force Fields, 4. van der Waals Parameters of Alkali Metal Cations and Halide Anions. J. Phys. Chem. A 1997, 101, 7243−7252. (34) Gramlich, V.; Meier, W. M. The Crystal Structure of Hydrated NaA: A Detailed Refinement of a Pseudosymmetric Zeolite Structure. Z. Kristallogr. 1971, 133−149, 139−149. (35) Kirkpatrick, S.; Gelatt, C. D.; Vecchi, M. P. Optimization by Simulated Annealing. Science 1983, 220, 671−680. (36) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge Densities. Theor. Chim. Acta 1977, 44, 129−138.
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work is supported by the Swedish Science Council (VR) and the EXSELENT Center for porous materials at Stockholm University. A.M. and A.L. wish to thank Swedish National Infrastructure for Computing (SNIC) for generous allocation of computing time and Andrea Gabrieli, Dr. Marco Sant, and Prof. Pierfranco Demontis for fruitful discussions. N.H. thanks the Swedish Energy Agency for a grant.
■
REFERENCES
(1) Liu, Q.; Mace, A.; Bacsik, Z.; Sun, J.; Laaksonen, A.; Hedin, N. NaKA Sorbents with High CO2-over-N2 Selectivity and High Capacity to Adsorb CO2. Chem. Commun. 2010, 46, 4502−4504. (2) Aaron, D.; Tsouris, C. Separation of CO2 from Flue Gas: A Review. Sep. Sci. Technol. 2005, 40, 321−348. (3) Lin, L.-C.; Berger, A. H.; Martin, R. L.; Kim, J.; Swisher, J.; Jariwala, K.; Rycroft, C. H.; Bhown, A. S.; Deem, M. W.; Haranczyk, M.; et al. In Silico Screening of Carbon-Capture Materials. Nat. Mater. 2012, 11, 1−9. (4) Pluth, J. J.; Smith, J. V. Accurate Redetermination of Crystal Structure of Dehydrated Zeolite A. Absence of Near Zero Coordination of Sodium. Refinement of Silicon, Aluminum-Ordered Superstructure. J. Am. Chem. Soc. 1980, 102, 4704−4708. (5) Crank, J. Zeolite Molecular Sieves; Wiley & Sons: New York, 1974. (6) Breck, D. W.; Eversole, W. G.; Milton, R. M.; Reed, T. B.; Thomas, T. L. Crystalline Zeolites. I. The Properties of a New Synthetic Zeolite, Type A. J. Am. Chem. Soc. 1956, 78, 5963−5971. (7) Yeh, Y. T.; Yang, R. T. Diffusion in Zeolites Containing Mixed Cations. AIChE J. 1989, 35, 1659−1666. (8) Myers, A.; Prausnitz, J. Thermodynamics of Mixed-Gas Adsorption. AIChE J. 1965, 11, 121−127. (9) Bernal, M. P.; Coronas, M.; Menéndes, M.; Santamaría, J. Separation of CO2/N2 Mixtures Using MFI-Type Zeolite Membranes. AIChE J. 2004, 50, 127−135. (10) Sirwardane, R. V.; Shen, M. S.; Fisher, E. Adsorption of CO2 on Zeolites at Moderate Temperatures. Energy Fuels 2005, 19, 1153− 1159. (11) Gardner, T. Q.; Flores, A. I.; Noble, R. D.; Falconer, J. L. Transient Measurements of Adsorption and Diffusion in H-ZSM-5 Membranes. AIChE J. 2002, 48, 1155−1167. (12) Harlick, P. J. E.; Tezel, F. H. Adsorption of Carbon Dioxide, Methane and Nitrogen: Pure and Binary Mixture Adsorption for ZSM5 with SiO2Al2O3 Ratio of 280. Sep. Pur. Technol. 2003, 33, 199−210. (13) Hasegaw, T.; Watenabe, K.; Kusakabe, K.; Morooka, S. The Separation of CO2 Using Y-Type Zeolite Membranes Ion-Exchanged with Alkali Metal Cations. Sep. Pur. Technol. 2001, 22−23, 319−325. (14) Kusakabe, K.; Kuroda, T.; Uchino, K.; Hasegawa, Y.; Marooko, S. Gas Permeation Properties of Ion-exchanged Faujasite-type Zeolite Membranes. AIChE J. 2004, 45, 1220−1226. H
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
(37) Becke, A. D. Density Functional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (39) Singh, C. U.; Kollman, P. A. An Approach to Computing Electrostatic Charges for Molecules. J. Comput. Chem. 1984, 5, 129− 145. (40) García-Sánchez, A.; Dubbeldam, D.; Calero, S. Modeling Adsorption and Self-Diffusion of Methane in LTA Zeolites: The Influence of Framework Flexibility. J. Phys. Chem. C 2010, 114, 15068−15074. (41) 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, 11376− 11385. (42) Raineri, F. O.; Friedman, H. L. Velocity Correlations in the Molecular Dynamics Ensemble: Computation of the Distinct Diffusion Coefficients. J. Chem. Phys. 1989, 91, 5642−5647. (43) Breck, D. W. The Mathematics of Diffusion; Oxford University Press: New York, 1975. (44) Shang, J.; Li, G.; Singh, R.; Gu, Q.; Nairn, K.; Bastow, T.; Medhekar, N.; Doherty, C. M.; Hill, A. J.; Liu, J. Z.; et al. Discriminative Separation of Gases by a “Molecular Trapdoor” Mechanism in Chabazite Zeolites. J. Am. Chem. Soc. 2012, 136, 19246−19253. (45) Ruthven, D.; C., R. S. Adsorptive Separations of Light Olefins from Parafins. Microporous Mesoporous Mater. 2007, 104, 59−66. (46) Mace, A.; Laasonen, K.; Laaksonen, A. Free Energy Barriers for CO2 and N2 in Zeolite NaKA: An Ab Initio Molecular Dynamics Approach. Phys. Chem. Chem. Phys. 2013, in press, DOI: 10.1039/ c3cp52821a. (47) Tunca, C.; M., F. D. A Transistion-State Theory Approach to Adsorbate Dynamics at Arbitrary Loadings. J. Chem. Phys. 1999, 111, 2751−2760. (48) Gabrieli, A.; Demontis, P.; Pazzona, F. G.; Suffritti, G. B. Speeding Up Simulation of Diffusion in Zeolites by a Parallel Synchronous Kinetic Monte Carlo Algorithm. Phys. Rev. E 2011, 83, 056705−1−8. (49) Abouelnasra, M. K. F.; Smit, B. Diffusion in Confinement: Kinetic Simulations of Self- and Collective-Diffusion Behavior of Adsorbed Gases. Phys. Chem. Chem. Phys. 2013, 14, 11600−11609. (50) Saravanan, C.; Auerbach, S. Theory and Simulation of Cohesive Diffusion in Nanopores: Transport in Subcritical and Supercritical Regimes. J. Chem. Phys. 1999, 110, 11000−11011.
I
dx.doi.org/10.1021/jp4048133 | J. Phys. Chem. C XXXX, XXX, XXX−XXX