ARTICLE pubs.acs.org/JPCC
Alcohol Adsorption onto Silicalite from Aqueous Solution Ruichang Xiong, Stanley I. Sandler,* and Dionisios G. Vlachos Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716, United States ABSTRACT: A methodology based on the combination of grand canonical Monte Carlo (GCMC) and expanded ensemble (EE) simulations has been used to study the adsorption of alcohols from dilute aqueous solution onto the silicalite zeolite. The chemical potential of the guest alcohol molecules in bulk aqueous solutions is calculated by the EE method, and the adsorption isotherms of the alcohols are calculated using GCMC simulations. This approach results in adsorption isotherms that relate the loadings to external concentrations without using an equation of state, experimental data, or performing computationally expensive simulations such as two-phase Gibbs ensemble Monte Carlo (GEMC). Also, the method is “force-field-consistent”, in that the same force field is used for the bulk fluid and the adsorption calculations, and the resulting isotherms can be directly compared with experimental data. We have established the validity of the method used by comparing the calculated adsorption isotherms with experimental data for methanol and ethanol on silicalite, both from aqueous solution and from the vapor phase. Calculated heats of adsorption are also in reasonable agreement with experimental data. In addition, these simulations provide information on the location of the alcohol adsorption sites. We also observe that the presence of water increases the adsorption of alcohols at dilute concentrations onto silicalite; therefore, an accurate estimate of amount of water adsorbed in the nanoporous substrate is important to the prediction of alcohol adsorption from aqueous solution. The method presented here can be used to study the adsorption of any molecular species on any nanoporous adsorbent.
1. INTRODUCTION The effective separation of organic molecules from dilute aqueous solutions is of importance in many industry applications. Two examples are fermentation for chemicals production from biomass1 and the removal of hazardous organic materials from wastewater.2,3 Distillation is widely used for separations, but it is energy-intensive and is not always applicable.4 One energetically efficient alternative is adsorption onto porous materials. Zeolite-based adsorption plays an important role in separations and catalysis.5,6 Zeolites are nanoporous crystalline materials with well-defined pore structures. They have been considered as adsorbents, catalysts, ion-exchange agents, and sensors due to their excellent adsorption and molecular sieve properties.7,8 They selectively adsorb molecules from mixtures based on shape selectivity and also by their amphoteric nature, hydrophobicity, or hydrophilicity, which depends on the type of zeolite and the alumina content in the zeolite framework. Their polar nature makes them a good choice for separating organic/water mixtures. Over the past decade, several experimental studies have examined the adsorption of alcohols and polyols on zeolites. Milestone and Bibby9 reported the adsorption and concentration of C1C5 alcohols from dilute aqueous solutions and found that ethanol was concentrated from 2 to 35% and butan-1-ol was concentrated from 0.5% to 98% by adsorption on silicalite from aqueous solution, and they suggested that silicalite could be used instead of distillation. Cekova et al.10 investigated the removal of organic molecules from aqueous solutions in a series of experiments with lower alcohols and showed that silicalites could be used for the removal and recovery of these organic molecules. Nomura et al.11 studied the transport of ethanol/water mixtures r 2011 American Chemical Society
through a silicalite membrane by pervaporation and vapor permeation and reported that silicalite has high ethanol selectivity. Other studies included Dubinin et al.,12 Long et al.,13 and Mallon et al.14 These experimental studies have provided useful information on the adsorption and separation characteristics of zeolites for alcohols. Molecular simulations are useful in providing data and information on adsorption thermodynamics, on how the adsorption sites change with increased loading, and on diffusion. Typically, grand canonical Monte Carlo (GCMC) simulation is used to compute adsorption isotherms and heats of adsorption, and molecular dynamics (MD) simulation is used to study the diffusion. Gupta et al.15 studied pure methanol adsorption on silicalite using the all-atom OPLS (OPLS-AA) force field16 for the methanolmethanol interaction and adjusting both the van der Waals (VDW) and Coulombic interactions between methanol and silicalite to match loadings at low pressures. Blanco and Auerbach17 investigated methanol diffusion in NaY, DAY, and silicalite zeolites. Plant et al.18 computed methanol diffusion in the NaY zeolite using a modified Blanco and Auerbach17 force field. Kristof and co-workers19 predicted the adsorption of watermethanol and waterethanol mixtures in the NaA zeolite using modified united-atom OPLS (OPLS-UA) force fields.20 Kuhn et al.21,22 performed both experiments and simulations of the adsorption and diffusion of water, methanol, and ethanol on the hydrophobic all-silica decadodecasil 3R (DD3R) zeolite. An Received: June 6, 2011 Revised: August 23, 2011 Published: August 23, 2011 18659
dx.doi.org/10.1021/jp205312k | J. Phys. Chem. C 2011, 115, 18659–18669
The Journal of Physical Chemistry C interesting observation from their study is that the loadings of both methanol and ethanol at low pressures (fugacities) are slightly enhanced by the presence of water. Krishna and van Baten23 adopted the same force field as that used by Kuhn et al.21 to study the hydrogen bonding effects in the adsorption of wateralcohol mixtures in several zeolites and found that hydrogen bonding effects are stronger in wateralcohol mixtures than in the pure constituents. In GCMC simulations, the chemical potential or fugacity of the adsorbent is fixed in the bulk phase that at equilibrium is equated with that of the adsorbed phase. Therein lies the problem that we address in this paper, which is that GCMC simulations give loading as a function of imposed chemical potential while what is measured experimentally is loading versus solution composition or pressure. All of the GCMC simulations mentioned above were performed for low-pressure gas-phase adsorption, where the chemical potential of the adsorbate is easily related to the pressure of the pure component or its partial pressure in the bulk gas by the ideal gas law, so that the chemical potential set in the simulations can be directly related to the pressure measured in experiment. However, for adsorption from an aqueous or other liquid phases, or from a high-pressure gas phase, the interrelationship between the measured composition (concentration) and the chemical potential set in GCMC simulation is more complicated. Consequently, indirect methods have been used that bear no relation to the force field used in the simulations. For example, for their GCMC simulations, Chempath and Snurr24 used an equation of state (EOS) to convert the composition of each component in a liquid binary system of aromatics to the partial pressure of each component in the gas phase. Yazaydin25 also used EOS to calculate the fugacity of each component in the gas phase for a liquid methyl tertiary butyl ether (MTBE)water and 1,4-dioxanewater binary system. Narasimhan et al.26 used experimental data to estimate the infinite-dilution chemical potential of p-cresol in an ideal dilute solution and then used this to fix the bulk chemical potential to be used in GCMC simulations. However, the needed experimental data are not always available, and EOS predictions can be inaccurate for the calculation of the partial pressures from liquid compositions, especially of aqueous alcohol mixtures. Further, there is an inherent inconsistency of such calculations because bulk fluid simulations with the same force field as that used for the adsorption simulations are unlikely to lead to the same chemical potentialcomposition relationship as that obtained from an EOS or infinite-dilution measurements. The use of different interrelations between the chemical potential and concentration or pressure (that is, a relationship obtained with the force field used in simulation and a different relationship obtained from a macroscopic model such as an EOS or an activity coefficient model extrapolated from infinite-dilution behavior) will result in a displacement of the simulated adsorption isotherms with respect to concentration or pressure depending on which interrelationship is used. For internally consistent predictions of adsorption isotherms, the same force field for the adsorbing molecules should be used for all phases and in all parts of the calculations. In principle, one way to avoid this inconsistency is to use the Gibbs ensemble Monte Carlo (GEMC) simulation method developed by Panagiotopoulos27 to calculate the equilibrium for an adsorbate (either as a pure component or a component in a mixture) in the bulk and adsorbed phases. This method has been applied successfully to simulate the vaporliquid phase diagram for
ARTICLE
many compounds, including alcohols.28 However, it is computationally difficult for simulating the sorption of organic molecules from aqueous solution into a zeolite because molecule transfer between dense liquid solution and zeolite is extremely inefficient, even when combined with the configurational-bias technique29 to enhance the efficiency of sampling. Further, for adsorbing solutes from dilute aqueous solution, large systems of the dense phases are needed for both the aqueous and zeolite-adsorbed phases, making equilibration even more difficult and time-consuming. Here, we present a technique that can be used to predict internally consistent adsorption thermodynamics and illustrate its use with application to the adsorption from vapor and from dilute aqueous solution of two alcohols onto the silicalite zeolite. This method does not use an EOS or experimental data to relate the bulk chemical potential to concentration and is force-fieldconsistent in that the same force field model is used throughout. The object of this paper is to outline and illustrate the method that can be used for all adsorption problems of the type considered. We use the adsorption of methanol and ethanol from aqueous solutions as illustrative examples and validate the method by comparing the results of simulation with experimental data. Isodensity surfaces are also calculated to identify the adsorption sites of alcohols in silicalite. The method can be readily extended to the adsorption of other compounds from the dilute solution onto other adsorbents.
2. METHODOLOGY 2.A. Thermodynamic Background. GCMC simulates a system at constant chemical potential, volume, and temperature. During the simulation, the number of adsorbed molecules fluctuates, and one calculates the average loading. A requirement for equilibrium between adsorbed and bulk solution phases is the equality of adsorbate chemical potentials sol μads A ¼ μA
ð1Þ
where μads A is the chemical potential of the adsorbate molecule in the adsorbent framework and μsol A is the chemical potential of the adsorbate (solute) molecule in the bulk solution. As already mentioned, the information needed to compare the results of simulation with experiment is the relation between the chemical potential used in GCMC simulations and the bulk phase pressure or concentration, which are the measurable variables. The chemical potential of the solute in solution is30 3 1 id μsol A ¼ μA þ μ~A ¼ kb T lnðFA ΛA qA Þkb T lnÆexp½βUðR0 Þæ0
ð2Þ where T is the temperature, kb is the Boltzmann constant with β = (kbT)1, FA is the number density of the solute, ΛA is the de Broglie wavelength of the solute, qA is the internal degrees of freedom partition function of the solute, U(R0) is the intermolecular energy of a solute molecule at a fixed position R0 with the rest of the system, Ææ0 is the ensemble average, μid A is the chemical potential in the ideal gas state, and μ ~A is the excess chemical potential or the work (free energy) of coupling a solute molecule to the other molecules in solution. Here, we are interested in adsorption from dilute solutions, which is how experiments are often performed. The mole fraction, xA, in the dilute solution is given by F ð3Þ xA ¼ A F 18660
dx.doi.org/10.1021/jp205312k |J. Phys. Chem. C 2011, 115, 18659–18669
The Journal of Physical Chemistry C
ARTICLE
where F is the number density of solution taken to be the solvent number density for a solution dilute in the solute. Substituting eq 3 into eq 2, we obtain 3 1 ∞ μsol A ¼ μA þ kb T lnðFΛA qA Þ þ kb T ln xA
ð4Þ
~A, which is the excess chemical potential at where μ∞ A = limxAf0 μ infinite dilution of the solute. In the limit of dilute solutions, a solute molecule is surrounded by solvent molecules. Thus, the excess chemical potential at infinite dilution is the coupling work of placing a solute molecule into the pure solvent. Equation 4 provides the chemical potential of the adsorbate (solute) molecule at composition xA in dilute solutions. At higher solute concentrations, one would have to include an activity coefficient based on the Henry’s law standard state in the last term in eq 4. Further, assuming that the internal degrees of freedom do not change between the solution and gas phases, the chemical potential of the solute in any mixture of the real gas is ! 0 Λ3A q1 0 ads id A μA ¼ μA þ μ~A ¼ kb T ln þ kb T ln PA þ μ~A kb T
solution (in which there are solutesolvent interactions but no solutesolute interactions). 2.B. Direct Calculation of the Infinite-Dilution Excess Chemical Potential: Expanded Ensemble (EE) Method. The chemical potential of the adsorbed species in the bulk fluid in eq 4 can be obtained at a specified mole fraction if the excess chemical potential at infinite dilution is known and the concentration is sufficiently dilute that the value of the Henry’s Law activity coefficient is unity. The solvation free energy is the excess chemical potential of transferring a solute from the ideal gas to the real solution at the same temperature and concentration. For a reference state at constant temperature and volume, the excess chemical potential can be computed as follows μex ðT, V Þ ¼ kb T lnÆexp½βUðR0 ÞæN, V , T
In dilute solution, U(R0) is the intermolecular potential energy between a solute at a fixed position R0 and the pure solvent. Alternatively, if the reference state is constant temperature and pressure, the excess chemical potential is computed from μex ðT, PÞ ¼ kb T ln
ð5Þ where PA is the ideal gas pressure contribution of the adsorbate (solute) to the total pressure of the system, and μ ~A0 is the work of turning an ideal gas adsorbate molecule into an interacting molecule. Using μ ~A0 = kbT ln ϕA and fA = PAϕA, where fA is the fugacity of the adsorbate (solute) and ϕA is the fugacity coefficient, we then obtain the fugacity-based chemical potential of the adsorbate in any real gas mixture ! Λ3A q1 ads A ð6Þ μA ¼ kb T ln þ kb T ln fA kb T If the gas mixture density is sufficiently low (e.g., at low pressure), there is no coupling work between gas molecules, and thus, ϕA = 1, and the fugacity equals the pressure. In molecular simulations for adsorption from the gas phase, we specify the chemical potential at each pressure, while for adsorption from solution, we need to relate the chemical potential to composition. From eqs 4 and 6, the composition is related to the fugacity by ∞ μ ð7Þ fA ¼ Fkb T exp A xA kb T The same solute force field should be used in the gas and solution phases; otherwise, there will be an inconsistency in the relationship between the chemical potential and composition of the phases. The Henry’s Law constant in the solution phase (HLCS), KHS, is defined as the infinite dilution limit of the fugacity/concentration ratio of a solute in a solution31 and has the dimensions of pressure fA KH ¼ lim ð8Þ xA f 0 xA By comparing eqs 7 and 8, we obtain KHS ¼ Fkb T expðβμ∞ AÞ
ð9Þ
So that the HLCS is related to the solvation free energy of transfer of a solute from the ideal gas state to an infinitely dilute
ð10Þ
ÆV exp½βUðR0 ÞæN, P, T ÆV æN, P, T
ð11Þ
In liquids at low pressure, the volume fluctuations are small; thus, the calculation of chemical potential in either the canonical (NVT) ensemble or the isobaricisothermal (NPT) ensemble is essentially the same. The distinction is that it is simpler to use the NVT ensemble, but the experiments correspond to using the NPT ensemble. Equation 10 is essentially the formula for the Widom test particle insertion (TPI) method.32 For simulation of a dense liquid system with a complex solute molecule, random insertion attempts often result in particle overlaps, giving a very large potential energy and little or no contribution to the ensemble average. Therefore, the Widom TPI method is inefficient at high densities. To overcome this difficulty, different simulation techniques have been devised, but they are generally based on two approaches;29 one is computing the difference in the free energy between two systems as in the thermodynamic integration (TI) method,33,34 and the other is from the ratio of the partition functions determined using the conventional Metropolis algorithm,35 such as the free-energy perturbation (FEP)36 or Bennett acceptance ratio (BAR) methods.37 [Note that the Widom TPI method can be considered to be a single-step scheme of the FEP method in the forward direction (insertion).36] An overview and evaluation of these methods for the determination of free energies can be found in a review by Beveridge and Dicapua,38 and a more detailed description of these methods is given by Chipot and Pohorille.39 The accuracy of the TI method is affected by the choice of numerical integration method and the number of intermediate points, and that of the ratio method depends on the acceptance efficiency. A number of ratio methods have been developed to increase the efficiency of acceptance. One is the expanded ensemble (EE) method proposed by Lyubartsev et al.,40 which has been used to calculate the free energy and chemical potential of alkanes, aromatic, and organic compounds in pure liquids and in liquid solutions4146 and to determine the solvation free energies for aqueous ionic solutions.47 The solubility of a series of drug-related compounds has also been studied using the EE method, and the solvation free energies were calculated with a precision of ∼2 kJ/mol.48 The EE method increases the acceptance efficiency by calculating the ratio of partition functions in a 18661
dx.doi.org/10.1021/jp205312k |J. Phys. Chem. C 2011, 115, 18659–18669
The Journal of Physical Chemistry C
ARTICLE
number of subensembles, each of which is associated with a value of a coupling parameter. In this work, we have used the EE method, which is especially suited for calculating the chemical potential of molecular solutes at moderately high solvent densities. We now briefly describe the theoretical background of the EE method. Consider a solution with N solvent molecules into which one solute molecule is gradually inserted. The configurational partition function for a single NPT subensemble is Z
Q ðN, P, T, λi Þ ¼
Z
dV expðβPV Þ
þ λi hNþ1 ðqNþ1 ÞÞÞ
dqN dqNþ1 expðβðHN ðqN Þ
ð12Þ
N
where q are the internal degrees of freedom of the N solvent particles, qN+1 are the internal degrees of freedom of the partially inserted N + l solute particle, HN is the Hamiltonian of the system of N solvent particles including both intramolecular and intermolecular interactions, hN+1 is the intermolecular interaction energy of the solute particle interacting with all of the solvent particles, and λi is an additional coupling parameter, varying from 0 (for the noninteracting or ghost solute molecule) to 1 (the solute molecule fully interacting with the solvent molecules) as i changes from 0 to M (the total number of subensembles). The partition function of the expanded ensemble is given by ΩNPT ¼
M
∑ expðηi ÞQ ðN, P, T, λi Þ i¼0
ð13Þ
where ηi is the balancing factor in the ith subensemble introduced to facilitate transitions between subensembles. The probability of visiting the ith subensemble is pðλi Þ expðηi ÞQ ðN, P, T, λi Þ
ð14Þ
A number of Monte Carlo steps based on the conventional Metropolis rule are used to produce a probability distribution over the subensembles. The probability ratio of the two extremes, i = 0 and M, gives the free-energy difference between the solute in the ideal gas and that in the real solution, which is the excess chemical potential μex ¼ ln
Q ðN, P, T, λi ¼ M Þ pðλi ¼ 0 ¼ 0Þ ¼ ln þ ηi ¼ M ηi ¼ 0 Q ðN, P, T, λi ¼ 0 Þ pðλi ¼ M ¼ 1Þ
ð15Þ The intermediate ensembles are nonphysical states and are not used in the calculation of the free energy or chemical potential but are important because they provide a reversible path to connect the initial and the final states. The balancing factors are adjusted in a series of trial runs to obtain an approximately uniform distribution of subensembles. We adopted the adaptive algorithm of Åberg et al.47 and revised by Chang et al.46 to iteratively adjust the balancing factors such that the acceptance ratios of forward and backward transitions between two adjacent subensembles are equal to each other.
3. SIMULATION DETAILS The MFI-type silicalite crystal structure was used in this study. Silicalite has three crystalline forms,49 a monoclinic form (MONO), an orthorhombic form with Pnma symmetry (ORTHO), and an orthorhombic form with P212121 symmetry (PARA). We performed simulations using the ORTHO form. The silicalite structure was constructed according to International Zeolite
Association (IZA) specifications,50 with lattice parameters of (a) 20.090, (b) 19.738, and (c) 13.142 Å, and has two types of channels, a straight channel parallel to the y direction and zigzag (sinusoidal) channels along the x direction. A single unit cell is made up of 96 silicon and 192 oxygen atoms and contains four straight channel segments, four sinusoidal channels, and four channel intersections. The pore size of each channel is about 5.5 Å. Following others,15,24,26,51 we assumed that the zeolite crystal is rigid in our simulations to increase computational efficiency. We note that the framework flexibility could have impact on the extent of adsorption. In some studies, the influence of framework flexibility has been considered. Nair and Tsapatsis52 found that the adsorption of o-xylene or m-xylene causes distortion of the straight channels of silicalite and suggested that a flexible framework should be considered in modeling. Sorenson et al.53 observed crystal expansion due to alkane adsorption but no significant change with benzene. Vlugt and Schenk54 found the influence of the flexible framework on the adsorption of hydrocarbons to be quite small at low loadings, but for large molecules such as isobutane and heptane, there is an effect at high loading. Demontis and Suffritti55 discussed in detail the effect of the framework flexibility in MD simulation of zeolites and related materials and stated that a flexible framework approach should be compulsory when doing simulations for systems in which the diameters of the largest openings and those of the diffusing molecules are similar. A silicalite structural transition has been found with large aromatic molecules24 and tetrachloroethene,56 and in both cases, a structural transition occurs at higher loadings where more pore space is required. As our study focuses on methanol and ethanol adsorption, which are relatively small molecules, and on the basis of the shape of the experimental isotherms, no structural transition or chemical changes appear to occur in silicalite even up to saturation loadings at room temperature;12 therefore, flexibility of the framework would not be expected to significantly affect amount of adsorption. In this work, GCMC simulations were used to calculate the extent of adsorption, and EE simulations were used to calculate the infinite-dilution excess chemical potential of the solutes in the aqueous solution. Simulations were performed at 298 or 303 K, depending on the temperatures at which the experiments were performed. Each GCMC simulation was performed at fixed fugacity (chemical potential), volume, and temperature using at least 2 2 2 unit cells of silicalite. In the low loading regime, we increased the number of unit cells to have a sufficient number of molecules adsorbed to obtain statistically reliable results. We performed GCMC simulations with four types of moves, (1) molecular translation based on center of mass (COM) positions, (2) molecular rotation based on COM, (3) molecular insertion with random position and random orientation, and (4) deletion of a molecule in the system. Each GCMC move was attempted with a probability ratio for displacement, rotation, insertion, and deletion of 4:2:3:3, respectively. The configurations of the inserted molecules were sampled from a larger pool of pregenerated configurations from ideal gas MD simulations at the same temperature. The energy bias scheme devised by Snurr et al.57 was used to increase the insertion efficiency. To increase the equilibration efficiency, the initial configurations for each run were taken from the equilibrated configurations of the previous run. The simulations were allowed to equilibrate for at least 1 107 Monte Carlo (MC) steps before data production with another 5 107 MC steps to sample the thermodynamic 18662
dx.doi.org/10.1021/jp205312k |J. Phys. Chem. C 2011, 115, 18659–18669
The Journal of Physical Chemistry C
ARTICLE
Table 1. Interaction Parameters for the Potentials Used in This Studya model
interacting site 16
methanol
ethanol16
σ (Å)
ε (kcal/mol)
q (e)
C
3.50
0.066
0.145
O H (CH3)
3.12 2.50
0.170 0.030
0.683 0.040
H (OH)
0.0
0.0
C (CH3)
3.50
0.066
0.418 0.180
C (CH2OH)
3.50
0.066
0.145
O
3.12
0.170
0.683
H (CHx)
2.50
0.030
0.060
H (OH)
0.0
0.0
water66
O H
3.165492 0.0
0.1554253 0.0
0.8476 0.4238
water67
O
3.150574
0.15210
0.834
H
0.0
0.0
O
3.165492
0.1554253
water68
H
0.0
0.0
water69
O
3.165492
0.1554253
H
0.0
0.0
silicalite
O Si
3.20 0.0
0.18678 0.0
0.418
0.417 0.82 0.41 0.82 0.41 1.025 2.050
adsorbed phase, the same OPLS-AA force field16 was used for the alcohols in both the GCMC and EE simulations because in simulations, this force field can reproduce the thermodynamic properties of the pure liquids reasonably well. To test the effect of the water model in the EE free-energy calculations, we used the following three-site models for water: the extended SPC (SPC/E) model,66 the TIP3P model,67 the flexible SPC model (SPC/Ft) of Toukan and Rahman,68 and the flexible SPC model (SPC/Fw) of Wu et al.69 Silicalite was modeled using the force field parameters of Kuhn et al.21 as a starting point. We slightly modified the LJ parameters to match the experimental vapor-phase adsorption loadings of pure methanol (shown later). For the water and silicalite interaction, we used the potential parameters reported in Yazaydin25 without modification. The force field parameters used in this study are listed in Table 1. The geometric combining rules commonly used in the OPLS-AA force field were used for unlike atom pair interactions not mentioned. The intermolecular interactions were truncated at a cutoff distance of 12.0 Å, and longrange corrections were applied. The electrostatic interactions were computed using the Ewald summation method.58 Standard periodic boundary conditions and the minimum image convention were employed in all three dimensions. Unless otherwise noted, the uncertainties in the simulation results are smaller than the symbol sizes in the figures presented in the next section.
a
Note: The intramolecular potential models for methanol, ethanol, and water are not listed but can be found in the original references. The interaction parameters for the unlike pairs are estimated using geometric rules for both the σ and ε parameters, except for water and silicalite LJ 12-6 interaction parameters (σOwOz = 4.24 Å and εOwOz = 0.21858 kcal/mol) that were taken from Yazaydin.25
properties of interest. More detailed information about the GCMC method can be found elsewhere.29,58 Our in-house parallelized GCMC code was used previously for the adsorption simulations of complex molecules in metalorganic frameworks (MOFs).5961 The bulk phase EE simulations were performed in a NPT ensemble at 1 atm of pressure and 298 K temperature. In the EE simulations, MD simulations in the NPT ensemble are used to provide phase space trajectories and the transitions between subensembles that were chosen by the conventional MC Metropolis rule. More detailed information about EE simulation can be found elsewhere.47 We used 512 water molecules and a single solute molecule in each run. This system size has been found to be large enough to eliminate size effects in the calculations. The double-time step method was used with the r-RESPA algorithm of Tuckerman et al. in MD62 to integrate the equations of motion. The short time step of 0.2 fs was used for the intramolecular degrees of freedom, and the long time step of 2.0 fs was used for the other interactions. Fixed bond constraints with the RATTLE algorithm63 were used for water. The NPT ensemble was maintained using the NoseHoover method.64,65 The total number of subensembles used in each simulation was 25 to obtain average transition probabilities of 5060%. Attempts to make transitions among the subensembles were performed every 20 fs. The balancing factors were updated after every attempted transition between subensembles in the first 2 ns. The total simulation time for each run was 12 ns. The nonbonded intramolecular and intermolecular interactions were described by the LJ 126 potential plus a Coulombic point-charge interaction. As one of the goals of our method is force field consistency in simulations of the bulk fluid and the
4. RESULTS AND DISCUSSION In Figure 1, we show the calculated adsorption isotherms of pure methanol and ethanol vapors in silicalite and compare them with the experimental data of Dubinin et al.12 at 303 K. Our GCMC simulations were performed at both 303 and 298 K. The LJ parameters for silicalite were first optimized for the adsorption of hydrocarbons by Pascual et al.70 Recently, Kuhn et al.21 used them to model methanol and ethanol adsorption in all-silica DD3R zeolite, with methanol and ethanol modeled with the TraPPE force field.28 The simulation results qualitatively predicted the trend of experimental isotherms. Gupta et al.15 used the OPLS-AA force field16 to model methanol in silicalite but increased the ε parameter for the zeolite oxygen57 by 20% and varied the partial charges on the atoms of methanol only in its interaction with silicalite for better agreement of the simulated and measured Henry’s constants. (The original OPLS-AA empirical charges were still used for the methanolmethanol interaction.) Except for an overestimation at higher loadings (corresponding to higher pressures) and a slightly underestimated heat of adsorption at zero occupancy, the results predicted by Gupta et al.15 were in reasonably good agreement with experiment. This provides evidence that good predictions can be made using bulk fluid-phase potentials for the adsorbents and published lattice oxygen parameters. This procedure of adjusting only the interaction parameters between the adsorbent and framework to reproduce the experimental data has also been used for water adsorption in silicalite.7173 On the basis of those published results, we used the OPLS-AA force field16 for methanol with the original empirical charges for both methanolmethanol and methanolsilicalite Coulombic interactions. In our simulations, we increased the value of the σ parameter of the zeolite oxygen by about 7% from the value given by Kuhn et al.21 The resulting calculated vapor-phase adsorption isotherms for pure methanol are in good agreement with experimental data, as shown in Figure 1a. At 303 K, the predictions generally match the experimental data very well, though at higher 18663
dx.doi.org/10.1021/jp205312k |J. Phys. Chem. C 2011, 115, 18659–18669
The Journal of Physical Chemistry C
ARTICLE
Figure 1. Calculated pure vapor-phase adsorption isotherms of (a) methanol and (b) ethanol in silicalite compared to the experimental data of Dubinin et al.12 The insets show the low loading regime corresponding to low pressures.
Figure 2. Calculated heats of adsorption for pure methanol and pure ethanol from the vapor phase compared to the experimental results of Dubinin et al.12
loadings, corresponding to higher external pressures, there is a slight underestimation of the loading. At low loadings (
