Article pubs.acs.org/JPCC
Prediction of Na+/NH4+ Exchange in Faujasite Zeolite by Molecular Dynamics Simulation and Thermodynamic Integration Method Lin Wang and Huai Sun* School of Chemistry and Chemical Engineering, and Ministry of Education Key Laboratory of Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai, 200240, China ABSTRACT: Molecular dynamics simulations and thermodynamic integration method are used to calculate free energy differences when replacing Na+ with NH4+ in hydrated faujasite (Y-type) zeolite and in an aqueous solution of these cations. Force field parameters are optimized by quantum mechanics density functional theory calculations and then validated by calculating cation distributions in dehydrated zeolite and differences in hydration free energies of the cations. Using force field and Monte Carlo simulations, we predict cation distributions at different sites of hydrated zeolite with different mole fractions of NH4+. The free energy changes upon replacement of Na+ with NH4+ at different sites are calculated based on the predicted distributions. The statistically weighted average of the free energy differences of Na+/NH4+ exchange in zeolite is compared with that in an aqueous solution of the cations. The equilibrium isotherm is predicted and found to agree well with experimental data. transition metal and NH4+ ions in Na−clinoptilolite, and their results showed a good selectivity of NH4+. Lim et al.11 investigated the dehydration of NH4+-exchanged zeolite Y and obtained NH4+ distributions. Although the reported ionexchange profiles are not exactly the same because of differences in experimental conditions, all of these aforementioned results indicate that Na−Y zeolite exhibits a high selectivity for NH4+ ion exchange. However, complete NH4+ exchange in Na−Y zeolite cannot be directly achieved. Predicting ion exchanges in zeolites is a challenging task for computational studies. Ion exchange is essentially a chemical equilibrium process governed by free energy changes. Predicting free energy change is usually difficult because it is much smaller than the total free energies of relevant states. Calculation must be sufficiently accurate to reach a precision of roughly 1.5 silicon/aluminum ratio (Si/Al) that has excellent structural stability as well as a large and accessible pore volume. Thus, zeolite Y has wide-ranging industrial applications, including ion exchangers, sorption agents, molecular sieves, and catalysts.1−6 For these applications, Na+/NH4+ exchange is an important industrial process. Synthesized zeolite Y is usually in its sodium form (Na−Y). To produce its hydrogen form (H−Y), ammonium solution is commonly used to replace Na+ with NH4+. Then, water and NH3 are removed by heating (deamination). However, completely replacing Na+ with NH4+ is difficult and an excessive amount of ammonia is usually used, which has serious environmental consequences. Thus, researchers are seeking improved technologies to reduce the usage of ammonia. To design more effective cationexchange agents, the exchange mechanisms must be understood and the exchange isotherms must be predicted. The isotherms of Na+/NH4+ exchange in different zeolites have been experimentally measured.7−10 Sherry et al.7 presented ion-exchange isotherms describing the exchange of alkali metal cations (Li+, K+, Rb+, and Cs+) and NH4+ in Na−Ytype zeolite. The results showed that 32% of Na+ ions cannot be replaced by NH4+ at room temperature. Wang et al.9 compared the ion-exchange profile of NH4+ in three zeolites, namely, Na−clinoptilolite, Na−Y, and Na−P. They also calculated the thermodynamic parameters for the exchange processes. Langella et al.10 reported the exchange isotherm of © XXXX American Chemical Society
Received: April 4, 2013 Revised: June 7, 2013
A
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Figure 1. Structure of FAU-type zeolite and typical cation interaction sites. Blue balls denote typical cation sites.
Crystallographic structure data are obtained from the experimental neutron diffraction studies of Fitch et al.15 The FAU structure belongs to the Fd3m space group of symmetry and has a cubic length of 24.85 Å. Figure 1 shows the porous structure of a FAU-type zeolite. The structure consists of supercages and sodalite cages interconnected by hexagonal prisms. The typical sites for cations are shown in Figure 1. Site I is at the center of the hexagonal prism. Site I′ is located in the sodalite cage near the hexagonal prism. Site II is located in the supercage near the sodalite cage. Site III′ is located in the 12ring window of the supercage. In one unit cell, the number of site I is 16 and the number of sites I′, II, and III′ is 32. In the simulation model, 56 aluminum atoms are randomly distributed in the framework following the Löwenstein rule,16 which prevents adjacent Al−O−Al linkages, and Dempsey’s rule,17 which minimizes the number of adjacent Al−O−Si−O− Al linkages. Korányi et al.18 compared four different periodic unit models for FAU and found that only the body-diagonal model describes the distribution of Al in accordance with the experimental observation. The body-diagonal model, which fulfills the Löwenstein and Dempsey rules, is used in the current work. 2.2. DFT Calculations. The cluster model is used in the DFT calculations of the binding energies between cations and zeolite fragments. The 6-T model representing the zeolite is cut from the crystal structure, and the dangling bonds are saturated with hydrogen atoms. Four models are constructed to represent the different cation sites, as illustrated in Figure 2. Figure 2a carries a single six-membered ring (S6R) consisting of 6T (Si or Al) atoms with different numbers of aluminum atoms. One Al atom is in S6R (S6R-1Al; Figure 2a), and two Al atoms are in
molecular dynamics (MD) simulations. Assuming that the configurations are stable or metastable (during the simulation period) in phase spaces (given concentration, temperature, and pressure), free energy changes can be directly calculated and the results can be used to identify equilibrium points.12,13 Because the free energy changes associated with phase separations are usually very small, the configurations in the phase space are fairly stable for the MD simulations; this approach could be quite efficient, especially for predictions of dense phase separations such as liquid−liquid equilibriums. In addition, this approach is advantageous for its generality as no special modification of the simulation software is required. In this work, we conduct MD/TI simulations to calculate the free energy differences of Na+/NH4+ exchange in zeolite Y and predict the Na+/NH4+-exchange isotherms. The free energy is calculated from the configurations (distributions of different cations at different sites) generated using normal MC and MD simulations. Regarding the force field, we develop force field parameters based on quantum mechanics (QM) density functional theory (DFT) calculations, and then the force field is validated against QM calculation data and experimental data to determine cation distributions. The section 2 presents our calculation methods, section 3 contains the results and discussion, and section 4 summarizes our main findings.
2. MODELS AND METHODS 2.1. FAU Zeolite. FAU is an aluminosilicate zeolite with the general chemical formula MmAlmSi192−mO384, where M is the cation (Na+ or NH4+) and m is the number of cations in the unit cell, which also determines the Si/Al ratio. The model used in this work is Na56Y, which corresponds to Si/Al = 2.42. B
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Figure 2. Optimized zeolite 6T models: (a) S6R-1Al, (b) S6R-2Al-m, (c) S6R-2Al-p, and (d) D4R. (e) Side view of the S6R structure. (f) Side view of the D4R structure. The color codes are as follows: pink, aluminum; gray, silicon; red, oxygen; white, hydrogen; and purple, cation (Na+ or NH4+).
(LJ) 12−6 and Coulombic terms are used to represent the interactions.
meta-position (S6R-2Al-m; Figure 2b). Two Al atoms are in para-position (S6R-2Al-p; Figure 2c). Figure 2d is a 6T model with two adjacent four-membered rings (D4R), which represents the surrounding of site III′. Parts e and f of Figure 2 are the side views of the S6R and D4R structures, respectively. The clusters are optimized using the B3LYP/6-31+G(d,p) method with the T atoms and O atoms fixed. Using the optimized structures, the binding energies (BEs) are calculated by BE = Ecluster − Ecation − Eframe
⎡⎛ ⎞12 ⎛ ⎞6 ⎤ qiqj σij σij Uij = 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + r rij ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠
(2)
The LJ parameters are expressed as pairwise terms. The LJ term is used for the oxygen atoms for zeolite (Oz) only. The charge parameters (atomic partial charges) of the zeolite framework and pairwise LJ parameters between the Na+ and Oz were estimated from the work of Jaramillo et al.20 and optimized by fitting to the DFT binding energies and structures. The charge parameters of Oz atom type, and the LJ parameters, epsilon (εij) and sigma (σij) between the cations and Oz atom types, were optimized manually to fit the DFT energy data. The interaction parameters between zeolite Oz and water oxygen (Ow) were obtained from the work of Abrioux et al.21 The TIP3P model22 is used to describe the water molecules. The cations are treated as particles, each bearing a formal charge of 1. The LJ parameters between
(1)
where Ecluster denotes the optimized energy of the whole cluster, Ecation denotes the cation energy, and Eframe denotes the optimized energy of the zeolite framework. The DFT calculations are performed using the Gaussian package.19 2.3. Force Field. The interactions considered in this work include three parts, namely, zeolite, cation, and water. Considering that the rigid model is used, only Lennard-Jones C
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
From the equilibrated configurations, NVT MD simulations are conducted for the free energy calculations. The “multiconfiguration” TI formulation24 is used to calculate the free energy change (ΔG) when replacing Na+ with NH4+ in multiple steps.
cations and water were obtained from the work of Jensen et al.23 These parameters, together with the TIP3P parameters for water, accurately predict the relative hydration free energy of Na+ and NH4+. Table 1 lists the final force field parameters used in this work.
λ=1
ΔG =
Table 1. Atomic Partial Charges and Interaction Potential Parameters for FAU Zeolite, Na+, NH4+, and Water in This Worka charge (e)
Al Si Os Oa Na NH4 Ow Hw
1.75 1.65 −0.825 −1.1 1.0 1.0 −0.834 0.417 σ (Å)
ε (kJ/mol)
Na−Oz NH4−Oz Oz−Ow Na−Ow NH4−Ow Ow−Ow
3.376 3.797 3.077 3.581 4.101 3.151
0.038 56 0.038 56 0.709 0.036 47 0.036 47 0.636 4
(3)
P(λ) = (1 − λ)P Na+ + λP NH4+
(4)
where P represents the parameters. At each step, equilibrium runs are performed to calculate the ensemble-averaged partial derivative ⟨∂H/∂λ⟩, where H is the Hamiltonian. The quantities are numerically integrated to calculate the free energy difference (eq 3). In the MD simulations, the Nosé−Hoover method25,26 is used to control temperature, and the time step is 1 fs. The zeolite framework is fixed at the original position obtained from experimental crystallographic structure data. Water molecules are constrained using the SETTLE algorithm. In the MC and MD simulations, the LJ interactions are evaluated using a 12.0 Å cutoff and compensated with tail corrections. The electrostatic interactions are computed using the Ewald method for MC simulations and the particle mesh Ewald method27,28 for MD simulations. The cutoff distance of real space is 12.0 Å in both cases. The MC simulations are performed using Towhee7.0.2.29 The MD simulations are conducted using Gromacs4.5.5.30
Lennard-Jones Interaction atom pair
dλ λ
A coupling parameter (λ) is used to transform the system from the initial (λ = 0, Na+ ion) to the end (λ = 1, NH4+ ion) state. The charge parameters (+1) are kept constant; the LJ parameters and the masses of cations are coupled with linear interpolation:
Atomic Charges atom type
∂H ∂λ
∫λ=0
a
Oz, oxygen in the zeolite framework; Os, zeolite oxygen connected by two silicon atoms; Oa, zeolite oxygen connected by one silicon and one aluminum; Ow, oxygen atom in water.
2.4. Simulations. Standard MC and MD simulations are conducted. In these simulations, the zeolite framework and internal degree of freedom for water are fixed (the rigid model), and the positions of water molecules and cations are fully relaxed. The NVT (constant particle number N, volume V, and temperature T) MC simulations are performed to predict the equilibrium distribution of the cations in dehydrated and hydrated zeolite models. Each simulation contains 1 × 108 moves. Different moves are randomly selected with fixed probabilities. In dehydrated zeolite, the trail moves are translation (20%) and a superdisplacement move (80%) which take the molecular out of the box and place it back into the same box randomly. With water added, the trail moves are translation (20%), rotation (20%), and superdisplacement (60%). For systems containing both Na+ and NH4+, the trail moves include a particle position exchange move between Na+ and NH4+, which also provides fast equilibration for different types of cation distribution. The maximum translational and rotational displacements are adjusted to achieve an acceptance probability of 0.5. Using the MC generated configurations, NVT MD simulation annealing is conducted to further equilibrate the system and calculate cation distributions in dehydrated and hydrated zeolite models. The simulated annealing is carried out three times, and each has the following steps: heat to 1000 K in 50 ps, run for 500 ps at 1000 K, cool to 298 K in 50 ps, and then run for 500 ps at 298 K. Using the equilibrated configurations, the trajectory was collected every 1 ps to calculate the cation distributions. A cation is considered within a site if its position is S6R-1Al > D4R. This indicates that S6R structures (sites I′ D
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Table 3. Cation−Oxygen Distances r1, r2, and r3 (Å) and Distance of Cation from the 6T Plane d (Å) Calculated Using Quantum Mechanics and Force Field Methods (Figure 2) Quantum Mechanics Na+ S6R-1Al S6R-2Al-m S6R-2Al-p D4R
NH4+
r1
r2
r3
d
r1
r2
r3
d
2.19 2.20 2.21 2.29
2.31 2.23 2.29 2.29
2.45 2.42 2.35
0.95 0.78 0.78
2.61 2.63 2.64 2.67
2.87 2.75 2.81 2.70
3.19 2.94 2.82
2.02 1.83 1.81
Force Field Na S6R-1Al S6R-2Al-m S6R-2Al-p D4R
NH4+
+
r1
r2
r3
d
r1
r2
r3
d
2.17 2.19 2.21 2.26
2.35 2.19 2.31 2.27
2.76 2.52 2.31
1.03 0.74 0.74
2.47 2.48 2.48 2.58
2.70 2.52 2.56 2.63
3.09 2.77 2.62
1.70 1.47 1.41
and II) have stronger binding sites for Na+ or NH4+ than the D4R structures (site III′). These results are consistent with the DFT calculations conducted by Vayssilov et al.31 with smaller clusters whose Si and Al atoms are saturated with hydrogen atoms. The energy data also indicate that Na+ binds to the framework with a stronger force than NH4+. Table 3 compares the most important atom−atom distances calculated using force field and DFT methods. The optimized structures of the Na+ clusters according to the distances listed in Table 3 are illustrated in Figure 2. In the S6R structure, r1, r2, and r3 (Figure 2) are the distances between the cation and three oxygen atoms, and d is the distance between the cation and plane formed by the six T sites. The distances between the cations and oxygen are 2.2−2.4 Å and 2.6−2.9 Å for Na+ and NH4+, respectively. The d values range from 0.78 to 0.95 Å for Na+ and from 1.81 to 2.02 Å for NH4+. In the D4R structures, the cation is stabilized at the position above the Al atom, bisecting the O−Al−O angle. The variables r1 and r2, which represent the distances between the closest oxygen atoms, are listed in Table 3. The values are 2.29 Å for Na+ and 2.67−2.70 Å for NH4+. Overall, NH4+ is slightly farther from the framework and has weaker binding energies than Na+. The force field distances of all Na+ structures well agree with the QM results. However, the distances of NH4+ in the S6R structures are approximately 0.2 Å shorter than the DFT results, which is ascribed to the optimum balanced parameters for the binding energies and structures. Considering the simple spherical model used for NH4+, this deviation is still reasonable. The force field is validated by predicting the cation distributions in dehydrated zeolite, for which adequate experimental data can be used for comparison. The results are listed in Table 4. The calculated Na+ distributions are consistent with the experimental data measured by different groups.32−35 The finding that sites I′ and II are more populated by Na+ is consistent with the energy data showing that the S6R structure exhibits stronger binding energies than the D4R structure. The distributions of NH4+ in zeolite are also calculated and listed in Table 3. For NH4+, no experimental data are available for direct comparison. The experimental results of Lim et al.11 are those measured for a zeolite with a formula of (NH4)59H5Na5K2Y, which contains five sodium and two potassium ions. Nevertheless, our calculated data show the same trend. Previous experimental data have suggested that NH4+ ions are concentrated only at sites I′ and II, and our
Table 4. Calculated Cation Distributions in Dehydrated Zeolite Compared with the Experimental Dataa Na+ I I′ II III′ total a
NH4+
calcd
ref 32
ref 33
ref 34
ref 35
calcd
ref 11
5.9 13.5 31.0 5.6 56
7.04 13.76 29.44 3.76 54.00
8.00 18.88 30.08 0.04 57.00
4.00 17.60 32.00 1.4 55.00
9.30 13.70 25.30 3.50 51.80
10.0 8.9 29.0 8.1 56
0 30.7 28.8 0 59
Experimental data are labeled with references.
results show that sites I, I′, and II are more populated than site III′. The size of NH4+ is suitable for site I in the hexagonal prism. The cation interacts with six oxygen atoms in the zeolite framework. Thus, the results show that more NH4+ than Na+ ions are at site I. However, the experimental results show that all potassium ions remain at site I, which may have caused the difference between the experimental and simulation results. Table 5 lists the calculated cation−oxygen distances in dehydrated zeolite Y and experimental data for comparison. For Na−Y, the Na+−oxygen distance well agrees with the experimental data.15,33,36−38 For Na+ in site I, experimental data suggest that Na+ is at the center of the hexagonal prism. However, detailed analysis of the simulated data shows that Na+ is slightly away from the center of site I with a Na(I)−O(3) distance of 2.25−2.50 Å. This result is in accordance with other simulations20 and experiment studies.38 The distance between NH4+ and oxygen is slightly lower than the experiment result, which is consistent with the shorter distances predicted for the clusters discussed earlier. 3.2. Cation Distributions in Hydrated Zeolite Y. The distributions of cations in hydrated zeolite Y are predicted using five models representing different amounts of NH4+-exchanged z e o l i t e Y : ( N a 5 6 Y, N a 4 2 (NH 4 ) 1 4 Y, N a 2 8 ( NH 4 ) 2 8 Y , Na14(NH4)42Y, and (NH4)56Y). These compositions correspond to NH4+ mole fractions of 0, 0.25, 0.5, 0.75, and 1.0. For the hydrated model, the water amount was fixed at 250 per unit cell according to the experimentally deduced amount of 240− 260 for zeolite Y.39 In principle, the equilibrium water amount should be different between Na+ and NH4+ types of zeolites. Because the total amount is an approximated value, and the differences between the Na+ and NH4+ are small, it was assumed that the amount of water was unchanged. In fact, E
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
Table 5. Calculated Cation−Oxygen Distances (Å) in Dehydrated Zeolite Y Compared with Experimental Dataa Na(I)−O(3) Na(I′)−O(3) Na(II)−O(2)
calcd
ref 15
ref 33
ref 36
ref 37
ref 38
2.25−2.50 2.20−2.39 2.22−2.45
2.71 2.24 2.39
2.71 2.44 2.33 calcd
2.71 2.32 2.34
2.71 2.25 2.39
2.3−2.4
NH4(I′)−O(3) NH4(II)−O(2) a
2.50−2.78 2.52−2.69
ref 11 2.75 2.81
Experimental data are labeled with references.
The radial distribution functions (RDFs) between a single type of cation (NH4+ mole fraction = 1 or 0) and water oxygen atoms at different sites of hydrated zeolite and in aqueous solution are calculated. The results are shown in Figure 3 for
Jeffroy et al. predicted the water amounts for different alkali metal types of zeolites and found that the water amount was nearly constant.14 Table 6 shows the calculated distributions of Na+ and NH4+ in hydrated zeolite Y at different mole fractions of NH4+. At Table 6. Cation Distributions in Hydrated (250 Water Molecules) Zeolite Y with Different Mole Fractions of NH4+ x(NH4+) 0
0.25
0.5
0.75
1
I I′ II III′ total
0.3 17.0 23.4 15.3 56
0.0 12.0 20.5 9.5 42
0.0 10.0 12.3 5.7 28
0.0 6.3 5.7 2.0 14
0 0 0 0 0
I I′ II III′ total
0 0 0 0 0
4.0 2.0 5.3 2.7 14
6.0 4.0 12.7 5.3 28
6.0 6.7 20.6 8.7 42
7.0 11.7 29.2 8.1 56
Na+
NH4+
each composition, the numbers of NH4+ and Na+ are obtained. The results show that, at a low mole fraction of NH4+ (0.25), NH4+ prefers to be at sites I and II. With increased concentration of NH4+, the number of NH4+ at site I slightly increases, which suggests that site I is saturated by NH4+ at low mole fractions. At a high concentration of NH4+ (0.75), NH4+ ions are more concentrated at sites II and III′. At this concentration, the remaining Na+ ions are mostly concentrated at sites I′ and II. The influences of water on the distributions of pure Na+ and NH 4 + in hydrated zeolite are compared. The cation distributions in dehydrated zeolite (Table 4) show that Na+ ions at site I in the dehydrated cage move to site I′ in the hydrated cage because Na+ interacts with more water molecules at site I′ than that at site I. In addition, some Na+ ions at site II move to site III′ when water is present, which indicates the Na+ at site III′ is more favorable when water is present. As discussed earlier, the binding energy between the S6R-1Al structure (one possible configuration at site II) and Na+ is similar to the binding energy between the D4R structure (site III′) and Na+. However, Na+ ions at site III′ interact with more water molecules. Conversely, the strong binding energy between the S6R-2Al-p structure (another possible configuration at site II) and Na+ retains many Na+ ions at site II. For NH4+, the distributions are similar to that in dehydrated zeolite. Only a few NH4+ ions at site I move to site I′ for the same reason as Na+. The cation radius of NH4+ well fits the size of the hexagonal prism, so most NH4+ ions remain at site I.
Figure 3. Radial distribution functions between the cation and oxygen of water in aqueous solution and hydrated zeolite.
comparison. For all binding sites, the position of the first peak is consistent with the RDF of cation−Ow in aqueous solution (0.26 nm for Na−Ow and 0.30 nm for NH4−Ow). The cations and water have similar contact radii in zeolite and aqueous solution. However, the peak heights differ because of space limitation. The coordination numbers of water near the cations at different sites are calculated by integrating the first peak of the F
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
chemical potentials in zeolite and solution are equal at given compositions. Using x to represent the mole fraction of NH4+, the Gibbs free energy of eq 5 can be described as
corresponding RDF. The results are listed in Table 7. For cations in aqueous solution, the coordination number of Na+ is Table 7. Calculated Coordination Numbers of Water Near Cations in Hydrated Zeolite and Aqueous Solution aq soln I′ II III′
Na+
NH4+
5.9 2.3 2.4 4.0
7.2 2.5 4.5 5.4
zeo sol + ΔΔG(xzeo , xsol) = ΔG Na → NH4 +(xzeo) − ΔG Na +→ NH4 +(xsol)
(6)
Equation 6 shows the relationship between the mole fractions in zeolite and solution. To determine the free energy differences, we calculate the free energy changes when one Na+ is replaced with NH4+ in zeolite and solution. The free energy change calculated by changing one Na+ to NH4+ is equivalent, by finite difference, to the change in chemical potential. In aqueous solution, the free energy difference of changing one Na+ to NH4+ can be expressed by
5.9, whereas that of NH4+ is 7.2 because of the larger ionic radii of NH4+. In zeolite, the coordination numbers vary at different sites. In the sodalite cage at site I′, both Na+ and NH4+ coordinate with ∼2.4 water molecules. In site II, Na + coordinates with 2.4 water molecules, whereas NH 4 + coordinates with 4.5 water molecules. This result is due to the greater distance of NH4+ from the ring and its larger contact area for water in the supercage. At site III′, the water coordination numbers for both Na+ and NH4+ are slightly smaller than those in aqueous solution. 3.3. Free Energy Calculations. The simulation protocol for TI calculations is established and then validated by calculating the hydration free energy difference resulting from the replacement of Na+ by NH4+ in aqueous solutions and hydrated zeolites. The TI calculations are conducted using 26 decoupling λ points that are equally distributed from 0 to 1. The simulation time is 1 ns each for equilibration and data collection at each λ point. The typical TI energy profile in Figure 4 shows a fairly smooth curve suitable for numerical
sol + ΔG Na → NH4 + = v Na +μ Na + + v NH4 +μ NH + = μ NH + − μ Na + 4
4
(7)
where v is the stoichiometric coefficient. vNH4+ = 1, vNa+ = −1, and μ is the chemical potential which can be expressed as a simple relation: 0 μ Na+ = μ Na (1 − xsol) + + RT ln γ Na +
(8)
0 μ NH + = μ NH x + + RT ln γ NH + sol 4
4
4
Combining eqs 7 and 8, we apply the following equation to calculate the free energy differences at different mole fractions of NH4+. γNH + 4 sol 0,sol + ΔG Na + RT → NH4 +(xsol) = ΔG Na +→ NH4 + + RT ln γNa+ xsol ln 1 − xsol (9) where γ is the activity coefficient of the cation. The reference hydration free energy ΔG0 must be calculated. The reference concentration we use is 0.1 M, which corresponds to a ratio of 1 cation to approximately 500 water molecules. As the concentration is close to that of the reference state, we assume that the activity coefficients are the same for Na+ and NH4+. Thus, the free energy change is simplified as xsol sol 0,sol + ΔG Na → NH4 +(xsol) = ΔG Na +→ NH4 + + RT ln 1 − xsol (10) +
The free energy change when one Na is replaced with NH4+ is calculated using eq 10 and shown in Figure 5. The free energy differences when one Na+ is replaced with NH4+ in solution range from ca. 70 to 95 kJ/mol depending on the mole fraction of NH4+. In zeolite, multiple simulations with different mole fractions of NH4+ are implemented to calculate free energy differences. For the free energy differences of replacing one Na+ with NH4+ in zeolite, we explicitly calculate the free energy changes for each site at a given mole fraction of NH4+ in the simulation boxes. In each simulation, only one Na+ is replaced with NH4+. For each binding site, five cations are randomly selected to calculate the average free energy change. If the number of Na+ ions is less than five in a binding site, all Na+ ions at this site are included in the calculation. These calculations were repeated using three independent initial configurations. The results are used to calculate the average values and errors based on the assumption that the configurations for the TI calculations are
Figure 4. Calculated free energy profiles for the exchange of Na+ with NH4+ in aqueous solution and hydrated zeolite.
integration. The calculated hydration free energy difference for the replacement of Na+ with NH4+ under constant conditions (0.1 M, 298 K, 1 atm) is 81.4 kJ/mol, which is consistent with Jensen’s result23 of 80.3 kJ/mol and experimental data of 79.8 kJ/mol.40 The Na+ and NH4+ ion-exchange equilibrium between zeolite and solution can be described by the following relation: Na +zeo + NH4 +sol ⇔ Na sol + + NH4 +zeo
(5)
where “zeo” and “sol” mean hydrated zeolite and aqueous solution, respectively. Equilibrium is achieved when the G
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
where Pi is the probability of the exchange at site i, and Pi is calculated by Pi =
Ni e−ΔGi / RT ∑i Ni e−ΔGi / RT
(12) +
where Ni is the number of Na ions at site i. The statistically averaged results are listed in the last column of Table 8. The average free energy increases with increased NH4+ mole fraction. The range of changes from 74.6 (x = 0.0) to 84.2 (x = 0.75) is small (9.6 kJ/mol). At the equilibrium point between solution and zeolite, ΔΔG in eq 6 is equal to 0. Thus, using the free energy difference calculated in zeolite at different mole fractions, we can calculate the corresponding mole fraction of NH4+ in aqueous solution as follows: xsol zeo 0,sol + ΔG Na → NH4 +(xzeo) − ΔG Na +→ NH4 + = RT ln 1 − xsol
Figure 5. Free energy difference of exchanging Na+ with NH4+ in aqueous solution at different NH4+ mole fractions.
(13)
NH4+
The corresponding mole fractions in aqueous solution are calculated using eq 13, and the results are shown in Figure 6. The diagonal line indicates the equal partition of zeolite and
stable. This assumption is fairly solid because when the coupling parameter is changed from 0 to 1 (Na+ to NH4+), only LJ parameters are changed and the charge parameter (+1) is constant, which helps stabilize the configurations. The ΔG values for different sites at different mole fractions of NH4+ in hydrated zeolite are listed in Table 8. At a low mole Table 8. Calculated Free Energy Difference ΔG (kJ/mol) of Exchanging Na+ with NH4+ for Each Site at Different NH4+ Mole Fractions in Hydrated Zeolitea xzeo
I′
0 0.25 0.5 0.75
105.6 (4.2) 89.5 (0.9) 95.4 (1.6) 102.8 (2.4)
II 74.3 76.2 78.5 84.4
(1.4) (1.1) (0.8) (1.3)
III′ 81.2 81.2 82.0 84.0
(0.2) (0.9) (0.5) (0.9)
av 74.6 76.6 78.9 84.2
a
Uncertainties (error bars) are enclosed in parentheses. Statistically averaged ΔG values are listed in the last column.
fraction of NH4+ (xzeo = 0), the Na+ ions at site II are most favorably exchanged. With increased xzeo, the exchange free energy at site I′ decreases from 105.6 to 89.5 kJ/mol. This result is due to the movement of some NH4+ ions to site I and the decreased number of cations in the sodalite cage from 17 to 14, which reduce the repulsion among the cations in the sodalite cage. However, the exchange free energy at site I′ increases to 95.4 kJ/mol at xzeo = 0.5 because more NH4+ ions occupy site I′, and the repulsion among the cations in the sodalite cage is greater because NH4+ has a larger radius than Na+. The exchange free energy at sites II and III′ increases with increased xzeo, which indicates that the exchange process in the supercage becomes more difficult as more NH4+ is exchanged. The results suggest that site II is the most favorable site for NH4+ exchange throughout the entire exchange process from the thermodynamic point of view. The exchange process at site I′ is difficult as indicated by the free energy of ∼90−105 kJ/mol throughout the entire process. To predict the overall Na+/NH4+-exchange isotherm in zeolite Y, the averaged ΔG for zeolite is calculated as follows: ΔGav =
∑ PiΔGi site
Figure 6. Isotherm of Na+/NH4+ exchange in zeolite Y. The gray dashed line is the diagonal for a clear view. The dots are the calculated values in this work, and the solid lines are based on the experimental data of Sherry,7 Theng et al.,8 and Wang et al.9.
solution. The data above the diagonal line indicate that the mole fraction of NH4+ in zeolite is higher than that in solution; thus, the exchange is effective. Overall, the simulated isotherm shows the same trend as that of the experimental data.7−9 At low NH4+ mole fractions, the calculated data well agree with the three measured data, indicating that the exchange is thermodynamically effective. With increased mole fraction of NH4+, maximum effectiveness is achieved at approximately 0.2−0.3 mole fraction of the solution, and then the effectiveness decreases. The calculated maximum effectiveness is very close to all three measurements, and the actual data agree with those of Wang et al., but are higher than those reported in two earlier works.7−9 The exchange is effective until the mole fraction reaches ca. 0.75 instead of 0.6 based on our data and experimental data.7−9
(11) H
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
(5) Sato, K.; Nishimura, Y.; Honna, K.; Matsubayashi, N.; Shimada, H. Role of HY Zeolite Mesopores in Hydrocracking of Heavy Oils. J. Catal. 2001, 200, 288−297. (6) Frising, T.; Leflaive, P. Extraframework Cation Distributions in X and Y Faujasite Zeolites: A Review. Microporous Mesoporous Mater. 2008, 114, 27−63. (7) Sherry, H. S. The Ion-Exchange Properties of Zeolites. I. Univalent Ion Exchange in Synthetic Faujasite. J. Phys. Chem. 1966, 70, 1158−1168. (8) Theng, B. K. G.; Vansant, E.; Uytterhoeven, J. B. Ion Exchange in Synthetic Zeolites. Part 1. Ammonium and Some of Its Alkyl Derivatives in Linde Sieves X and Y. Trans. Faraday Soc. 1968, 64, 3370−3382. (9) Wang, Y.; Lin, F.; Pang, W. Ion Exchange of Ammonium in Natural and Synthesized Zeolites. J. Hazard. Mater. 2008, 160, 371− 375. (10) Langella, A.; Pansini, M.; Cappelletti, P.; de Gennaro, B.; de’Gennaro, M.; Colella, C. NH4+, Cu2+, Zn2+, Cd2+ and Pb2+ Exchange for Na+ in a Sedimentary Clinoptilolite, North Sardinia, Italy. Microporous Mesoporous Mater. 2000, 37, 337−343. (11) Lim, W. T.; Seo, S. M.; Wang, L.; Lu, G. Q.; Heo, N. H.; Seff, K. Single-Crystal Structures of Highly-Exchanged, Fully Deaminated, and Fully Tl+-Exchanged Zeolite Y (FAU, Si/Al = 1.56), All Fully Dehydrated. Microporous Mesoporous Mater. 2010, 129, 11−21. (12) Liu, Y.; Li, X.; Wang, L.; Sun, H. Prediction of Partition Coefficients and Infinite Dilution Activity Coefficients of 1-ethylpropylamine and 3-methyl-1-pentanol Using Force Field Methods. Fluid Phase Equilib. 2009, 285, 19−23. (13) Cheng, T.; Li, F.; Dai, J.; Sun, H. Prediction of the Mutual Solubility of Water and Dipropylene Glycol Dimethyl Ether Using Molecular Dynamics Simulation. Fluid Phase Equilib. 2012, 314, 1−6. (14) Jeffroy, M.; Boutin, A.; Fuchs, A. H. Understanding the Equilibrium Ion Exchange Properties in Faujasite Zeolite from Monte Carlo Simulations. J. Phys. Chem. B 2011, 115, 15059−15066. (15) Fitch, A. N.; Jobic, H.; Renouprez, A. Localization of Benzene in Sodium-Y-Zeolite by Powder Neutron Diffraction. J. Phys. Chem. 1986, 90, 1311−1318. (16) Löwenstein, W. The Distribution of Aluminum in The Tetrahedra of Silicates and Aluminates. Am. Mineral. 1954, 39, 92−96. (17) Dempsey, E. Calculation of Madelung Potentials for FaujasiteType Zeolites. I. J. Phys. Chem. 1969, 73, 3660−3668. (18) Korányi, T. I. Distribution of Aluminum in the Periodical Building Units of Faujasites. J. Phys. Chem. C 2007, 111, 2520−2524. (19) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; et al. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (20) Jaramillo, E.; Auerbach, S. M. New Force Field for Na Cations in Faujasite-Type Zeolites. J. Phys. Chem. B 1999, 103, 9589−9594. (21) 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. (22) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (23) Jensen, K. P.; Jorgensen, W. L. Halide, Ammonium, and Alkali Metal Ion Parameters for Modeling Aqueous Solutions. J. Chem. Theory Comput. 2006, 2, 1499−1509. (24) Kollman, P. Free Energy Calculations: Applications to Chemical and Biochemical Phenomena. Chem. Rev. 1993, 93, 2395−2417. (25) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (26) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-Space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (27) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N*log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092.
4. CONCLUSION Using MD and TI methods, we calculate the free energy differences of replacing Na+ with NH4+ in zeolite Y and aqueous solution. The predicted isotherm of the Na+/NH4+ exchange in zeolite Y well agrees with the experimental data. The successful prediction is ascribed to the quality of the underlying force field. We optimize the force field parameters obtained from the literature by comparing them to QM DFT data. The force field accurately describes the differences in bonding energies between the two cations within water and zeolite frameworks. The agreement with the experimental data on the distribution of Na+ and NH4+ in dehydrated zeolites validates the force field. Another critical step in this work is the calculation of free energy differences for exchanging Na+ with NH4+ in zeolite Y with high precision. The calculation is conducted based on extensive simulations of the distributions of different cations in hydrated zeolites with different mole fractions of NH4+. Although no experimental data are available for direct comparison, the distributions of cations are reasonable as shown by the analysis of the binding energies of cations with water and zeolite, as well as the geographic features of different sites of the zeolite. With the obtained distributions, the free energy changes at different sites with different mole fractions of NH4+ are calculated. Using the calculated free energy changes at different sites, the statistically weighted average values are determined and compared with the free energy differences of exchange in aqueous solutions to predict the equilibrium isotherm. The free energy differences of replacing Na+ with NH4+ in solution are directly calculated. The free energy difference of replacing Na+ by NH4+ at the reference state well agrees with the experimental data, which further validates the TI calculation protocol and force field. Overall, our findings indicate that the MD/TI method can be used to predict ion-exchange isotherms using an accurate force field. The configurations of cations in zeolite can be fully explored, and the configurations are stable with respect to the disturbance of exchange. The last condition may not always be true but may be a valid assumption, at least for cations with similar sizes and charges.
■
AUTHOR INFORMATION
Corresponding Author
*Tel.: +86-21-5474-8987. Fax: +86-21-5474-1297. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported in part by the National Science Foundation of China (Nos. 21073119 and 21173146).
■
REFERENCES
(1) van Santen, R. A.; Kramer, G. J. Reactivity Theory of Zeolitic Brønsted Acidic Sites. Chem. Rev. 1995, 95, 637−660. (2) Barthomeuf, D. Basic Zeolites: Characterization and Uses in Adsorption and Catalysis. Catal. Rev. 1996, 38, 521−612. (3) Zecchina, A.; Arean, C. O. Diatomic Molecular Probes for Mid-IR Studies of Zeolites. Chem. Soc. Rev. 1996, 25, 187−197. (4) Corma, A. From Microporous to Mesoporous Molecular Sieve Materials and Their Use in Catalysis. Chem. Rev. 1997, 97, 2373− 2420. I
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX
The Journal of Physical Chemistry C
Article
(28) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (29) Martin, M. G. MCCCS Towhee. http://towhee.sourceforge.net (accessed 2006). (30) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701−1718. (31) Vayssilov, G. N.; Staufer, M.; Belling, T.; Neyman, K. M.; Knözinger, H.; Rösch, N. Density Functional Studies of AlkaliExchanged Zeolites. Cation Location at Six-Rings of Different Aluminum Content. J. Phys. Chem. B 1999, 103, 7920−7928. (32) Van Dun, J. J.; Dhaeze, K.; Mortier, W. J. TemperatureDependent Cation Distribution in Zeolites. 2. Dehydrated NaxHY (x = 13, 23, 42, 54), Ca15HY, and Sr27Y. J. Phys. Chem. 1988, 92, 6747− 6754. (33) Eulenberger, G. R.; Shoemaker, D. P.; Keil, J. G. Crystal Structures of Hydrated and Dehydrated Synthetic Zeolites with Faujasite Aluminosilicate Frameworks. I. The Dehydrated Sodium, Potassium, and Silver Forms. J. Phys. Chem. 1967, 71, 1812−1819. (34) Jirák, Z.; Vratislav, S.; Bosácě k, V. A Neutron Diffraction Study of H, Na-Y Zeolites. J. Phys. Chem. Solids 1980, 41, 1089−1095. (35) Marra, G. L.; Fitch, A. N.; Zecchina, A.; Ricchiardi, G.; Salvalaggio, M.; Bordiga, S.; Lamberti, C. Cation Location in Dehydrated Na−Rb−Y Zeolite: An XRD and IR Study. J. Phys. Chem. B 1997, 101, 10653−10660. (36) Grey, C. P.; Poshni, F. I.; Gualtieri, A. F.; Norby, P.; Hanson, J. C.; Corbin, D. R. Combined MAS NMR and X-ray Powder Diffraction Structural Characterization of Hydrofluorocarbon-134 Adsorbed on Zeolite NaY: Observation of Cation Migration and Strong Sorbate− Cation Interactions. J. Am. Chem. Soc. 1997, 119, 1981−1989. (37) Lievens, J. L.; Mortier, W. J.; Verduijn, J. P. Influence of The Framework Composition on The Cation-Site Energy: The Structures of Dehydrated NaxHGaY (x = 54, 36, and 21) Zeolites. J. Phys. Chem. 1992, 96, 5473−5477. (38) Engelhardt, G. Cation Location in Dehydrated Zeolite NaY Revisited: SI Position is Displaced from the Center of the Hexagonal Prism. Microporous Mater. 1997, 12, 369−373. (39) Breck, D. W. Zeolite Molecular Sieves; Wiley: New York, 1974. (40) Marcus, Y. A Simple Empirical Model Describing the Thermodynamics of Hydration of Ions of Widely Varying Charges, Sizes, and Shapes. Biophys. Chem. 1994, 51, 111−127.
J
dx.doi.org/10.1021/jp403326n | J. Phys. Chem. C XXXX, XXX, XXX−XXX