J. Phys. Chem. C 2009, 113, 12225–12235
12225
Potential Model Development Using Quantum Chemical Information for Molecular Simulation of Adsorption Equilibria of Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4 E´va Csa´nyi,† Tama´s Kristo´f,*,† and Gyo¨rgy Lendvay*,‡,§ Department of Physical Chemistry, UniVersity of Pannonia, H-8201 Veszpre´m, P.O. Box 158, Hungary, Department of General and Inorganic Chemistry, UniVersity of Pannonia, H-8201 Veszpre´m, P.O. Box 158, Hungary, and Chemical Research Center, Hungarian Academy of Sciences, P.O. Box 17, H-1525 Budapest, Hungary ReceiVed: March 20, 2009; ReVised Manuscript ReceiVed: May 11, 2009
A theoretical examination of two simple Lennard-Jones + Coulomb type potential models of zeolite NaA-4 was performed by quantum chemical calculations using density functional theory. A comparison of potential curves obtained from model parameters with those from quantum chemical calculations as a function of the position of a Na+ ion with respect to an appropriately chosen fragment of zeolite was used as a primary indicator in testing how realistic the behavior of the applied models is. This comparison provides a good starting point for further parameter tuning of a potential model using experimental data. The potential models were tested in different aspects: structural and adsorption equilibrium properties (quantum chemical structures vs pair correlation functions as well as equilibrium amount of adsorption) were calculated with pure water, methanol, and ethanol. Tuning the potential parameters resulted in a new model that satisfactorily reproduces the quantum chemical equilibrium location of the three adsorbates with respect to the zeolite framework. From the pair correlation functions, the behavior of the new potential model was found to be between those of the two investigated earlier models. Monte Carlo simulation of adsorption for water-methanol and water-ethanol mixtures showed that, at higher pressures (above p ) 1 kPa), the newer model predicts greater adsorption selectivity of water to alcohols than the earlier models, which fits to available membrane permeation experiments better. 1. Introduction Zeolites are porous materials consisting of a silicon-aluminum-oxygen network with the positive charge of the Al ions being supplemented by mobile cations. They have excellent adsorption and molecular sieve properties, making them good catalysts and ion-exchange agents. They are stable in harsh chemical environments and at high temperatures. The size of the pores in different zeolites covers a wide range, and accordingly, the type of molecules that can penetrate them and be sorbed varies widely for different types of zeolites. A synthetic variant, zeolite NaA-4, is characterized by small pore diameters (∼0.8 nm). This type of zeolite can be used for dehydration of organic solvents. An important application of zeolite NaA-4 is the dehydration of ethyl alcohol, being used industrially in bioethanol production. The performance of zeolite NaA-4 as a drying agent can, in principle, be determined experimentally. However, due to the complexity of such experiments, only a limited set of experimental data is available. A possible way to collect the missing information is a theoretical approach, such as molecular simulation. Using properly parametrized molecular simulations, one can determine properties, such as the equilibrium amount of adsorption, heat of adsorption, and kinetic parameters, for example, diffusion coefficients. * Corresponding author. Email:
[email protected] (T.K.); lendvay@ chemres.hu (G.L.). † Department of Physical Chemistry, University of Pannonia. ‡ Department of General and Inorganic Chemistry, University of Pannonia. § Hungarian Academy of Sciences.
The most important criterion for an atomistic-level molecular simulation to produce successful and realistic information on adsorption is that an appropriate potential function characterizing the interaction of the atoms within the adsorbent as well as between atoms of the adsorbent and adsorbate be available. For condensed-phase systems, in principle, one can use the methods of electronic structure theory to derive the potential energy of an atomistic system, but the efficiency of the computational implementations does not yet allow productive simulation of a zeolite-adsorbent system. Instead, model potentials are used that are expected to properly describe the interactions between adsorbate and adsorbent and adequately reflect the polar properties of the adsorbent and the structure of the adsorbent as well as that of the adsorbent/adsorbate system. The model potentials are most often constructed from atom-atom pair potentials depending only on the distance between the respective atoms. For these pair potentials, appropriate functional forms are chosen and their parameters are determined empirically-with proper guidelines-to ensure that the experimental data on the properties to be modeled are reproduced satisfactorily well. The pair potentials often differ significantly from those measured or calculated ab initio for a single pair of particles (as it is the case with some of the most widely used potentials, such as SPC/E for water1). This is not surprising because (1) they are expected to characterize a condensed-phase system where the location of the numerous neighboring particles influences the apparent two-body interactions so that the latter can be described by pair potentials only if proper effectiVe parameters are selected and (2) our target is the description of one or more macroscopic
10.1021/jp902520p CCC: $40.75 2009 American Chemical Society Published on Web 06/17/2009
12226
J. Phys. Chem. C, Vol. 113, No. 28, 2009
observables that, reflecting the average behavior of a large number of particles, only indirectly depend on the interaction of individual atom pairs. According to the generally applied approach, the pair potentials are constructed from terms corresponding to dispersive and Coulomb interactions. The dispersion interaction is usually described by the Lennard-Jones (LJ) potential. To model the Coulomb interactions, the charge of the atoms/ions may differ from the integer values normally considered in chemistry, again to produce an effective, average interaction. Model potentials for simulation of adsorption of molecules in zeolites have to be able to describe the interaction of the atoms of the zeolite framework with each other as well as with those of the adsorbate. There are several models for the zeolite framework that can be classified as flexible,2 rigid,3-5 and semiflexible.6 The rigid model implies that the atomic positions are fixed in the calculations. In a semiflexible model, the locations of the framework atoms are fixed but the sodium ions can move. The flexible model allows the framework atoms also to move. Furthermore, there are oxygen atoms in various positions in the framework that can be characterized by different parameter sets. Recently we have studied the adsorption characteristics of zeolite NaA-4 by molecular simulations.7,8 We used several force fields for the zeolite framework taken from the literature, but none of them proved to be completely satisfactory. To be able to provide results on the selectivity of adsorption of water with respect to alcohols, we found that it is necessary to develop a new zeolite framework model. We decided to develop such a potential model by modifying an earlier force field by varying the parameters using a new paradigm. According to this paradigm, we sequentially assess the performance of the actual model potential in two respects: (1) how do the potential curves calculated with the model along some selected lines in the configuration space of the system match those obtained using electronic structure theory and (2) how well does the molecular simulation reproduce macroscopic properties, such as the amount of adsorption measured experimentally. We perform the molecular simulation only if the first criterion is satisfied. We expect that the comparison with quantum chemical potential data ensures that the potential parameters are realistic, and so they provide a good starting point for parameter tuning using experimental data. This way, we can get a reasonably accurate model that enables us to make predictions. In principle, we could fit the parameters of the model so that the quantum chemical potential curves are reproduced as well as possible. However, our own experience9,10 and those of others11,12 with potential surface fitting show that the number of parameters in a simplified model potential consisting of pairwise interaction functions is too little to guarantee that a complex potential surface, with all its minima, steepness of walls, etc., could be uniformly well reproduced. An additional reason making the attempt of a full potential surface fit less desirable is that quantum chemical calculations can be performed for only a subset of the complete system, and if we fitted the parameters of the model potential to the quantum chemical potential surface for the subset, we could easily face the contradiction mentioned above, namely, that the accurate potential that is valid for a limited number of particles separated from the rest of the system does not match the model potential that best reproduces the macroscopic data. An advantage of our paradigm is that it makes the search for good parameters of the pair potentials more efficient. By checking whether a parameter combination proposed during
Csa´nyi et al. parameter search is realistic, we can discard the ones that will probably produce unacceptable results in a complete molecular dynamics or Monte Carlo (MC) simulation, this way reducing the number of attempts that actually involve the rather timeconsuming step, simulation. In our search for good potential parameters for zeolite NaA4, we supplemented the parameter improvement procedure we reported earlier (and sketch below) with additional steps according to the strategy described above. The first model we used is the zeolite model proposed by Lee et al.6 (in the following, model A). This model predicted too high equilibrium loadings in our simulations for pure methanol and ethanol with respect to the corresponding experiments.13 To correct this weakness of the model, we developed a new model (model B7) in which we kept the semiflexible concept of the framework. The atomic positions of the zeolite NaA-4 framework were taken from X-ray diffraction experiments, and the point charges as well as the Lennard-Jones interaction sites were assigned to point masses corresponding to the framework atoms. For simplicity, all types of oxygen atoms were characterized by the same parameter set and the dispersion (as well as the soft-sphere repulsion) interactions of framework atoms were taken into account only through the effective LJ potential of O. Model B shows better agreement with the experimental data for the loadings of water, methanol, and ethanol. Furthermore, in simulations of mixtures, the selectivity of water to alcohols is predicted to be higher and, in contrast to model A, model B does not exhibit any inversion of selectivity in a wide range of pressure.7 The main advantage of model B is its transferability: it can also be applied to zeolites with a variable Si/Al ratio. However, in simulations of the adsorption of CO2 in zeolite NaA-4, the experimental equilibrium loadings were significantly underestimated.14 As CO2 is a nonpolar gas under normal conditions, we conjectured from this result that the model is too polar. The body of information listed above served as a starting point for improvement according to the above-mentioned paradigm. We calculated potential curves, that is, cuts of the potential energy surface, using the methods of electronic structure theory. These calculations, in addition to providing information on the interaction of the Na+ ion and the zeolite framework, also supported our assumption that the charge parameter of the Na+ ion needs to be reduced. As described above, comparison of the model potential curves with the quantum chemical ones provided a quick way of selecting the parameter domain that produces realistic interaction potentials, and the time-consuming molecular simulations were performed only within this domain. Using only two-body (Lennard-Jones and Coulomb) potentials, we cannot expect that fitting to the quantum chemical potential surface would make sense. Instead, the target of parameter optimization-within the domain approved by the quantum chemical calculations-was the best performance of the model potential in classical MC simulations of adsorption of pure water, methanol, and ethanol on zeolite NaA-4. In this paper, we first describe the methods and procedures used in the electronic structure calculations, followed by the summary of the methods of molecular simulations. The basic principles and the parameters of the group of models we developed are presented in section 3. In section 4, we compare the performance of models B and C in simulation of adsorption of pure adsorbates in zeolite NaA-4 and test how realistic the atom-atom pair potentials are by comparing them with the DFT potential curves. Finally, the results of modeling simultaneous
Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4
Figure 1. Zeolite fragment (a corner of the supercage) used in QC calculations. The Si, Al, and O atoms are denoted by yellow, gray, and red, respectively. The reference coordinate system (axes x, y, z) and the different types of O atoms are also shown.
adsorption, the main goal of the development of the new model, are presented. 2. Computational Details Density Functional Theory Calculations. We used firstprinciples electronic structure methods, namely, density functional theory (DFT)15-18 with the B3-LYP combination of functionals19,20 and the SVP basis set.21 Since these calculations are relatively time-consuming, they were performed using a fragment of the zeolite that can still adequately represent the behavior of the zeolite adsorbent. Adsorption takes place in the interior of the supercage of the zeolite framework. After some preliminary calculations, we have chosen a corner of the supercage, which consists of one 8R, two 6R, and two 4R windows connected to each other, as seen in Figure 1 (the reference coordinate system used is also shown). The terminal O atoms were capped by H atoms. The geometry of the fragment was frozen during the calculations because the zeolite framework was kept fixed during the classical simulations, too. We have checked and found that minor modifications of the geometry of the fragment do not cause any significant difference in the results of the DFT potential curve calculations. All quantum chemical calculations were carried out with the Gaussian 03 suite of programs.22 The geometry corresponding to the minima on the potential energy surface of the Na+ ion + zeolite system was calculated by restricted geometry optimization, keeping the zeolite fragment frozen. The optimum positions were determined for one as well as for two Na+ ions. The obtained positions agree with the results of X-ray diffraction experiments.23 Note that at least one of the Na+ ions is very mobile. One Na+ ion is located close to the center of the 6R ring. If a second Na+ ion is also present, the potential energy is the lowest if the first ion is at the center of the 6R ring and the second is near the center of the 8R ring. Later, when we optimized the location of the adsorbates, the sodium atoms were initially placed at these positions but were allowed to move. We also calculated DFT potential energy curves by systematically varying the position of the Na+ ions, which served as a reference for testing the pair potential models.
J. Phys. Chem. C, Vol. 113, No. 28, 2009 12227 The Na+ ion was translated along various straight lines anchored to the zeolite fragment, which are expected to provide satisfactory information for fitting of the parameters of the Na+-zeolite potential. The adsorption of water, methanol, and ethanol molecules in the supercage of zeolite NaA-4 was simulated in DFT calculations including the model zeolite fragment, two sodium ions, as well as one adsorbate molecule. The potential minima corresponding to the lowest-energy position of the adsorbate with respect to the Na+ + zeolite fragment were located by restricted geometry optimization. The framework of the zeolite fragment was frozen; the Na+ ions and the adsorbate molecule were allowed to move. The initial positions of Na+ ions were at the positions obtained in previous optimization in the absence of the adsorbate and remained essentially there during the optimization. The structure of the adsorbate molecule in the DFT calculations was set initially identical to the structure used in the classical molecular simulations (as it is shown in Table 2), but after the optimum location of the adsorbate molecule was found, this restriction has been relaxed and a full optimization of the adsorbate with respect to the (frozen) zeolite framework was performed. Monte Carlo Simulations. The molecular simulation of adsorption was carried out using the grand canonical MC method. The partial pressure values of the adsorbate molecules in the gas phase were given indirectly by specifying the component’s chemical potential that was calculated from the ideal gas law. In our earlier work, the suitability of the ideal gas law was verified by test simulations.7 The LJ interactions between sites of the zeolite and the adsorbate molecules, such as water, methanol, and ethanol, were calculated by the Lorentz-Berthelot combining rules. The Wolf method24 was used to treat the Coulomb interactions. The advantage of this method is that it lacks the inherent periodicity of Ewald summation,25 and consequently, the computation time can be reduced. It was shown earlier26 that the Wolf method is an efficient alternative to the Ewald summation, and in most systems, the best agreement can be achieved with rc ) L/2 and R ) 2/rc (where L is the box length, rc is the cutoff radius, and R is the convergence parameter). In our system, these parameters are the following: L ) 2.4555 nm, rc ) 1.22775 nm, and R ) 1.6290 1/nm. The grand canonical simulations with pure components consisted of an equilibration period (in the order of 108 MC moves) and a subsequent averaging period (at least 2 × 108 MC moves), where the ratio of insertion and deletion steps was 70-80%. In these simulations, we used the configurationalbias technique27 to increase the sampling efficiency. In mixture simulations, the sampling efficiency was further increased by identity change attempts28 and, in the case of water-ethanol mixtures, a special, so-called energy biased identity change method was used.29 The creation of molecules inside the sodalite cages was prevented artificially by placing purely repulsive dummy atoms at the center of these cages because the standard random insertion of molecules cannot take into account the physical diffusion pathways in the zeolite. The accessibility of the sodalite cages for adsorbate molecules is still an open question; there is some disagreement in the literature concerning it.30,31 Sodalite cages are very probably not accessible for molecules such as methanol and ethanol. However, in contrast to some experimental reports (e.g., ref 32), our previous simulation results suggest that smaller molecules such as water are also unable to pass through the windows of the sodalite cage.
12228
J. Phys. Chem. C, Vol. 113, No. 28, 2009
Csa´nyi et al.
TABLE 1: Lennard-Jones Energy (ε) and Size (σ) Parameters and Partial Charges (q) of Models A,6 B,7 and C model A sites
Na+
Si
Al
Oa
Ob
Oc
(ε/k)/K σ/nm q/(electron charge)
2507.3 0.1776 0.5502
64.18 0.4009 0.6081
64.18 0.4009 0.6081
78.02 0.289 -0.4431
78.02 0.289 -0.4473
78.02 0.289 -0.438
TABLE 2: Lennard-Jones Energy (ε) and Size (σ) Parameters, Partial Charges (q), and Geometry Data (Bond Length δ, Bond Angle r and Dihedral Angle θ) for Water, Methanol, and Ethanol model
sites
(ε/k)/K
σ/nm
q/(electron charge)
water34 O 78.197 0.3166 0.8476 H 0.4238 δ(O-H) ) 0.1 nm, R(H-O-H) ) 109.47°
model B sites
Na+
(ε/k)/K σ/nm q/(electron charge)
20.0 0.32 1.0
Si
Al
O
3.7
2.7
100.0 0.34 -1.85
Si
Al
O
1.7
200.0 0.33 -1.2
model C sites
Na+
(ε/k)/K σ/nm q/(electron charge)
100.0 0.25 0.7
a
Member of 4R and 8R. 4R and 6R.
b
2.4
Member of 6R and 8R. c Member of
When we calculated the potential curves using various model potentials for comparison with the corresponding reference DFT results, no periodic boundary conditions and long-range corrections were applied: each atom-atom pair energy was calculated from the terms of the direct pair (Lennard-Jones and Coulomb) interactions. Separate classical molecular simulations were also performed for comparison of the most probable positions of adsorbate molecules with the DFT optimized geometries. The most probable positions in the simulations were determined from site-site pair correlation functions. The latter were obtained from simulations with one adsorbate particle (in the grand canonical simulation, the chemical potential was adjusted to ensure that only one particle is present, on average, in the simulation box). We calculated the distances between the sites of the adsorbent and adsorbate molecules from the first (and possibly second) peaks of the pair correlation functions, and these were used in the comparison with the DFT results. 3. Model Potentials Our model for the zeolite framework is semiflexible: the framework atoms are fixed at the atomic positions taken from X-ray diffraction experiments,23 and the nonframework Na+ ions are allowed to move. Zeolite NaA-4, as other zeolites, correspond to the general formula Mx/n[(AlO2)x(SiO2)y], where the nonframework cations (the valence of which is n) are denoted by M. The primary unit of the zeolite structure is the TO4 tetrahedron, which can have either Si or Al in its center. These units are generally organized into secondary units called rings. Zeolite NaA-4 is of framework type LTA (Linde type A) and has three types of rings: 4R, 6R, and 8R containing 4, 6, and 8 O atoms, respectively. The interconnection of 4R and 6R rings forms nearly spherical cages, which are called sodalite cages. The O bridges connect the sodalite cages to each other so that they form supercages23 with a diameter of about 1.2 nm. The crystal structure of zeolite NaA-4 belongs to the fm3c space group with a lattice parameter of 2.4555 nm. The MC simulation box consisted of 576 framework atoms (96 Si, 96 Al, 384 O atoms) and 96 Na+ ions. In the model, each AlO4 tetrahedron is connected to a SiO4 tetrahedron to fulfill the Lo¨wenstein rule that prohibits Al-O-Al linkages. The potential models are defined by the positions of the
methanol36 CH3 105.2 0.374 0.265 O 86.5 0.303 -0.700 H 0.435 δ(O-H) ) 0.0945 nm, δ(O-CH3) ) 0.1425 nm R(H-O-CH3) ) 108.53° ethanol35 CH3 104.17 0.3775 CH2 59.38 0.3095 0.265 O 85.55 0.3070 -0.700 H 0.435 δ(O-H) ) 0.0945 nm, δ(O-CH2) ) 0.1430 nm δ(CH2-CH3) ) 0.1530 nm R(CH3-CH2-O) ) 108°, R (H-O-CH2) ) 108.5° θ(H-O-CH2-CH3) ) 180°
interaction sites and their parameters, namely, the Lennard-Jones energy (ε) and size (σ) parameters and the point charges (q). Earlier, we used the semiflexible potential model proposed by Lee et al.6 (model A). In all simulations of adsorption, we found that this model predicts higher equilibrium loadings for pure adsorbates, such as methanol and ethanol, than the corresponding experimental results. Therefore, we changed to an other literature model, the potential model of Faux et al.,2 and the model was modified by optimizing its parameters to the equilibrium adsorption of pure water, methanol, and ethanol (model B). The parameters of models A and B are given in Table 1. The main changes in model B with respect to model A are (1) the assumption that all oxygen atoms in the framework are of the same type (i.e., each O atom has the same parameter set), (2) the significantly larger charge parameter of Na+ ions (being set to 1), and (3) reduction of the LJ parameters for Si and Al to 0 (i.e., these atoms participate only in Coulomb interactions). Although the charge parameters were adopted from ref 2 only with slight modifications, it is important to emphasize that the difference between the charge parameters of Si and Al in model B was set equal to the charge parameter of Na+ ions The model parameters were optimized for the equilibrium amount of adsorption of water, methanol, and ethanol (normal liquids) and of the much less polar hydrogen sulphide (normal gas). As a result, the simulations of adsorption of water-alcohol mixtures produced better results for both the amount of adsorption and the selectivity data using model B than with model A. We tested model B for other adsorbates. For carbon dioxide adsorption, we found that the results do not match the corresponding experimental data.14 This suggests that model B is too polar in the simulations with nonpolar adsorbate molecules, suggesting that the model needs to be further improved. The main problem is the selection of the charge parameter of the Na+ ion because the polarity of the model strongly depends on it, significantly affecting the adsorption behavior of the zeolite. Following a detailed structure examination of models A and B, we assumed that model B could be improved by decreasing the charge parameter of Na+ ions (which calls forth some changes in the LJ parameters). We found that this parameter should be between 0.55 and 1, and the charge parameters of the framework atoms need to be changed accordingly to fulfill electroneutrality. We
Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4
J. Phys. Chem. C, Vol. 113, No. 28, 2009 12229
TABLE 3: Simulated and Experimental Data13 for the Amount of Adsorption (na)a
peaks and a slight dissimilarity (0.01-0.03 nm) in the position of peaks, especially in the O1-H(adsorbate) interactions. For model C, the adsorbent-adsorbate interactions seem to be weaker than those for model B. In all three zeolite models, the positions of the first O1-H(adsorbate) peaks indicate the existence of hydrogen bonding between the adsorbent and the adsorbates. Even though the LJ distance parameter of O is smaller in model C, these curves suggest that the adsorbates cannot approach the zeolite framework atoms as close as in model B. The comparatively featureless O3(zeolite)-H(adsorbate) pair correlation functions are not presented, as they do not provide any new information. From the Na(zeolite)-O(adsorbate) pair correlation functions plotted in Figure 3, we can observe that the positions of the first peaks are more or less identical in the case of models A and C, and the relative intensity of these interactions is nearly the same for models B and C. These results justify what we have expected, namely, that the performance of model C is between those of models A and B. Comparison of the Model Potentials with DFT. As we mentioned above, in models B and C there is only one parameter set for O instead of three parameter sets as in model A, and it is necessary to examine in detail the influence of this dissimilarity. We compared the potential energy curves as a function of the location of the Na+ ion calculated with the model potentials and with DFT. The quantum chemical results can be considered as a reliable reference because the minimum energy positions of Na+ ions according to the DFT calculations-one close to the center of the 6R ring and one near the center of the 8R ring-are the same as those in the X-ray diffraction experiments.23 Three characteristic straight lines were chosen along which the Na+ ions are moved in the zeolite fragment, starting from the potential energy minima corresponding to the Na+ ions being in the 8R or in the 6R ring. We have not considered positions of Na+ ions in the 4R ring or very near to it. Accordingly, the lines along which the potential energy is sampled run (1) perpendicular to the 8R plane, passing through the center of the 8R ring (parallel to the x axis in Figure 1), (2) along one of the medians of the 8R ring (the z axis in Figure 1), and (3) perpendicular to the 6R plane, passing through the center of the ring. These potential curves are the ones that we used to guide our potential parameter optimization, as described in the Introduction. In Figure 4a-c, we plotted the potential energy curves obtained with the DFT method and with the three model potentials. The energy along each curve is referred to the respective energy minimum. We found that the energy minima are obtained at the same position with each method along the line perpendicular to the 8R plane through the center of the 8R ring (see Figure 4a). It is seen that the curve of model C runs the closest to the DFT results. If the Na+ ion goes along the diameter of the 8R ring (see Figure 4b), the shapes of the curves are very similar but the minimum of model A is closer to the QC minimum than those of the other two models (for which the location of the minimum is nearly the same). The reason for this difference is that the LJ size parameters of the O atoms and Na+ ions are larger in models B and C because, in these models, the consequence of the omission of the LJ interactions for the framework atoms, Al and Si, had to be counterbalanced by the LJ size parameters of O atoms. This difference appears in Figure 4c, too, where the Na+ ion passes through the center of the 6R ring. There is a local maximum on the curves at the DFT minimum for all models; however, the curve of model A fits the DFT results better. The appearance of the maximum is an artifact: the diameter of the 6R ring is
number of adsorbate molecules per unit cell (na)b water
methanol
ethanol
p
1 kPa
0.1 kPa
1 kPa
0.1 kPa
1 kPa
0.1 kPa
experimental model A model B model C
195 203(2) 197(2) 187(2)
184 182(2) 185(2) 173(2)
76 95(1) 77(1) 77(1)
70 93(1) 73(1) 73(1)
46 63(1) 49(1) 48(1)
44 60(1) 48(1) 47(1)
a Simulated data for models A and B are taken from ref 7. Numbers in parentheses represent statistical uncertainties of the simulations in the last reported digits. b
selected for the new model the Na+ and Al parameters proposed by Maurin et al.:33 qNa ) 0.7 and qAl ) 1.7 (Table 1). We then varied the LJ parameters, and with each, the potential curves were calculated at the locations of the Na+ ion used in the DFT calculations (the Na+ ion moving along straight lines). This is a very quick calculation (requires only fragments of a second) in contrast to the full MC simulations (requiring several days). By comparing to the DFT data, it was easy to discard unrealistic combinations of the charge and LJ distance parameters that poorly match the corresponding DFT data. When a promising parameter set was found, we performed classical test simulations and calculated the equilibrium adsorption data for pure water, methanol, and ethanol. During model optimizations, the LJ size parameters of model B decreased and the energy parameters increased to some extent. The model with the resulting parameter set will be called model C in this paper. In the simulations, we used the well-known SPC/E34 model, the OPLS35 model, and an OPLS-type model36 for the representation of the water, ethanol, and methanol molecules, respectively. These models are suitable for the reproduction of thermodynamic properties of the pure substances. The interaction parameters as well as the geometry parameters of the adsorbates are summarized in Table 2. 4. Results and Discussion Performance of Model C in Simulations of Adsorption of Pure Substances in Zeolite NaA-4. For the final parametrization of model C, the experimental adsorption data for pure water, methanol, and ethanol at T ) 298 K13 were used. The equilibrium adsorption loadings at p ) 1 and 0.1 kPa are shown in Table 3 (some data are taken from our previous work7). Comparing the three models, we can say that models B and C yield a reasonably good reproduction of the experimental loadings for all three adsorbates. The results for models B and C are quite similar, especially for methanol and ethanol. The quality of reproduction of the experimental data for model A is not as good as that of models B and C and becomes worse in the order of water, methanol, and ethanol. From pair correlation functions, we can get information about the relative strength of the different interactions and the features of hydrogen bonding. As the main differences between models B and C are associated with the Na+ ions and the O atoms, the O(zeolite)-adsorbate and the Na(zeolite)-adsorbate interactions are the key properties to be examined. Three types of O atoms are determined by site symmetries: O1 is a member of a 4R and an 8R ring, O2 is a member of a 6R and an 8R ring, and O3 is a member of a 4R and a 6R ring. As expected, the pair correlation functions plotted in Figure 2 for the adsorbate’s H atom and the framework O atoms are quite similar for models B and C. There is a significant difference in the intensity of
12230
J. Phys. Chem. C, Vol. 113, No. 28, 2009
Csa´nyi et al.
Figure 2. O1(zeolite)-H(adsorbate) and O2(zeolite)-H(adsorbate) pair correlation functions with one adsorbent molecule of water, methanol, and ethanol at 378 K. The dashed, dotted, and solid lines correspond to models A, B, and C, respectively.
relatively small, and even a slight enhancement of the size parameters can cause a large difference in the potential energies if the framework atoms of the 6R ring are close to the Na+ ion.
The consequence is that, in contrast to the DFT and X-ray results, the most probable position of the Na+ ions in the classical simulations is not exactly the center of 6R ring but
Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4
J. Phys. Chem. C, Vol. 113, No. 28, 2009 12231
Figure 3. Na(zeolite)-O(adsorbate) pair correlation functions with one adsorbent molecule of water, methanol, and ethanol at 378 K. The dashed, dotted, and solid lines correspond to models A, B, and C, respectively.
slightly outside of its plane. This is, in a sense, the price we have to pay to get a much better simulation performance of a macroscopic system. We also tested how realistic are the zeolite models combined with the zeolite-adsorbate interaction by finding the most probable positions of the adsorbate molecules, water, methanol, and ethanol. The reference QC calculations were performed with the selected zeolite fragment (Figure 1) supplemented by two Na+ ions, as described in section 2. After the initial constrained optimizations with rigid adsorbates, the geometry of water, methanol, or ethanol was relaxed, and the optimized structures showed negligible deviations from the structures used in simulations (the largest change was 2 pm in the bond lengths, 3° in bond angles, and 2° in torsional angle). The position of the Na+ ions did not change significantly either. One of them always stayed in its most favorable position, that is, in the center of the 6R ring. The other Na+ ion followed the movement of the adsorbate molecule, but it stayed essentially on the center line of the 8R ring. The optimized positions of the adsorbates are shown in Figure 5. The Na+ ion helps the adsorption of the adsorbate molecules through an orientation effect: the O atom of the adsorbate tends to take a position close to the Na+ ion in the 8R ring in such a way that the H atoms of the adsorbate’s OH group point toward the O atoms of the framework. The calculations show clearly that, unlike methanol and ethanol molecules, the water molecule is bound to the zeolite framework by two hydrogen bonds. The presence of the additional hydrogen bond is reflected by the binding energies of the adsorbates to
the zeolite, which are 168.5, 115.2, and 111.2 kJ mol-1 for water, methanol, and ethanol, respectively. Although we do not think these numbers can be considered very accurate, the tendency is certainly correct. The significantly larger binding energy of water as compared to that of the alcohols justifies why water adsorption is preferred on this type of zeolite, which causes improved selectivity of water to alcohols. For the zeolite models A, B, and C, the most probable positions of adsorbate molecules were determined from grand canonical simulations using one adsorbate particle, where the first peaks of the calculated pair correlation functions (for all pairs of interaction sites) provided the distances between different adsorbent-adsorbate atom pairs. Shown in Figure 6 is the comparison of the simulation versus the DFT results for the distances between the O and H atoms of the adsorbate and the nearest O, Si, and Al atoms and the Na+ ion (located in the 8R ring) of the zeolite. The figures show the deviation of the distance obtained in molecular simulations from that found by quantum chemical geometry optimization. Figure 6a shows that models B and C describe well the water-zeolite interactions; the overall performance of model A is somewhat poorer. The distances obtained with model C are quite similar to those from model B. In the case of the Na+-O(water) interaction, model B is somewhat better. It is seen that the results for model A are also in good agreement with the DFT results, except for the O(zeolite)-H(water) and Na(zeolite)-O(water) interactions. Figure 6b illustrates that model C provides the best overall agreement with the DFT calculations in the case of methanol.
12232
J. Phys. Chem. C, Vol. 113, No. 28, 2009
Csa´nyi et al.
Figure 5. Optimized positions of water (a), methanol (b), and ethanol (c). The Si, Al, O, H, and Na+ atoms are denoted by yellow, gray, red, white, and blue, respectively.
Figure 4. Potential energies for different positions of one Na+ ion: (a) when the ion is translated along a line perpendicular to the plane of ring 8R (the line parallel to the x axis in Figure 1), passing through the center of 8R (the negative direction corresponds to locations outside of the supercage, i.e., the negative z axis in Figure 1), (b) when the ion is translated along one of the medians of the 8R ring (the z axis in Figure 1), and (c) when the ion is translated along a line perpendicular to the plane of ring 6R passing through the center of 6R (a line in the x-z plane in Figure 1). The energy is measured from the lowest energy of each curve. The solid blue line corresponds to the quantum chemical calculation; the dashed, dotted, and solid lines correspond to models A, B, and C, respectively.
The results of model B are quite similar to those of model C, except for the Na(zeolite)-O(methanol) distances, where model B produces a more than 0.1 nm worse overestimate. It is visible that model A gives the worst results for the different interatomic distances, except that between Na(zeolite) and O(methanol). In
the case of ethanol (Figure 6c), models B and C describe the interactions between the framework atoms and the adsorbate molecules quite well; only the Na(zeolite)-O(ethanol) distance is too large for model B. From Figure 6, we can make the conclusion that, overall, model C predicts more realistically the behavior of the system for the adsorption of water, methanol, and ethanol. 5. Simulation of Adsorption of Water-Alcohol Mixtures on Zeolite NaA-4 Our main purpose with the development of model C is an adequate modeling of the adsorption selectivity for wateralcohol mixtures. We performed such simulations at T ) 378 K. Because this temperature is above the normal boiling points of the components, the maximum pressure was set to p ) 100 kPa, the pressure of the available gas permeation experiments for water-alcohol mixtures.13 A detailed comparison of models A and B was made in our earlier work, simulating at different pressures and mole fractions of bulk phase.7 For model C, we carried out simulations at equimolar composition of the bulk
Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4
J. Phys. Chem. C, Vol. 113, No. 28, 2009 12233
Figure 6. Distances between the atoms of the zeolite framework and the O (Ow, Om, Oe) as well as H (Hw, hydroxyl-Hm, hydroxyl-He) atoms of a (a) water, (b) methanol, or (c) ethanol molecule at the lowest energy position with respect to the zeolite, as obtained from MC simulations compared with the corresponding DFT values. The distances for the MC simulations are equal to the location of the first or highest maximum of the pair correlation functions. The dashed, dotted, and solid lines correspond to models A, B, and C, respectively.
water-alcohol mixtures. We present here the comparison of the simulation results obtained with the three model potentials under the latter conditions. The comparison of the results of mixture adsorption (shown in Figure 7) clearly indicates the difference between the three
models. Model A predicts a reversal of equilibrium selectivity at relatively high pressures (around 20 kPa) in both water/ methanol and water/ethanol mixtures: at higher pressures, water adsorption is preferred, whereas at lower pressures, methanol or ethanol adsorption is preferred. In the case of model B, the
12234
J. Phys. Chem. C, Vol. 113, No. 28, 2009
Csa´nyi et al.
Figure 7. Simulation results for the amount of adsorption (na) for equimolar mixtures of water-methanol and water-ethanol on zeolite NaA-4 at 378 K. Dashed, dotted, and solid lines represent models A, B, and C, respectively. The circles and squares correspond to water and methanol or ethanol, respectively. The statistical uncertainties of the calculated loadings do not exceed the symbol size. The results for model A are taken from ref 8.
inversion of selectivity cannot be observed even at very low pressure (see ref 7). The results of model C differ from those of the other two models; namely, at low pressures (around 1 Pa), the amount of adsorbed water and methanol/ethanol is approximately equal. In addition, model C shows the best selectivity rate at higher pressures (above p ) 1 kPa). Note that zeolite NaA-4 is applied industrially at higher pressures (around p ) 140 kPa), for example, for ethanol dehydration during the high-efficiency pressure swing adsorption process. Experimental data are available at ambient pressure only from membrane permeation experiments,13 the results of which depend not only on the difference in adsorptivity but also on the diffusivity of the components. As a result, no direct comparison can be made with our simulation results. Nevertheless, we can assume that model C predicts the amount of adsorption better than the earlier models used for water/methanol and water/ethanol mixtures because the experiments show very high selectivity of water to alcohols. The observation mentioned in section 4, namely, that, according to the DFT calculations, the binding energy of water to the zeolite is much larger than that of alcohols, supports that the high selectivity predicted by model C is realistic. 6. Conclusion We developed a new potential model for zeolite NaA-4 based on a new paradigm. During parameter optimization, we utilized potential curves obtained in quantum chemical calculations using density functional theory for judging the quality of the model potentials. In contrast to an often used approach, instead of fitting the model potential to these curves, we used the DFT results to screen parameter sets that occurred during model development and discarded the ones that produced unrealistic potential curves. The quantum chemical reference calculations compared with the results of classical MC simulations of two published models (models A6 and B7) suggested a possible way to improve the potential model for zeolite NaA-4. Namely, the charge of the Na+ ion should be set lower than 1 (and, of course, such a change involves the modification of the rest of the pair potential parameters). Considering the results of structure examinations, we proposed a new parameter set starting from the later model (model B). When the calculated equilibrium amounts of adsorption obtained in simulations using the three models at T
) 298 K are compared with available experimental data, it can be established that the new model provides a reasonably good reproduction of the experimental loadings for water, methanol, and ethanol. The simulation results for water-methanol and water-ethanol mixtures showed that, at higher pressures (above p ) 1 kPa), the newer model (model C, which, in all other aspects, is more realistic than the others) predicts a larger relative amount of adsorption of water than that of the earlier models. This result may be of great importance considering that zeolite NaA-4 is applied industrially at higher pressures. The higher selectivity of water to alcohols at higher pressures agrees better with the membrane permeation experiments.13 The higher selectivity is in good agreement with the larger binding energy found in the DFT calculations for the water to zeolite versus the alcohol to zeolite connections. Our earlier studies concerning the adsorbent-adsorbate pair correlation functions suggested that model B is more realistic than model A,7 so here, we compare the results of model C to those of model B. We found the pair correlation functions of models B and C to be quite similar, except that the intensity of the peaks differs. The lower peak intensities for model C indicate weaker adsorbent-adsorbate interactions. The O1(zeolite)H(adsorbate) pair correlation functions of the newer model suggest that the adsorbates, even though the LJ distance parameter of O(zeolite) is smaller in model C, cannot approach the zeolite framework atoms as close in this model as in the case of the other two models. As it can be expected on the basis of the model parameters, the Na(zeolite)-O(adsorbate) pair correlation functions (especially the first peaks) also show that the performance of model C lies between those of models A and B. Examination of the potential curves calculated as a function of a Na+ ion with respect to a fragment of zeolite proved that the newer model describes the framework-Na+ interactions better than our previous models and even coincides, almost perfectly, with one of the quantum chemical curves. These results also confirm that, in the course of model development, we managed to partially counterbalance the unrealistic components of the model potential (namely, that the LJ parameters for Si and Al were set to 0). The comparison of the optimum positions of the adsorbate molecules water, methanol, and
Water-Methanol (Ethanol) Mixtures in Zeolite NaA-4 ethanol obtained in quantum chemical calculations with those derived from molecular simulations also supports that the new model presented in this paper behaves more realistically than all earlier ones. Acknowledgment. This work was supported by the Hungarian Scientific Research Fund (Grant Nos. OTKA K 75132 and K 77938). References and Notes (1) Abascsal, J. L. F.; Vega, C. J. Phys. Chem. C 2007, 111, 15811– 15822. (2) Faux, D. A.; Smith, W.; Forester, T. R. J. Phys. Chem. B 1997, 101, 1762. (3) DeLara, E. C.; Kahn, R.; Goulay, A. M. J. Chem. Phys. 1989, 90, 7482. (4) Bandyopadhyay, S.; Yashonath, S. Chem. Phys. Lett. 1994, 223, 363. (5) Jameson, C. J.; Jameson, A. K.; Baello, B. I.; Lim, H. M. J. Chem. Phys. 1994, 100, 5965. (6) Lee, S. H.; Moon, G. K.; Choi, S. G.; Kim, H. S. J. Phys. Chem. 1994, 98, 1561. (7) Rutkai, G.; Csa´nyi, E´.; Kristo´f, T. Microporous Mesoporous Mater. 2008, 114, 455–464. (8) Kristo´f, T.; Csa´nyi, E´.; Rutkai, G.; Mere´nyi, L. Mol. Simul. 2006, 32, 869–875. (9) Wu, G.-S.; Schatz, G. C.; Lendvay, G.; Fang, D.-C.; Harding, L. B. J. Chem. Phys. 2000, 113, 3150. (10) Pederson, L. A.; Schatz, G. C.; Ho, T.-S.; Hollebeek, T.; Rabitz, H.; Harding, L. B.; Lendvay, G. J. Chem. Phys. 1999, 110, 9091. (11) Higgins, C. J.; Chapman, S. J. Phys. Chem. A 2004, 108, 8009. (12) Bernshtein, V. B.; Oref, I. J. Phys. Chem. A 2000, 104, 706. (13) Okamoto, K.; Kita, H.; Horii, K.; Tanaka, K.; Kondo, M. Ind. Eng. Chem. Res. 2001, 40, 163. (14) Ahn, H.; Moon, J.; Hyun, S.; Lee, C. Adsorption 2004, 10, 111– 128. (15) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (16) Hohenberg, P.; Kohn, W. Phys. ReV. 1964, 136, B864. (17) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (18) Veszpre´mi, T.; Fehe´r, M. A kVantumke´mia e´s alkalmaza´sai (in Hungarian); Mu˝szaki Ko¨nyvkiado´: Budapest, Hungary, 2002.
J. Phys. Chem. C, Vol. 113, No. 28, 2009 12235 (19) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (20) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (21) Schaefer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (22) 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.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.05; Gaussian, Inc.: Wallingford, CT, 2004. (23) Pluth, J. J.; Smith, J. V. J. Am. Chem. Soc. 1980, 102, 4704. (24) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8254. (25) de Leeuw, S. W.; Perram, J. W.; Smith, E. R. Proc. R. Soc. London, Ser. A 1980, 27, A373. (26) Demontis, P.; Spanu, S.; Suffritti, G. B. J. Chem. Phys. 2001, 114, 7980. (27) Siepmann, J. I.; Frenkel, D. Mol. Phys. 1992, 75, 59. (28) Kofke, D. A.; Glandt, E. D. Mol. Phys. 1988, 64, 1105. (29) van′t Hof, A.; de Leeuw, S. W.; Peters, C. J. J. Chem. Phys. 2006, 124, 054905. (30) Furukawa, S.; Goda, K.; Zhang, Y.; Nitta, T. J. Chem. Eng. Jpn. 2004, 37, 67. (31) Jaramillo, E.; Chandross, M. J. Phys. Chem. B 2005, 108, 20155. (32) Mu¨ller, J. C. M.; Hakvoort, G.; Jansen, J. C. J. Therm. Anal. 1998, 53, 449. (33) Maurin, G.; Llewellyn, Ph.; Poyet, B.; Kuchta, B. J. Phys. Chem. B 2005, 109, 125. (34) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (35) Jorgensen, W. L. J. Phys. Chem. 1986, 90, 1276. (36) van Leeuwen, M. E.; Smit, B. J. Phys. Chem. 1995, 99, 1831.
JP902520P
