PCCP View Article Online
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 166
View Journal | View Issue
Free energy barriers for CO2 and N2 in zeolite NaKA: an ab initio molecular dynamics approach Amber Mace,a Kari Laasonenb and Aatto Laaksonen*a
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Ab initio Molecular Dynamics (AIMD) is used with spatial constraints to estimate the free energy barriers of diffusion for CO2 and N2 gas molecules in zeolite NaA and KA. We investigate the extent to which the diffusion of these gas molecules is hindered, in the two separate cases of a smaller Na+ ion or a larger K+ ion blocking the 8-ring pore window. In contrast to classical Molecular Dynamics, AIMD performs these computations accurately and unbiased in the absence of empirical parameterization. Our work has resulted in stable and reliable force profiles. The profiles show that the larger K+ ion effectively blocks Received 5th July 2013, Accepted 25th October 2013
the passage of both CO2 and N2 molecules while the smaller Na+ ion will allow both molecules to pass.
DOI: 10.1039/c3cp52821a
effect, which the size of the respective cation occupying the pore window has on diffusive properties of
www.rsc.org/pccp
each gas molecule. Hence, this effect can be altered through ion exchange to fine-tune the functionality of a specific zeolite as a molecular sieve.
These results are a quantitative demonstration of the concept of pore blocking where we compute the
Introduction Certain zeolites with narrow window apertures, such as those of zeolite A (framework type code LTA), coincide with the approximate size of small gaseous molecules such as CO2 and N2. These may be interesting candidates for adsorbents with swing adsorption technologies due to their potential molecular sieving capabilities and otherwise attractive properties. Zeolite A is a threedimensional cage type microporous aluminosilicate with an Al-toSi ratio of 1 resulting in a high cation content, with respect to the zeolite framework, which acts to neutralize the framework. This gives the framework neutralized by Na+ and K+ ions (zeolite NaKA), the chemical composition Na96xKxAl96Si96O394, 0 o x o 96. With the point group FM3c it has a cubic periodicity built up by so-called a- and b-cages where one unit cell consists of eight (2 2 2) a-cages as shown in Fig. 1. The a-cages intersect by 8-ring (8R) pore windows, where the number of oxygen atoms within the ring, confining the pore window, determines the window ring number. The framework also contains 6-ring (6R) and 4-ring (4R) pore windows, however it is only the 8R pore windows, which are large enough to allow the passage CO2 or N2. As shown by PXRD measurements,1,2 zeolite A has three different cation sites; site I coordinated with the centre of the 6R, site II in the plane of the 8R and site III coordinated with the centre of the 4R. When zeolite A solely a
Arrhenius Laboratory, Stockholm University, 10691 Stockholm, Sweden. E-mail: [email protected]; Fax: +46-8-152187; Tel: +46-8-162372 b Department of Chemistry, Aalto University, P.O. Box 16100, 00076 Espoo, Finland. E-mail: [email protected]; Tel: +358-40-5570044
166 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Fig. 1 An LTA unit cell (left) is built up by eight a-cages (bottom right). The smaller b-cages (top right) are comprised in the cubic intersect of the a-cages.
contains monovalent extra-framework cations, each of the 64 6R and 24 8R per unit cell will generally situate one cation, while the remaining eight cations are distributed among the 4R pore windows. Furthermore, the Na+ and K+ ions have slightly different positions relative to the pore window rings as illustrated in Fig. 2. Site II for K+ is centred in the 8R, while the smaller Na+ is approximately 1.3 Å off-centre.2 This off-centre placement results in four energetically and geometrically equivalent sites for Na+ in the 8R plane. Zeolite A neutralized with Na+ cations is referred to as zeolite NaA or 4A and with K+, zeolite KA or 3A. Here 4A and 3A denote the resulting pore window diameters of approximately 4 Å and
This journal is © the Owner Societies 2014
View Article Online
Paper
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Fig. 2
PCCP
Crystallographic K+ and Na+ positions in zeolite A.
3 Å respectively. This can be compared to the kinetic diameters of CO2 and N2, which are 3.3 Å and 3.6 Å, respectively.3 The kinetic diameter refers to the smallest effective dimension of the molecules, the critical molecular dimension. This dimension is determined from the separation distance, which corresponds to the minimum equilibrium separation of the Lennard-Jones potential. Hence, the stronger attractive forces of CO2 can rationalize this somewhat counter-intuitive relation that N2 is, in this aspect, larger than CO2. Following this, both CO2 and N2 should be able to pass through the 8R pore windows accommodating Na+, while a K+ ion in the same position will block the passage. In the recently published paper by Liu et al.,1 particularly appealing results were presented showing that through Na+-toK+ ion exchange in zeolite 4A an exceptionally high CO2-over-N2 selectivity could be reached while the total CO2 uptake still remained high at a specific K+/(K+ + Na+) ratio of 17 at%. Liu et al.1 predicted that this high selectivity was due to a sieving effect as a result of an effective pore diameter lying right in between the kinetic diameters of CO2 and N2 allowing the material to adsorb CO2 while effectively blocking out N2. In a separate work reported by two of the current authors,4 classical simulation methods have been utilized in order to model the separate contributions of thermodynamics and kinetics involved in the gas uptake in zeolite NaKA, and how these contributions are affected by the two different ion species. By combining results from classical Molecular Dynamics (MD) and Grand Canonical Monte Carlo (GCMC) simulations the experimental results for CO2 were qualitatively reproduced, indicating that the dynamics behind the high selectivity measured by Liu et al.1 was, in fact, a combination of these contributions. In this work we attempt to go further and look in detail at the dynamics involved, when the gas molecules diffuse from pore-to-pore, by using ab initio Molecular Dynamics (AIMD). For zeolites and other microporous materials, AIMD methods have been mainly used when investigating catalysis and chemisorption of guest molecules in the frameworks.5–10 While, typically, classical simulation methods such as GCMC and classical MD are used to model the physisorbed equilibrium uptakes11–15 and diffusion coefficients16–20 of gas molecules in zeolitic systems. These classical methods use predetermined force fields to model each inter-atomic interaction and are
This journal is © the Owner Societies 2014
commonly parameterized to empirical data sets, such as adsorption isotherms. A common procedure when fitting empirical force field parameters for zeolite–gas interactions is to use GCMC simulations to match an experimentally measured adsorption isotherm. This approach may be a good estimation when investigating gas uptake with zeolites without size restricting pores. However, when an unrestricted pore-to-pore diffusion of the adsorbent is, to some extent, hindered by the size or cation blocking of the pore windows, such as in the case of zeolite A, the uptake will be affected. In this case, determining a force field by appointing the full uptake to thermodynamics will create a bias. Also combinations of different force field parameters may give the same outcome, which is apparent when comparing the fairly large amount of zeolite NaA–CO2 force fields with widely varying parameters for the same interactions.11,13–15 This may not have a great effect on a larger scale, but when modelling a specific interaction, the outcome may vary greatly between different force fields; hence the results are not accurate. Further, recent work by Garcı´a´nchez et al.21 showed that the choice of force field does, in Sa fact, strongly affect the outcome of diffusion simulations. With this approach based on periodic density functional theory (DFT) and empirical van der Waals corrections an accurate and unbiased description of the molecular interactions can be obtained. We can investigate how the forces act on the different gas molecules when moving through the pores as well as what effect the cation type has on this. Also the energy barriers for the passage through the pores can be computed. The DFT does not need specific parameterization for this system as is the case for classical molecular dynamics with empirical force fields and thus the DFT can really predict the results. In addition, it is easy to investigate the dynamics of the N2 molecules compared with CO2, which is problematic with classical MD due to the very slow diffusion of N2 molecules in the framework. The down side of this approach is the high computational cost forcing us to narrow down the computations to a smaller simulation cell and shorter simulation length compared to what is possible for classical MD. Also, it is necessary to use constraints in order to look at the specific dynamics of interest. From this point of view, this work is well needed and can be seen as pioneering in the aspect of using AIMD to investigate the dynamics of physisorbed guest molecules in microporous materials.
Computational details For the AIMD simulations, the CP2K software package is utilized. CP2K performs Born Oppenheimer Dynamics, which calculates DFT forces ‘‘on the fly’’ with the QUICKSTEP22 module using the GPW method.23 Here, the atom centred Gaussian basis sets are used together with augmented plane wave basis sets. For the DFT force calculations we used the PBE functional and the DZVP-MOLOPT-SR basis set together with GTH pseudopotential for all atom types.24 For all atoms, only the valence electrons were treated explicitly with the DFT. Both
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 167
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
PCCP
Na+ and K+ have only one explicit valence electron. The plane wave kinetic energy cutoff was set to 280Ry. Also, the empirical vdW correction, DFT-D3, by Grimme25 was used. The simulations were run with a 1.5 fs time-step where the equations of motion were propagated using the Velocity-Verlet algorithm. Further, the simulations were run in an NVT ´–Hoover thermostat with the ensemble employing the Nose temperature set to 298.15 K. The free energy was calculated by using the constraint method.26 This method allows a direct estimation of the mean force from the time average of the force on the geometric constraint by applying the blue-moon ensemble method.27 The free energy is then estimated by integrating the mean force curve. The used constraint will be explained later. The zeolite NaA structure was created in the Accelrys Materials Studio Suite.28 The coordinates were taken from the structure database and were originally presented by Gramlich and Meier.29 The Sorption module was used to position the 96 extra-framework Na+ ions per unit cell in the framework using a Canonical Monte Carlo ensemble simulated annealing energy minimizing procedure.30,31 As a simulation cell in the AIMD computations 1/4 of a unit cell was used, corresponding to two intersecting a-cages on the x-axis, as shown in Fig. 3, and periodic boundary conditions were implemented with cell parameters a, b, c = 24.6 Å, 12.3 Å, 12.3 Å. One gas molecule is placed in the centre of the left a-cage. As we attempt to measure and compare the energy barriers for the gas molecules to pass through the pore window from one a-cage to another we use a constraint to steer the gas molecule through the intersecting pore window. This is implemented by using a point-to-plane constraint where the plane is defined by three oxygen atoms in the 8R and the point is one atom in the gas molecule, for CO2 the oxygen atom and for N2 one of the nitrogen atoms. Hence, the entropy contribution comes from the molecular movement around the fixed atom. The point-to-plane distance is visualized in Fig. 3. With this constraint the gas molecule is free to move within a plane at the defined distance from the window plane in the AIMD simulation run as well as rotate around the constrained point-atom.
Fig. 3 The simulation cell consists of two of the zeolite A building units, a-cage. One gas molecule (here CO2) is placed in the left a-cage and a constraint is set confining the distance between the oxygen of the CO2 molecule and the left 8R window plane.
168 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Paper
Simulation runs are set up with a subsequently decreasing point-to-plane distance. Near the force maxima a typical distance step was 0.25 Å and around the centre of the cage the distance step was 1.0 Å. For each constraint value an AIMD simulation of length of at least 10 ps is run to ensure good convergence of the constraint force, since the force has large fluctuations. To check the stability of the system the same procedure was performed in the opposite direction, where the gas molecule was steered towards the 8R in the centre of the simulation cell. The next step was then to look at how the ion exchange affects the system. We exchanged the Na+ ion in the left 8R in the simulation cell for a K+ ion and proceeded to perform the simulations in the same manner as for the previous system for both the CO2 and N2 molecules. Finally, for a comparison the same simulations were performed in the system where the ion was removed completely from the 8R pore window.
Results and discussion The AIMD computed forces for each distance step simulation run are presented as a function of the constrained point-toplane distance in Fig. 4 and 5. In Fig. 4 the force data are presented for the movement of CO2 and N2 molecules through the whole length of the Na+ containing zeolite cage, from the left 8R to the centre 8R. This picture displays the clear force differences for the 8R-approaching CO2 compared with N2 where the latter experiences a larger repulsive force. This is in line with the hypothesis that it is more difficult for the larger N2 molecules to pass through an 8R, which situates an ion. The energy barriers are calculated by integrating the forces in the point-to-plane distance range 0–4 Å where the approximated barrier for the N2 diffusion is 27 kJ mol1 and 21 kJ mol1 for CO2, giving a barrier difference of 6 kJ mol1. From the Arrhenius equation, this will correspond to a transition state relation of around 1 : 11. From the trajectories of CO2 and N2 at point-to-plane distances o1 Å, we observe
Fig. 4 The absolute values of the averaged forces are, for CO2 and N2 in the pure Na+ zeolite A, plotted as a function of the distance between the 8R window plane and the COM of the gas molecule.
This journal is © the Owner Societies 2014
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Paper
Fig. 5 The absolute values of the averaged forces for CO2 and N2 in the pure Na+ zeolite A as well as for the K+ ion exchange structure are plotted as a function of the constrained distance between the 8R window plane and Oxygen of CO2 and one N-atom of N2. The energy barrier values are presented in the parenthesis and are defined by the integral of the force in the range 0–4.0 Å point-to-plane distances.
differences in the behaviour of the gas molecules that can help to rationalize the difference in the energy barriers. N2 tends to point towards the cation, resulting from an induced dipole interaction with the Na+ ion. This effect of strong steering counteracts the N2 molecule from penetrating the pore. For the CO2, on the other hand, the molecule interacts both with the Na+ and the cage. The CO2 is thus oriented more perpendicular with the pore plane due to the balancing of CO2 interactions between the Na+ ion and the cage. Overall, the CO2 is affected less than N2 by the Na+ ion, and this will result in a decreased transition state for CO2 to pass through the pore window. The other notable difference between the CO2 and N2 curves is the behaviour in their centre regions. The CO2 potential is not zero at the centre of the zeolite cage, which is due to the stronger interaction of the CO2 molecule with the zeolite cage. This stronger interaction is due to the larger quadrupole moment of CO2 and the bending vibrations of CO2, causing the molecule to have an instantaneous dipole moment. Even when the CO2 molecule is in the centre of the pore the environment is not exactly symmetric, as the Na+ ion is situated off the centre of the 8R and the centred CO2 molecule generally attracts an 8R ion in the xy or the xz plane. Hence, the molecule will be drawn more to the left or the right depending on which side of the 8R the Na+ is located. This causes the non-zero force even for the CO2 around the 6 Å position. The force integral in the centre region of the zeolite cage (3–9.6 Å) is 19 kJ mol1 for CO2 compared to 6.0 kJ mol1 for N2. This gives us a crude estimation of the binding energy difference between CO2 and N2. The real force curve is antisymmetric, so the numbers above need to be divided by 2 and the energy difference will be ca. 6.5 kJ mol1. At room temperature this will lead to an adsorption selectivity of 1 : 15 favouring CO2. This is a good agreement with the GCMC simulations of Mace
This journal is © the Owner Societies 2014
PCCP
et al.,4 where the selectivity was in the range of 1 : 10–1 : 20, depending on the K+ content. Another important observation from Fig. 4 is that the force curves are symmetric and the data scattering is relatively small. This indicates that the simulations are long enough to yield well-converged forces with error bars varying between 0.2–1 103 a.u. These force fluctuations are especially large around 2–3 Å from the planes (0.5–1 103 a.u.), hence the points at distances 2, 2.5, 3, 9, 10 and 10.5 Å. This can be explained by the gas molecule being at a critical distance from the pore where the rotations have a large effect on the force and the difference between the gas molecules being either parallel or perpendicular with the plane is quite large. As the distance to the plane decreases, the gas molecule stops rotating since the space becomes more constricted, and in the point-to-plane distance region 0–2 Å the force fluctuations are comparably small (0.2–0.3 103 a.u.). If the distance to the plane increases, the molecules will rotate, but will not cause large fluctuations of the force (0.3–0.5 103 a.u.) as the interaction with the plane is smaller. The next set of simulations focuses on the role of Na+ and K+ in the gas diffusion. In Fig. 5 the force curves for all the systems are plotted up to a point-to-plane distance of 4 Å. Within the parentheses the energy barrier for each respective force curve is presented. We have not calculated the force data on the other end of the cage since the Na+ loaded zeolite showed that the calculations are essentially symmetric. This figure clearly demonstrates the role of the extra framework cations in the 8R pore window and strongly supports the idea of the possibility of fine-tuning the functionality of a material as a molecular sieve through ion exchange. The difference in energy barriers for CO2 and N2 is apparent. Also the simulations clearly demonstrate the difference in the blocking efficiency of Na+ and K+ ions. In the case of a cation-free pore window the energy barrier is negligible, around 4 kJ mol1, leaving a, more or less, free range for CO2 and N2 to pass through the 8R pore window. The presence of the Na+ ion in the pore window increases the energy barrier with 17 kJ mol1 and 23 kJ mol1 for CO2 and N2, respectively, while the presence of the K+ ion on the other hand increases the energy barrier with 41 kJ mol1 and 44 kJ mol1 for respective gases. From the free energy barriers we determine the transition state relations between the various gas–cation combinations, which we wish to compare. We calculate the relative transition rate constants from the Arrhenius equation according to, ka Aa ðEa Eb Þ kT ¼ e kb A b
(1)
where the collision rate factors Aa and Ab refer to that of respective gases. These relative transition rate constants calculated for each gas–ion combination are displayed in Table 1. The gas phase N2-to-CO2 collision rate ratio of 1.05 has been taken into account. From Table 1 we see relative Na+-over-K+ transition rate constants in the range of 104 and 107 for K+-over-empty pore.
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 169
View Article Online
PCCP
Paper
Table 1 Relative transition rate constants ka/kb for gas–ion combinations a and b
a +
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
CO2–Na CO2–K+ CO2–Na+ N2–Na+ N2–K+ N2–Na+ CO2–no ion CO2–Na+ CO2–K+
b
ka/kb [d.u.]
CO2–no ion CO2–no ion CO2–K+ N2–no ion N2–no ion N2–K+ N2–no ion N2–Na+ N2–K+
1.0 2.0 1.7 1.7 5.1 5.0 1.1 11 3.5
102 107 104 104 107 104
Thus the rate of gas molecules passing through a K+-containing 8R is negligible. This is in line with the conclusion drawn from the results presented by Liu et al.1 and Mace et al.,4 where the K+ ion effectively blocks the passage of CO2 and N2 through the pores and, in turn, hinders the uptake of each respective gas. When further analysing the shape of the force curves, there is another prominent distinction between the ion types; the K+ force curves begin to increase approximately 1 Å earlier than corresponding curves of Na+, around 4 Å point-to-plane distance compared to 3 Å for Na+. This can be explained at least partly by the fact that the K+ ion is larger and the interactions between the gas molecule and the K+ ion begins at a longer point-to-plane distance compared to Na+. However the van der Waals radii for the two ions differ by approximately 0.5 Å, as the radius is 2.3 Å for Na+ and 2.8 Å for K+, so a hard sphere model is not sufficient to explain this difference in the force curves. The remaining difference can be explained by the fact that the Na+ ion is situated more off center in the 8R pore window plane compared to the larger K+ ion. This will allow the gas molecule to come closer to the plane in the case of Na+, before the gas molecule– ion interaction begins. From the shape of the force curve it is obvious that the K+ ion results in a drastic increase of the diffusion barriers
Fig. 7 Visualization of the spatial behaviour of Na+ and K+, respectively, as the CO2 molecules is steered through the 8R pore window where the K+ ion is pushed out of the 8R plane while the Na+ ion is not.
compared to the case of the Na+ ion. The approximated diffusion barrier for N2 is 48 kJ mol1 and the CO2 barrier is close to that, 45 kJ mol1. These barriers are very high, and the diffusion will be very slow. Besides the much higher barriers the maxima are also shifted. For N2 the maximum is at 1 Å from the pore and for CO2 1.5 Å. A more detailed analysis of the kinetic mechanism reveals interesting details. In the case of Na+ the N2 and CO2 will pass through the pore without much disturbance to the position of the Na+ ion but in the case of K+ the ion will clearly be forced to depart from the 8R plane. The average deviations of the ions are shown in Fig. 6. Fig. 7 shows three snapshots of the CO2 passage top to bottom at 4 Å, 2 Å and 0 Å point-to-plane distance for Na+ on the left and K+ on the right. The fact that K+ is displaced but Na+ is not significantly affected is in agreement with the different sizes of these ions.
Interpretation and conclusions
Fig. 6 The mean ion displacements from the 8R-plane are plotted as a function of the point-to-plane distance. It is apparent from the graphs that both CO2 and N2 molecules are able to squeeze through the 8R when occupied by a Na+ ion opposed to the case of K+, which is forced out of the ring due to its large size.
170 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
This AIMD approach, to compute and compare the free energy barriers for the passage of both CO2 and N2 through Na+ and K+ blocked pore windows of zeolite A, has been successful in producing stable and reliable results. This approach has allowed us to investigate the effect of differential molecular sieving through ion exchange without the bias of empirical parameterization such as the case for empirical simulation methods. The results give strong support to the idea of a
This journal is © the Owner Societies 2014
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Paper
tuneable sieving capability through ion exchange in zeolite A as presented by Liu et al.,1 following experimental results. The computed free energy barriers have a strong qualitative correlation with the experimental results, showing Na+-to-K+ relative transition rate constants in the magnitude of 104, resulting in a negligible probability of one of the gas molecules passing through an 8R pore containing a K+ ion. The experiments show zero uptakes for both N2 and CO2 for zeolite KA as the pore windows are completely blocked for both N2 and CO2 molecules to pass through and enter the material. Our results are further supported by classical MD results,4 which show very small diffusion coefficients for CO2 in zeolite KA as compared to zeolite NaA where the simulated trajectories show that CO2 does not pass pores, which contain a K+ ion. Further, the higher energy barriers for N2 as compared to CO2 result in a CO2-to-N2 relative transition rate constant of 11, when passing a Na+-containing pore, as presented in Table 1. This number corresponds with the CO2-over-N2 selectivity measured experimentally for zeolite NaA by Liu et al.1 We evaluated the energy barriers for N2 and CO2 passing through the zeolite pore window. Due to the high computational cost of the AIMD simulations, the direct simulation of the diffusion coefficients of gas molecules in the zeolite NaKA framework was not feasible with AIMD. In the separate empirical MD work presented by us (Mace and Laaksonen4), we showed that a considerable contribution to the decrease in the CO2 diffusion, as a function of increasing K+ content, could be attributed to K+ having a sufficiently lower mobility than the comparably mobile Na+ ions. In the simulations, several spontaneous Na+ jumps between the equivalent sites within the 8R pore windows were observed, which are not observed for K+ as it sits in the centre of the 8R. Further, even in this short time period, a couple of spontaneous Na+ jumps between pore windows were also observed, which conforms to the observation of these occurrences from the classical simulations. We did not observe any K+ jump between the pores in the AIMD simulations. However, no conclusions can be drawn from this observation as there is only one K+-ion in the system, hence, the statistics are very poor. This work has provided additional information bringing us closer to understanding the full dynamics involved in what has shown to be a quite complex separation process where both thermodynamics and kinetics have significant contributions. Firstly, the thermodynamic selectivity originates in the fact that CO2 possesses a quadrupolar moment approximately three times that of N2, 14 1040 C m1 and 4.7 1040 C m1 respectively.32 This results in a CO2-over-N2 single component selectivity in the range of 10–20 as measured by GCMC simulations4 and can be compared with the single component selectivities of similar materials without size restricting pores measured experimentally.33–38 This also corresponds with the CO2-over-N2 binding energy difference of approximately 15 estimated from the AIMD data. Secondly, there is the kinetic contribution, the molecular sieving effect optimized through ion exchange, which is the source to the exceptionally high CO2-over-N2 selectivity making this material remarkable. Further in this work together with the
This journal is © the Owner Societies 2014
PCCP
separate MD work4 it is apparent that kinetic contribution can, in turn, be separated into two different dynamical processes. On the one hand, the difference in mobility of the ions as described by MD,4 and on the other hand, a contributing size effect, as demonstrated by this work. Hence, the larger K+ ion neither allows the gas molecule to pass nor does it diffuse away from the 8R position. While Na+ has a smaller size allowing the gas molecule to squeeze through a Na+ situated 8R without disturbing the position of the ion significantly. This size effect, together with the previous observation that Na+ has a comparably high mobility,4 explains the facilitated diffusion and uptake for zeolite NaA as compared to KA. The observation that a K+ ion situated in the 8R functions as a block, hindering the passage of either of the gas molecules should essentially be seen as a perturbation of the diffusion of the gas molecules in the system. At low K+ coverage (0–15 at%) the few 8R pore windows blocked by K+ ions do not strongly disturb the diffusion of the gas molecules as many paths still are open, hence the uptake is not noticeably decreased. This corresponds to the experimental uptake plateau for both CO2 and N2 presented by Liu et al.1 At medium K+ coverage the pore blocking starts to have a substantial effect on the diffusional possibilities and in turn decreases the uptake. At high K+ coverage the diffusion is almost blocked as the probability of finding a passable 8R is vastly decreased. The presented procedure for using constrained AIMD to measure energy barriers for gas molecules in porous materials is to our knowledge first of its kind. It is clear that the computations have provided solid data containing important and unbiased information allowing us to quantitatively differentiate the energy barriers for the different cases presented. Further, the results are fully in line with our own hypothesis as well as the widely acknowledged view that CO2 and N2 gases can pass an 8R window in zeolite A containing Na+ but not K+ due to the differences in size. Hence, this procedure should be fully applicable to other zeolites or similar porous materials where the effect on the energy barriers can be measured for a wider range of extra-framework ions, varying ring sizes, altered framework atom composition, framework rigidity or different small adsorbent molecules just to name a few possibilities. Mapping how these and other properties affect the diffusion and molecular sieving is likely to play an important role when finding candidates for potential materials for the purpose of gas separation. Further, energy barrier computations of this type enable us to determine accurately the transition rates. This information could potentially be used as input for Transition State Theory based coarse-graining methods such as kinetic Monte Carlo.39–42 If the transition state landscape is successfully parameterized, such methods could enable accurate diffusion simulations on up to a microsecond timescale and longer, even for N2, without the bias of empirically parameterized force fields.
Notes and references 1 Q. Liu, A. Mace, Z. Bacsik, J. Sun, A. Laaksonen and N. Hedin, Chem. Commun., 2010, 46, 4502. 2 J. J. Pluth and J. V. Smith, J. Am. Chem. Soc., 1980, 102, 4704. 3 D. W. Breck, Zeolite Molecular Sieves, Wiley, New York, 1974.
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 171
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
PCCP
4 A. Mace, N. Hedin and A. Laaksonen, J. Phys. Chem. C, 2013, DOI: 10.1021/jp4048133. 5 P. M. Zimmerman, D. C. Tranca, J. Gomes, D. S. Lambrecht, M. Head-Gordon and A. T. Bell, J. Am. Chem. Soc., 2012, 134, 19468. ´ngya ´n, J. Catal., 2011, 6 T. Bucˇko, L. Benco, J. Hafner and J. A 279, 220. 7 M. De Toni, R. Jonchiere, P. Pullumbi, F. Coudert and A. H. Fuchs, ChemPhysChem, 2012, 13, 3497. 8 D. P. Kling, E. S. A. Machado, H. C. Chagas, A. P. A. dos Santos, N. Rosenbach, J. Walkimar Carneirod and C. J. A. Mota, Chem. Commun., 2013, 49, 4480. 9 J. Hafner, J. Phys.: Condens. Matter, 2010, 22, 384205. 10 A. Pavlova, T. T. Trinh, R. A. van Santen and E. J. Meijer, Phys. Chem. Chem. Phys., 2013, 15, 1123. ´nchez, C. O. Ania, J. B. Parra, D. Dubbeldam, 11 A. Garcı´a-Sa T. J. H. Vlugt, R. Krishna and S. Calero, J. Phys. Chem. C, 2009, 113, 8814. 12 L.-C. Lin, A. H. Berger, R. L. Martin, J. Kim, J. A. Swisher, K. Jariwala, C. H. Rycroft, A. S. Bhown, M. W. Deem, M. Haranczyk and B. Smit, Nat. Mater., 2012, 11, 633. 13 G. Maurin, P. L. Llewellyn and R. G. Bell, J. Phys. Chem. B, 2005, 109, 16084. 14 E. Jaramillo and M. Chandross, J. Phys. Chem. B, 2004, 108, 20155. 15 E. D. Akten, R. Siriwardane and D. S. Sholl, Energy Fuels, 2003, 17, 977. 16 R. Krishna and J. M. van Baten, Microporous Mesoporous Mater., 2009, 125, 126. 17 R. Krishna and J. M. van Baten, Microporous Mesoporous Mater., 2011, 137, 83. 18 L. Madison, H. Heitzer, C. Russell and D. Kohen, Langmuir, 2011, 27, 1954. 19 M. Sant, J.-M. Leyssale, G. K. Papadopoulos and D. N. Theodorou, J. Phys. Chem. B, 2009, 113, 13761. 20 M. R. Hudson, W. L. Queen, J. A. Mason, D. W. Fickel, R. F. Lobo and C. M. Brown, J. Am. Chem. Soc., 2012, 143, 1970. ´nchez, J. van den Bergh, J. M. Castillo, S. Calero, 21 A. Garcı´a-Sa F. Kapteijn and T. J. H. Vlugt, Microporous Mesoporous Mater., 2012, 158, 12.
172 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Paper
22 www.cp2k.org (visited 26.11.2012). 23 G. Lippert, J. Hutter and M. Parinello, Mol. Phys., 1997, 92, 477. 24 S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703. 25 S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104. 26 M. Sprik and G. Ciccotti, J. Chem. Phys., 1998, 109, 7737. 27 E. A. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Chem. Phys. Lett., 1989, 156, 472. 28 www.accelrys.com/products/materials-studio (visited 22/01/ 2013). 29 V. Gramlich and W. M. Meier, Z. Kristallogr., Kristallgeom., Kristallphys., Kristallchem., 1971, 133, 134. 30 V. Cerny´, J. Optimiz. Theory App., 1985, 45, 41. 31 S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Science, 1983, 220, 671. 32 C. Graham, D. A. Imrie and R. E. Raab, Mol. Phys., 1998, 93, 49. ´ndes and J. Santamarı´a, 33 M. P. Bernal, M. Coronas, M. Mene AIChE J., 2004, 50, 127. 34 R. V. Sirwardane, M. S. Shen and E. Fisher, Energy Fuels, 2005, 19, 1153. 35 T. Q. Gardner, A. I. Flores, R. D. Noble and J. L. Falconer, AIChE J., 2002, 48, 1155. 36 P. J. Harlick and F. H. Tezel, Sep. Purif. Technol., 2003, 33, 199. 37 T. Hasegaw, K. Watenabe, K. Kusakabe and S. Morooka, Sep. Purif. Technol., 2001, 22–23, 319. 38 K. Kusakabe, T. Kuroda, K. Uchino, Y. Hasegawa and S. Marooko, AIChE J., 2004, 45, 1220. 39 C. Tunca and D. M. J. Ford, Chem. Phys., 1999, 111, 2751. 40 M. K. F. Abouelnasra and B. Smit, Phys. Chem. Chem. Phys., 2013, 14, 11600. 41 C. Saravanan and S. J. Auerbach, Chem. Phys., 1999, 110, 11000. 42 R. Krishna and J. M. van Baten, Chem. Eng. Technol., 2005, 28, 160.
This journal is © the Owner Societies 2014
PAPER
Cite this: Phys. Chem. Chem. Phys., 2014, 16, 166
View Journal | View Issue
Free energy barriers for CO2 and N2 in zeolite NaKA: an ab initio molecular dynamics approach Amber Mace,a Kari Laasonenb and Aatto Laaksonen*a
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Ab initio Molecular Dynamics (AIMD) is used with spatial constraints to estimate the free energy barriers of diffusion for CO2 and N2 gas molecules in zeolite NaA and KA. We investigate the extent to which the diffusion of these gas molecules is hindered, in the two separate cases of a smaller Na+ ion or a larger K+ ion blocking the 8-ring pore window. In contrast to classical Molecular Dynamics, AIMD performs these computations accurately and unbiased in the absence of empirical parameterization. Our work has resulted in stable and reliable force profiles. The profiles show that the larger K+ ion effectively blocks Received 5th July 2013, Accepted 25th October 2013
the passage of both CO2 and N2 molecules while the smaller Na+ ion will allow both molecules to pass.
DOI: 10.1039/c3cp52821a
effect, which the size of the respective cation occupying the pore window has on diffusive properties of
www.rsc.org/pccp
each gas molecule. Hence, this effect can be altered through ion exchange to fine-tune the functionality of a specific zeolite as a molecular sieve.
These results are a quantitative demonstration of the concept of pore blocking where we compute the
Introduction Certain zeolites with narrow window apertures, such as those of zeolite A (framework type code LTA), coincide with the approximate size of small gaseous molecules such as CO2 and N2. These may be interesting candidates for adsorbents with swing adsorption technologies due to their potential molecular sieving capabilities and otherwise attractive properties. Zeolite A is a threedimensional cage type microporous aluminosilicate with an Al-toSi ratio of 1 resulting in a high cation content, with respect to the zeolite framework, which acts to neutralize the framework. This gives the framework neutralized by Na+ and K+ ions (zeolite NaKA), the chemical composition Na96xKxAl96Si96O394, 0 o x o 96. With the point group FM3c it has a cubic periodicity built up by so-called a- and b-cages where one unit cell consists of eight (2 2 2) a-cages as shown in Fig. 1. The a-cages intersect by 8-ring (8R) pore windows, where the number of oxygen atoms within the ring, confining the pore window, determines the window ring number. The framework also contains 6-ring (6R) and 4-ring (4R) pore windows, however it is only the 8R pore windows, which are large enough to allow the passage CO2 or N2. As shown by PXRD measurements,1,2 zeolite A has three different cation sites; site I coordinated with the centre of the 6R, site II in the plane of the 8R and site III coordinated with the centre of the 4R. When zeolite A solely a
Arrhenius Laboratory, Stockholm University, 10691 Stockholm, Sweden. E-mail: [email protected]; Fax: +46-8-152187; Tel: +46-8-162372 b Department of Chemistry, Aalto University, P.O. Box 16100, 00076 Espoo, Finland. E-mail: [email protected]; Tel: +358-40-5570044
166 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Fig. 1 An LTA unit cell (left) is built up by eight a-cages (bottom right). The smaller b-cages (top right) are comprised in the cubic intersect of the a-cages.
contains monovalent extra-framework cations, each of the 64 6R and 24 8R per unit cell will generally situate one cation, while the remaining eight cations are distributed among the 4R pore windows. Furthermore, the Na+ and K+ ions have slightly different positions relative to the pore window rings as illustrated in Fig. 2. Site II for K+ is centred in the 8R, while the smaller Na+ is approximately 1.3 Å off-centre.2 This off-centre placement results in four energetically and geometrically equivalent sites for Na+ in the 8R plane. Zeolite A neutralized with Na+ cations is referred to as zeolite NaA or 4A and with K+, zeolite KA or 3A. Here 4A and 3A denote the resulting pore window diameters of approximately 4 Å and
This journal is © the Owner Societies 2014
View Article Online
Paper
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Fig. 2
PCCP
Crystallographic K+ and Na+ positions in zeolite A.
3 Å respectively. This can be compared to the kinetic diameters of CO2 and N2, which are 3.3 Å and 3.6 Å, respectively.3 The kinetic diameter refers to the smallest effective dimension of the molecules, the critical molecular dimension. This dimension is determined from the separation distance, which corresponds to the minimum equilibrium separation of the Lennard-Jones potential. Hence, the stronger attractive forces of CO2 can rationalize this somewhat counter-intuitive relation that N2 is, in this aspect, larger than CO2. Following this, both CO2 and N2 should be able to pass through the 8R pore windows accommodating Na+, while a K+ ion in the same position will block the passage. In the recently published paper by Liu et al.,1 particularly appealing results were presented showing that through Na+-toK+ ion exchange in zeolite 4A an exceptionally high CO2-over-N2 selectivity could be reached while the total CO2 uptake still remained high at a specific K+/(K+ + Na+) ratio of 17 at%. Liu et al.1 predicted that this high selectivity was due to a sieving effect as a result of an effective pore diameter lying right in between the kinetic diameters of CO2 and N2 allowing the material to adsorb CO2 while effectively blocking out N2. In a separate work reported by two of the current authors,4 classical simulation methods have been utilized in order to model the separate contributions of thermodynamics and kinetics involved in the gas uptake in zeolite NaKA, and how these contributions are affected by the two different ion species. By combining results from classical Molecular Dynamics (MD) and Grand Canonical Monte Carlo (GCMC) simulations the experimental results for CO2 were qualitatively reproduced, indicating that the dynamics behind the high selectivity measured by Liu et al.1 was, in fact, a combination of these contributions. In this work we attempt to go further and look in detail at the dynamics involved, when the gas molecules diffuse from pore-to-pore, by using ab initio Molecular Dynamics (AIMD). For zeolites and other microporous materials, AIMD methods have been mainly used when investigating catalysis and chemisorption of guest molecules in the frameworks.5–10 While, typically, classical simulation methods such as GCMC and classical MD are used to model the physisorbed equilibrium uptakes11–15 and diffusion coefficients16–20 of gas molecules in zeolitic systems. These classical methods use predetermined force fields to model each inter-atomic interaction and are
This journal is © the Owner Societies 2014
commonly parameterized to empirical data sets, such as adsorption isotherms. A common procedure when fitting empirical force field parameters for zeolite–gas interactions is to use GCMC simulations to match an experimentally measured adsorption isotherm. This approach may be a good estimation when investigating gas uptake with zeolites without size restricting pores. However, when an unrestricted pore-to-pore diffusion of the adsorbent is, to some extent, hindered by the size or cation blocking of the pore windows, such as in the case of zeolite A, the uptake will be affected. In this case, determining a force field by appointing the full uptake to thermodynamics will create a bias. Also combinations of different force field parameters may give the same outcome, which is apparent when comparing the fairly large amount of zeolite NaA–CO2 force fields with widely varying parameters for the same interactions.11,13–15 This may not have a great effect on a larger scale, but when modelling a specific interaction, the outcome may vary greatly between different force fields; hence the results are not accurate. Further, recent work by Garcı´a´nchez et al.21 showed that the choice of force field does, in Sa fact, strongly affect the outcome of diffusion simulations. With this approach based on periodic density functional theory (DFT) and empirical van der Waals corrections an accurate and unbiased description of the molecular interactions can be obtained. We can investigate how the forces act on the different gas molecules when moving through the pores as well as what effect the cation type has on this. Also the energy barriers for the passage through the pores can be computed. The DFT does not need specific parameterization for this system as is the case for classical molecular dynamics with empirical force fields and thus the DFT can really predict the results. In addition, it is easy to investigate the dynamics of the N2 molecules compared with CO2, which is problematic with classical MD due to the very slow diffusion of N2 molecules in the framework. The down side of this approach is the high computational cost forcing us to narrow down the computations to a smaller simulation cell and shorter simulation length compared to what is possible for classical MD. Also, it is necessary to use constraints in order to look at the specific dynamics of interest. From this point of view, this work is well needed and can be seen as pioneering in the aspect of using AIMD to investigate the dynamics of physisorbed guest molecules in microporous materials.
Computational details For the AIMD simulations, the CP2K software package is utilized. CP2K performs Born Oppenheimer Dynamics, which calculates DFT forces ‘‘on the fly’’ with the QUICKSTEP22 module using the GPW method.23 Here, the atom centred Gaussian basis sets are used together with augmented plane wave basis sets. For the DFT force calculations we used the PBE functional and the DZVP-MOLOPT-SR basis set together with GTH pseudopotential for all atom types.24 For all atoms, only the valence electrons were treated explicitly with the DFT. Both
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 167
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
PCCP
Na+ and K+ have only one explicit valence electron. The plane wave kinetic energy cutoff was set to 280Ry. Also, the empirical vdW correction, DFT-D3, by Grimme25 was used. The simulations were run with a 1.5 fs time-step where the equations of motion were propagated using the Velocity-Verlet algorithm. Further, the simulations were run in an NVT ´–Hoover thermostat with the ensemble employing the Nose temperature set to 298.15 K. The free energy was calculated by using the constraint method.26 This method allows a direct estimation of the mean force from the time average of the force on the geometric constraint by applying the blue-moon ensemble method.27 The free energy is then estimated by integrating the mean force curve. The used constraint will be explained later. The zeolite NaA structure was created in the Accelrys Materials Studio Suite.28 The coordinates were taken from the structure database and were originally presented by Gramlich and Meier.29 The Sorption module was used to position the 96 extra-framework Na+ ions per unit cell in the framework using a Canonical Monte Carlo ensemble simulated annealing energy minimizing procedure.30,31 As a simulation cell in the AIMD computations 1/4 of a unit cell was used, corresponding to two intersecting a-cages on the x-axis, as shown in Fig. 3, and periodic boundary conditions were implemented with cell parameters a, b, c = 24.6 Å, 12.3 Å, 12.3 Å. One gas molecule is placed in the centre of the left a-cage. As we attempt to measure and compare the energy barriers for the gas molecules to pass through the pore window from one a-cage to another we use a constraint to steer the gas molecule through the intersecting pore window. This is implemented by using a point-to-plane constraint where the plane is defined by three oxygen atoms in the 8R and the point is one atom in the gas molecule, for CO2 the oxygen atom and for N2 one of the nitrogen atoms. Hence, the entropy contribution comes from the molecular movement around the fixed atom. The point-to-plane distance is visualized in Fig. 3. With this constraint the gas molecule is free to move within a plane at the defined distance from the window plane in the AIMD simulation run as well as rotate around the constrained point-atom.
Fig. 3 The simulation cell consists of two of the zeolite A building units, a-cage. One gas molecule (here CO2) is placed in the left a-cage and a constraint is set confining the distance between the oxygen of the CO2 molecule and the left 8R window plane.
168 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Paper
Simulation runs are set up with a subsequently decreasing point-to-plane distance. Near the force maxima a typical distance step was 0.25 Å and around the centre of the cage the distance step was 1.0 Å. For each constraint value an AIMD simulation of length of at least 10 ps is run to ensure good convergence of the constraint force, since the force has large fluctuations. To check the stability of the system the same procedure was performed in the opposite direction, where the gas molecule was steered towards the 8R in the centre of the simulation cell. The next step was then to look at how the ion exchange affects the system. We exchanged the Na+ ion in the left 8R in the simulation cell for a K+ ion and proceeded to perform the simulations in the same manner as for the previous system for both the CO2 and N2 molecules. Finally, for a comparison the same simulations were performed in the system where the ion was removed completely from the 8R pore window.
Results and discussion The AIMD computed forces for each distance step simulation run are presented as a function of the constrained point-toplane distance in Fig. 4 and 5. In Fig. 4 the force data are presented for the movement of CO2 and N2 molecules through the whole length of the Na+ containing zeolite cage, from the left 8R to the centre 8R. This picture displays the clear force differences for the 8R-approaching CO2 compared with N2 where the latter experiences a larger repulsive force. This is in line with the hypothesis that it is more difficult for the larger N2 molecules to pass through an 8R, which situates an ion. The energy barriers are calculated by integrating the forces in the point-to-plane distance range 0–4 Å where the approximated barrier for the N2 diffusion is 27 kJ mol1 and 21 kJ mol1 for CO2, giving a barrier difference of 6 kJ mol1. From the Arrhenius equation, this will correspond to a transition state relation of around 1 : 11. From the trajectories of CO2 and N2 at point-to-plane distances o1 Å, we observe
Fig. 4 The absolute values of the averaged forces are, for CO2 and N2 in the pure Na+ zeolite A, plotted as a function of the distance between the 8R window plane and the COM of the gas molecule.
This journal is © the Owner Societies 2014
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Paper
Fig. 5 The absolute values of the averaged forces for CO2 and N2 in the pure Na+ zeolite A as well as for the K+ ion exchange structure are plotted as a function of the constrained distance between the 8R window plane and Oxygen of CO2 and one N-atom of N2. The energy barrier values are presented in the parenthesis and are defined by the integral of the force in the range 0–4.0 Å point-to-plane distances.
differences in the behaviour of the gas molecules that can help to rationalize the difference in the energy barriers. N2 tends to point towards the cation, resulting from an induced dipole interaction with the Na+ ion. This effect of strong steering counteracts the N2 molecule from penetrating the pore. For the CO2, on the other hand, the molecule interacts both with the Na+ and the cage. The CO2 is thus oriented more perpendicular with the pore plane due to the balancing of CO2 interactions between the Na+ ion and the cage. Overall, the CO2 is affected less than N2 by the Na+ ion, and this will result in a decreased transition state for CO2 to pass through the pore window. The other notable difference between the CO2 and N2 curves is the behaviour in their centre regions. The CO2 potential is not zero at the centre of the zeolite cage, which is due to the stronger interaction of the CO2 molecule with the zeolite cage. This stronger interaction is due to the larger quadrupole moment of CO2 and the bending vibrations of CO2, causing the molecule to have an instantaneous dipole moment. Even when the CO2 molecule is in the centre of the pore the environment is not exactly symmetric, as the Na+ ion is situated off the centre of the 8R and the centred CO2 molecule generally attracts an 8R ion in the xy or the xz plane. Hence, the molecule will be drawn more to the left or the right depending on which side of the 8R the Na+ is located. This causes the non-zero force even for the CO2 around the 6 Å position. The force integral in the centre region of the zeolite cage (3–9.6 Å) is 19 kJ mol1 for CO2 compared to 6.0 kJ mol1 for N2. This gives us a crude estimation of the binding energy difference between CO2 and N2. The real force curve is antisymmetric, so the numbers above need to be divided by 2 and the energy difference will be ca. 6.5 kJ mol1. At room temperature this will lead to an adsorption selectivity of 1 : 15 favouring CO2. This is a good agreement with the GCMC simulations of Mace
This journal is © the Owner Societies 2014
PCCP
et al.,4 where the selectivity was in the range of 1 : 10–1 : 20, depending on the K+ content. Another important observation from Fig. 4 is that the force curves are symmetric and the data scattering is relatively small. This indicates that the simulations are long enough to yield well-converged forces with error bars varying between 0.2–1 103 a.u. These force fluctuations are especially large around 2–3 Å from the planes (0.5–1 103 a.u.), hence the points at distances 2, 2.5, 3, 9, 10 and 10.5 Å. This can be explained by the gas molecule being at a critical distance from the pore where the rotations have a large effect on the force and the difference between the gas molecules being either parallel or perpendicular with the plane is quite large. As the distance to the plane decreases, the gas molecule stops rotating since the space becomes more constricted, and in the point-to-plane distance region 0–2 Å the force fluctuations are comparably small (0.2–0.3 103 a.u.). If the distance to the plane increases, the molecules will rotate, but will not cause large fluctuations of the force (0.3–0.5 103 a.u.) as the interaction with the plane is smaller. The next set of simulations focuses on the role of Na+ and K+ in the gas diffusion. In Fig. 5 the force curves for all the systems are plotted up to a point-to-plane distance of 4 Å. Within the parentheses the energy barrier for each respective force curve is presented. We have not calculated the force data on the other end of the cage since the Na+ loaded zeolite showed that the calculations are essentially symmetric. This figure clearly demonstrates the role of the extra framework cations in the 8R pore window and strongly supports the idea of the possibility of fine-tuning the functionality of a material as a molecular sieve through ion exchange. The difference in energy barriers for CO2 and N2 is apparent. Also the simulations clearly demonstrate the difference in the blocking efficiency of Na+ and K+ ions. In the case of a cation-free pore window the energy barrier is negligible, around 4 kJ mol1, leaving a, more or less, free range for CO2 and N2 to pass through the 8R pore window. The presence of the Na+ ion in the pore window increases the energy barrier with 17 kJ mol1 and 23 kJ mol1 for CO2 and N2, respectively, while the presence of the K+ ion on the other hand increases the energy barrier with 41 kJ mol1 and 44 kJ mol1 for respective gases. From the free energy barriers we determine the transition state relations between the various gas–cation combinations, which we wish to compare. We calculate the relative transition rate constants from the Arrhenius equation according to, ka Aa ðEa Eb Þ kT ¼ e kb A b
(1)
where the collision rate factors Aa and Ab refer to that of respective gases. These relative transition rate constants calculated for each gas–ion combination are displayed in Table 1. The gas phase N2-to-CO2 collision rate ratio of 1.05 has been taken into account. From Table 1 we see relative Na+-over-K+ transition rate constants in the range of 104 and 107 for K+-over-empty pore.
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 169
View Article Online
PCCP
Paper
Table 1 Relative transition rate constants ka/kb for gas–ion combinations a and b
a +
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
CO2–Na CO2–K+ CO2–Na+ N2–Na+ N2–K+ N2–Na+ CO2–no ion CO2–Na+ CO2–K+
b
ka/kb [d.u.]
CO2–no ion CO2–no ion CO2–K+ N2–no ion N2–no ion N2–K+ N2–no ion N2–Na+ N2–K+
1.0 2.0 1.7 1.7 5.1 5.0 1.1 11 3.5
102 107 104 104 107 104
Thus the rate of gas molecules passing through a K+-containing 8R is negligible. This is in line with the conclusion drawn from the results presented by Liu et al.1 and Mace et al.,4 where the K+ ion effectively blocks the passage of CO2 and N2 through the pores and, in turn, hinders the uptake of each respective gas. When further analysing the shape of the force curves, there is another prominent distinction between the ion types; the K+ force curves begin to increase approximately 1 Å earlier than corresponding curves of Na+, around 4 Å point-to-plane distance compared to 3 Å for Na+. This can be explained at least partly by the fact that the K+ ion is larger and the interactions between the gas molecule and the K+ ion begins at a longer point-to-plane distance compared to Na+. However the van der Waals radii for the two ions differ by approximately 0.5 Å, as the radius is 2.3 Å for Na+ and 2.8 Å for K+, so a hard sphere model is not sufficient to explain this difference in the force curves. The remaining difference can be explained by the fact that the Na+ ion is situated more off center in the 8R pore window plane compared to the larger K+ ion. This will allow the gas molecule to come closer to the plane in the case of Na+, before the gas molecule– ion interaction begins. From the shape of the force curve it is obvious that the K+ ion results in a drastic increase of the diffusion barriers
Fig. 7 Visualization of the spatial behaviour of Na+ and K+, respectively, as the CO2 molecules is steered through the 8R pore window where the K+ ion is pushed out of the 8R plane while the Na+ ion is not.
compared to the case of the Na+ ion. The approximated diffusion barrier for N2 is 48 kJ mol1 and the CO2 barrier is close to that, 45 kJ mol1. These barriers are very high, and the diffusion will be very slow. Besides the much higher barriers the maxima are also shifted. For N2 the maximum is at 1 Å from the pore and for CO2 1.5 Å. A more detailed analysis of the kinetic mechanism reveals interesting details. In the case of Na+ the N2 and CO2 will pass through the pore without much disturbance to the position of the Na+ ion but in the case of K+ the ion will clearly be forced to depart from the 8R plane. The average deviations of the ions are shown in Fig. 6. Fig. 7 shows three snapshots of the CO2 passage top to bottom at 4 Å, 2 Å and 0 Å point-to-plane distance for Na+ on the left and K+ on the right. The fact that K+ is displaced but Na+ is not significantly affected is in agreement with the different sizes of these ions.
Interpretation and conclusions
Fig. 6 The mean ion displacements from the 8R-plane are plotted as a function of the point-to-plane distance. It is apparent from the graphs that both CO2 and N2 molecules are able to squeeze through the 8R when occupied by a Na+ ion opposed to the case of K+, which is forced out of the ring due to its large size.
170 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
This AIMD approach, to compute and compare the free energy barriers for the passage of both CO2 and N2 through Na+ and K+ blocked pore windows of zeolite A, has been successful in producing stable and reliable results. This approach has allowed us to investigate the effect of differential molecular sieving through ion exchange without the bias of empirical parameterization such as the case for empirical simulation methods. The results give strong support to the idea of a
This journal is © the Owner Societies 2014
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
Paper
tuneable sieving capability through ion exchange in zeolite A as presented by Liu et al.,1 following experimental results. The computed free energy barriers have a strong qualitative correlation with the experimental results, showing Na+-to-K+ relative transition rate constants in the magnitude of 104, resulting in a negligible probability of one of the gas molecules passing through an 8R pore containing a K+ ion. The experiments show zero uptakes for both N2 and CO2 for zeolite KA as the pore windows are completely blocked for both N2 and CO2 molecules to pass through and enter the material. Our results are further supported by classical MD results,4 which show very small diffusion coefficients for CO2 in zeolite KA as compared to zeolite NaA where the simulated trajectories show that CO2 does not pass pores, which contain a K+ ion. Further, the higher energy barriers for N2 as compared to CO2 result in a CO2-to-N2 relative transition rate constant of 11, when passing a Na+-containing pore, as presented in Table 1. This number corresponds with the CO2-over-N2 selectivity measured experimentally for zeolite NaA by Liu et al.1 We evaluated the energy barriers for N2 and CO2 passing through the zeolite pore window. Due to the high computational cost of the AIMD simulations, the direct simulation of the diffusion coefficients of gas molecules in the zeolite NaKA framework was not feasible with AIMD. In the separate empirical MD work presented by us (Mace and Laaksonen4), we showed that a considerable contribution to the decrease in the CO2 diffusion, as a function of increasing K+ content, could be attributed to K+ having a sufficiently lower mobility than the comparably mobile Na+ ions. In the simulations, several spontaneous Na+ jumps between the equivalent sites within the 8R pore windows were observed, which are not observed for K+ as it sits in the centre of the 8R. Further, even in this short time period, a couple of spontaneous Na+ jumps between pore windows were also observed, which conforms to the observation of these occurrences from the classical simulations. We did not observe any K+ jump between the pores in the AIMD simulations. However, no conclusions can be drawn from this observation as there is only one K+-ion in the system, hence, the statistics are very poor. This work has provided additional information bringing us closer to understanding the full dynamics involved in what has shown to be a quite complex separation process where both thermodynamics and kinetics have significant contributions. Firstly, the thermodynamic selectivity originates in the fact that CO2 possesses a quadrupolar moment approximately three times that of N2, 14 1040 C m1 and 4.7 1040 C m1 respectively.32 This results in a CO2-over-N2 single component selectivity in the range of 10–20 as measured by GCMC simulations4 and can be compared with the single component selectivities of similar materials without size restricting pores measured experimentally.33–38 This also corresponds with the CO2-over-N2 binding energy difference of approximately 15 estimated from the AIMD data. Secondly, there is the kinetic contribution, the molecular sieving effect optimized through ion exchange, which is the source to the exceptionally high CO2-over-N2 selectivity making this material remarkable. Further in this work together with the
This journal is © the Owner Societies 2014
PCCP
separate MD work4 it is apparent that kinetic contribution can, in turn, be separated into two different dynamical processes. On the one hand, the difference in mobility of the ions as described by MD,4 and on the other hand, a contributing size effect, as demonstrated by this work. Hence, the larger K+ ion neither allows the gas molecule to pass nor does it diffuse away from the 8R position. While Na+ has a smaller size allowing the gas molecule to squeeze through a Na+ situated 8R without disturbing the position of the ion significantly. This size effect, together with the previous observation that Na+ has a comparably high mobility,4 explains the facilitated diffusion and uptake for zeolite NaA as compared to KA. The observation that a K+ ion situated in the 8R functions as a block, hindering the passage of either of the gas molecules should essentially be seen as a perturbation of the diffusion of the gas molecules in the system. At low K+ coverage (0–15 at%) the few 8R pore windows blocked by K+ ions do not strongly disturb the diffusion of the gas molecules as many paths still are open, hence the uptake is not noticeably decreased. This corresponds to the experimental uptake plateau for both CO2 and N2 presented by Liu et al.1 At medium K+ coverage the pore blocking starts to have a substantial effect on the diffusional possibilities and in turn decreases the uptake. At high K+ coverage the diffusion is almost blocked as the probability of finding a passable 8R is vastly decreased. The presented procedure for using constrained AIMD to measure energy barriers for gas molecules in porous materials is to our knowledge first of its kind. It is clear that the computations have provided solid data containing important and unbiased information allowing us to quantitatively differentiate the energy barriers for the different cases presented. Further, the results are fully in line with our own hypothesis as well as the widely acknowledged view that CO2 and N2 gases can pass an 8R window in zeolite A containing Na+ but not K+ due to the differences in size. Hence, this procedure should be fully applicable to other zeolites or similar porous materials where the effect on the energy barriers can be measured for a wider range of extra-framework ions, varying ring sizes, altered framework atom composition, framework rigidity or different small adsorbent molecules just to name a few possibilities. Mapping how these and other properties affect the diffusion and molecular sieving is likely to play an important role when finding candidates for potential materials for the purpose of gas separation. Further, energy barrier computations of this type enable us to determine accurately the transition rates. This information could potentially be used as input for Transition State Theory based coarse-graining methods such as kinetic Monte Carlo.39–42 If the transition state landscape is successfully parameterized, such methods could enable accurate diffusion simulations on up to a microsecond timescale and longer, even for N2, without the bias of empirically parameterized force fields.
Notes and references 1 Q. Liu, A. Mace, Z. Bacsik, J. Sun, A. Laaksonen and N. Hedin, Chem. Commun., 2010, 46, 4502. 2 J. J. Pluth and J. V. Smith, J. Am. Chem. Soc., 1980, 102, 4704. 3 D. W. Breck, Zeolite Molecular Sieves, Wiley, New York, 1974.
Phys. Chem. Chem. Phys., 2014, 16, 166--172 | 171
View Article Online
Published on 28 October 2013. Downloaded on 30/10/2014 10:07:49.
PCCP
4 A. Mace, N. Hedin and A. Laaksonen, J. Phys. Chem. C, 2013, DOI: 10.1021/jp4048133. 5 P. M. Zimmerman, D. C. Tranca, J. Gomes, D. S. Lambrecht, M. Head-Gordon and A. T. Bell, J. Am. Chem. Soc., 2012, 134, 19468. ´ngya ´n, J. Catal., 2011, 6 T. Bucˇko, L. Benco, J. Hafner and J. A 279, 220. 7 M. De Toni, R. Jonchiere, P. Pullumbi, F. Coudert and A. H. Fuchs, ChemPhysChem, 2012, 13, 3497. 8 D. P. Kling, E. S. A. Machado, H. C. Chagas, A. P. A. dos Santos, N. Rosenbach, J. Walkimar Carneirod and C. J. A. Mota, Chem. Commun., 2013, 49, 4480. 9 J. Hafner, J. Phys.: Condens. Matter, 2010, 22, 384205. 10 A. Pavlova, T. T. Trinh, R. A. van Santen and E. J. Meijer, Phys. Chem. Chem. Phys., 2013, 15, 1123. ´nchez, C. O. Ania, J. B. Parra, D. Dubbeldam, 11 A. Garcı´a-Sa T. J. H. Vlugt, R. Krishna and S. Calero, J. Phys. Chem. C, 2009, 113, 8814. 12 L.-C. Lin, A. H. Berger, R. L. Martin, J. Kim, J. A. Swisher, K. Jariwala, C. H. Rycroft, A. S. Bhown, M. W. Deem, M. Haranczyk and B. Smit, Nat. Mater., 2012, 11, 633. 13 G. Maurin, P. L. Llewellyn and R. G. Bell, J. Phys. Chem. B, 2005, 109, 16084. 14 E. Jaramillo and M. Chandross, J. Phys. Chem. B, 2004, 108, 20155. 15 E. D. Akten, R. Siriwardane and D. S. Sholl, Energy Fuels, 2003, 17, 977. 16 R. Krishna and J. M. van Baten, Microporous Mesoporous Mater., 2009, 125, 126. 17 R. Krishna and J. M. van Baten, Microporous Mesoporous Mater., 2011, 137, 83. 18 L. Madison, H. Heitzer, C. Russell and D. Kohen, Langmuir, 2011, 27, 1954. 19 M. Sant, J.-M. Leyssale, G. K. Papadopoulos and D. N. Theodorou, J. Phys. Chem. B, 2009, 113, 13761. 20 M. R. Hudson, W. L. Queen, J. A. Mason, D. W. Fickel, R. F. Lobo and C. M. Brown, J. Am. Chem. Soc., 2012, 143, 1970. ´nchez, J. van den Bergh, J. M. Castillo, S. Calero, 21 A. Garcı´a-Sa F. Kapteijn and T. J. H. Vlugt, Microporous Mesoporous Mater., 2012, 158, 12.
172 | Phys. Chem. Chem. Phys., 2014, 16, 166--172
Paper
22 www.cp2k.org (visited 26.11.2012). 23 G. Lippert, J. Hutter and M. Parinello, Mol. Phys., 1997, 92, 477. 24 S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 1703. 25 S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104. 26 M. Sprik and G. Ciccotti, J. Chem. Phys., 1998, 109, 7737. 27 E. A. Carter, G. Ciccotti, J. T. Hynes and R. Kapral, Chem. Phys. Lett., 1989, 156, 472. 28 www.accelrys.com/products/materials-studio (visited 22/01/ 2013). 29 V. Gramlich and W. M. Meier, Z. Kristallogr., Kristallgeom., Kristallphys., Kristallchem., 1971, 133, 134. 30 V. Cerny´, J. Optimiz. Theory App., 1985, 45, 41. 31 S. Kirkpatrick, C. D. Gelatt and M. P. Vecchi, Science, 1983, 220, 671. 32 C. Graham, D. A. Imrie and R. E. Raab, Mol. Phys., 1998, 93, 49. ´ndes and J. Santamarı´a, 33 M. P. Bernal, M. Coronas, M. Mene AIChE J., 2004, 50, 127. 34 R. V. Sirwardane, M. S. Shen and E. Fisher, Energy Fuels, 2005, 19, 1153. 35 T. Q. Gardner, A. I. Flores, R. D. Noble and J. L. Falconer, AIChE J., 2002, 48, 1155. 36 P. J. Harlick and F. H. Tezel, Sep. Purif. Technol., 2003, 33, 199. 37 T. Hasegaw, K. Watenabe, K. Kusakabe and S. Morooka, Sep. Purif. Technol., 2001, 22–23, 319. 38 K. Kusakabe, T. Kuroda, K. Uchino, Y. Hasegawa and S. Marooko, AIChE J., 2004, 45, 1220. 39 C. Tunca and D. M. J. Ford, Chem. Phys., 1999, 111, 2751. 40 M. K. F. Abouelnasra and B. Smit, Phys. Chem. Chem. Phys., 2013, 14, 11600. 41 C. Saravanan and S. J. Auerbach, Chem. Phys., 1999, 110, 11000. 42 R. Krishna and J. M. van Baten, Chem. Eng. Technol., 2005, 28, 160.
This journal is © the Owner Societies 2014
