2836
J. Phys. Chem. B 2007, 111, 2836-2844
Diffusion of Methanol in Zeolite NaY: A Molecular Dynamics Study David F. Plant,†,‡ Guillaume Maurin,‡ and Robert G. Bell*,† DaVy Faraday Research Laboratory, Royal Institution of Great Britain, 21 Albemarle Street, London W1S 4BS, United Kingdom, and Laboratoire LPMC, UMR CNRS 5617, UniVersite´ Montpellier II, Place E. Bataillon, 34095 Montpellier cedex 05, France ReceiVed: NoVember 10, 2006
Molecular dynamics simulations were performed in order to obtain a detailed understanding of the selfdiffusion mechanisms of methanol in the zeolite NaY system. We derived a new force-field term to describe the interactions between the methanol molecules and the extraframework cations. From the simulations, we show that diffusive behavior in the high-temperature range consists of a combination of both short- and long-range motions at low and intermediate loadings. This type of motion is characterized by an activation energy that decreases as the loading increases. At low loadings, we also observe short-range diffusive behavior based on a surface-mediated mechanism. The short-range behavior corresponds to motion only on the length scale of an FAU supercage, whereas the long-range behavior involves intercage diffusion. For the saturation loading corresponding to 96 methanol molecules per unit cell, only short-range motions within the same supercage predominate. Finally, the preferential arrangement of the adsorbate molecules around the extraframework cations are examined and compared with those previously deduced from experimental data.
1. Introduction Novel applications of methanol include DMFC fuel cells, which use a methanol/water solution as fuel, a clean technology that might be used to power vehicles in the future.1 Methanol is also involved in numerous catalytic processes. It can be produced via Fischer-Tropsch reaction from synthesis gas (carbon monoxide and hydrogen), instead of from oil,2 and its further transformation to hydrocarbons up to C10 is the basis of several industrially important reactions such as Mobil’s methanolto-gasoline (MTG) and methanol-to-olefin (MTO) processes.3,4 The key steps of these reactions rely on the adsorption of methanol at the surface of the catalyst and the formation of the first carbon-carbon bond. Zeolite materials play a crucial catalytic role as membrane reactors in this process either by interacting with methanol molecules to form an intermediate species during the chemical reaction5 or by increasing the rate of conversion via high separation selectivity of the products.1,6 Recently, a mechanism involving protonated formaldehyde and methane intermediates was proposed to explain this process in chabazite-type zeolites.7 Many reaction schemes and mechanisms have also been proposed to account for the catalytic behavior of HZSM-5 zeolite in such processes.8 Nevertheless, a complete understanding remains elusive. A great number of theoretical and experimental studies have focused specifically on the first stage of the MTG process, corresponding to the adsorption of methanol at a Brønsted acid site in the catalyst. The question of physisorption or chemisorption of methanol at this specific site has been the subject of much debate. Experimental results from infrared spectroscopy9 and nuclear magnetic resonance spectroscopy10 were interpreted in terms of proton transfer occurring during adsorption, from the bridging hydroxyl group of the zeolite to the adsorbate. Ab initio calculations performed by Shah et al.11 demonstrated this transfer * Corresponding author. E-mail:
[email protected]. † Royal Institution of Great Britain. ‡ Universite ´ Montpellier II.
to be unstable and suggested a mechanism in which the adsorbed methanol is physisorbed on the catalyst with the formation of hydrogen bonds with the Brønsted acid sites. Several later ab initio and DFT calculations have supported these findings.12,13 Methanol can also be used as an alkylating agent for aromatic compounds in zeolite-catalyzed reactions.14 The acido-basic nature of zeolites that can be modified using ion exchange, various substitutions, the introduction of dopants, and/or dealumination, allows them to act as efficient alkylating catalysts. HZSM-5 can thus be used selectively to alkylate toluene with methanol creating p-xylene, and much research effort has been focused on the modification of the ZSM-5 structure and acidity to improve product selectivity. In contrast, under basic conditions, side-chain (as opposed to ring) alkylation is favored.15 This reaction takes place in a number of basic zeolites, including alkali-metal-exchanged X and Y zeolites.16,17 In a recent study, Palomares et al.,17 carrying out in situ infrared spectroscopy experiments, investigated the selective alkylation of toluene in basic zeolites as a function of the alkali cation species by studying the surface chemistry of the reactants under working conditions. The optimization of each of the processes described above demands a sophisticated understanding of the interactions between the zeolite surfaces and the reactant methanol molecules, including the transport properties of the reagent within the catalyst micropores. Only a limited number of articles are available in the literature that deal with the interaction between methanol and basic zeolite materials. Infrared experimental studies performed by Rep et al.18 proposed various adsorbed methanol geometries in which methanol is coordinated to the extraframework sodium cations, depending on the structure (ZSM-5, MOR, and FAU) and composition of the zeolite framework (NaY and NaX). Different modes of adsorption were identified, with or without additional physical interactions between the hydroxyl hydrogen and the zeolite, and a further case where intermolecular hydrogen bonding occurs between
10.1021/jp0674524 CCC: $37.00 © 2007 American Chemical Society Published on Web 02/27/2007
Diffusion of Methanol in Zeolite NaY the adsorbed molecules, giving rise to methanol clusters. These experimental results were successfully compared with those obtained from DFT calculations.19 Experimental results obtained on Na-ZSM5, Na-MOR, or Na-FAU by temperature-programmed desorption (TPD),20 calorimetry,9 inelastic neutron scattering (INS),21 and Fourier transform infrared (FTIR) spectroscopy22 clearly show that the methanol molecules are not dissociated, even at high temperatures up to 623 K, and that they strongly interact via their oxygen atoms with the extraframework cations. Furthermore, it was reported that the adsorbate molecules form clusters characterized by either weak or strong hydrogen bonding. This latter finding was supported by a statistical thermodynamics study on Na-ZSM5 that showed that the pore filling proceeds via the formation of monomeric to tetrameric methanol/cation complexes located in the channel intersections.23 Here, for the first time, we report an investigation of both the microscopic self-diffusion mechanism and the preferential arrangement of methanol in zeolite NaY by means of molecular dynamics simulations. We use classical simulations based on interatomic potentials, as the adsorption processes are known to be entirely physical in nature under ambient conditions. Furthermore, we were able, in part, to use potential models developed by others that proved to be well-suited to our simulation systems, but we needed to characterize the interactions between the extraframework cations and the methanol molecules, as these had not previously been considered in such a simulation. In the first step of this work, we therefore derived a set of potential parameters, consistent with the other terms in our model, using ab initio cluster calculations to interrogate the interactions between the sodium cations and the oxygen of the methanol. The new interatomic potential is simple enough to be transferable between different cation-containing zeolite systems. We selected pair-potential parameters from the cvff force field24 to describe the interactions between the methanol molecules, which were then adjusted to reproduce the diffusion coefficient of methanol in the liquid phase. Additional potentials were used to represent the interactions between the methanol molecules and the zeolite and within the zeolite framework. Using these interatomic potentials, both the diffusion coefficients for methanol in a wide range of temperature (500-700 K) and the activation energies corresponding to the adsorbate motions within faujasite were evaluated as a function of the loading. These simulated data were then compared and contrasted with those obtained experimentally by Grenier et al.25 using pulsedfield gradient (PFG) NMR spectroscopy. In addition, we studied the preferential adsorption sites of methanol and compared them with the previously proposed sites. Overall, comparisons between simulation and experiment also allowed us to establish the validity of the force field, with its newly derived terms. 2. Computational Methodology A. Interatomic Potentials. I. Zeolite Framework. The successful simulation of methanol self-diffusion in the NaY system requires an accurate description of the potential between the methanol molecules and the zeolite framework, including the sodium cations, and among the methanol molecules. Furthermore, we needed a fully flexible model in which both the zeolite and the sorbate molecules were free to alter their internal configurations. The force field used to describe the zeolite framework included the partially ionic model of Ramsahye and Bell26 in which the atoms carry the following partial charges (in electron units): Si, +2.4; Al, +1.4; Oz, -1.2; and Na, +1.0. These values are consistent with those previously used by Kramer et al.27 The short-range interactions were
J. Phys. Chem. B, Vol. 111, No. 11, 2007 2837 TABLE 1: Pair-Potential Parameters and Charges Carried by the Atoms Used to Describe the Interactions within the Framework and between the Framework and the Extraframework Cations atom charges species
charge (e)
Al Si Na Oz
1.40 2.40 1.00 -1.20 Buckingham potential
ion pair
A (eV)
F (Å)
C (eV Å6)
Si-Oz Al-Oz Na-Oz Oz-Oz
30023.0 26998.0 8200.0 894.6
0.1621 0.1622 0.2180 0.3244
12.84 12.84 11.80 0.00
three-body angle potential angle
k (eV rad-2)
θ0 (deg)
Oz-Si-Oz Oz-Al-Oz
12.10 2.20
109.47 109.47
described by Buckingham potentials, including explicit Si-O and Al-O terms, with the corresponding parameters reported in Table 1.
V(r) ) A exp(-r/F) - C/r6 In addition, harmonic three-body terms were defined for the O-Si-O and O-Al-O intratetrahedral angles.
1 V(θ) ) kθ(θ - θ0)2 2 II. Intermolecular Terms. The polarizabilities of silicon and aluminum atoms, which are much lower than those of oxygen atoms, suggest that the repulsion dispersion contribution of the zeolite framework can be assigned only to oxygens of the framework (Oz) and extraframework cations (Na+). Intermolecular potentials between the zeolite framework and methanol molecules were based on those reported by Blanco and Auerbach.28 In these potentials, the short-range parts are given in the form of the Lennard-Jones 12-6 function.
V(r) ) 4
[(σr ) - (σr ) ] 12
6
Although the methanol partial charge distribution adopted by Blanco and Auerbach accurately reproduces the dipole moment, we used different values for the atoms of the methyl group, making the C-H bonds slightly more polarized (see Table 3 below), which is consistent with the ratio of charges we obtained from a Mulliken analysis of DFT calculations, which are described more fully in a later section. For this reason, and also because of the greater ionicity of the NaY framework compared to the siliceous models (DAY and silicalite) investigated by Blanco and Auerbach,28 the methanol-framework potential parameters were scaled slightly. Furthermore, the potential between the oxygen of methanol (Om) and the oxygen of the zeolite framework (Oz) was found to be too attractive in our model and was corrected by scaling the corresponding parameter (see Table 2). The Lennard-Jones parameters that represent the intermolecular interactions between the methanol molecules were similarly slightly adjusted from those reported in the cvff force field.24
2838 J. Phys. Chem. B, Vol. 111, No. 11, 2007
Plant et al.
TABLE 2: Parameters to Describe the Interactions between the Framework and the Extraframework Cations and among the Methanol Molecules Lennard-Jones potential ion pair
ij (eV)
σij (Å)
Om-Oz Hm-Oz C-Oz H-Oz Na-Hm Na-H Na-C
0.00762 0.00512 0.00591 0.00387 0.00878 0.00878 0.01340
4.104 3.374 4.310 3.732 2.6285 2.6285 3.0722
Buckingham potential ion pair
A (eV)
F (Å)
C (eV Å6)
Na-Om
4521.60
0.2177
0.00
Because of the presence of extraframework sodium cations, the simulations required an additional contribution able to reproduce the interactions between the sodium ions and the methanol. The interactions between the sodium ions and both the methyl group and the hydrogen of the hydroxyl group were taken from the cvff force field.24 A new potential between the extraframework cations and the oxygen of the hydroxyl group was then derived from ab initio calculations that adopted the following procedure: A potential energy curve was first generated from single-point energy calculations performed on a gas-phase cluster of one methanol molecule and one sodium ion (see Figure 1). We know that such a model introduces an overestimation of the electric field generated by the bare sodium ion compared to the zeolite system. However, the main features of the adsorbate/cation interaction, such as the equilibrium distance, are reproduced quite satisfactorily. Hence, this method constitutes a first approximation to develop a suitable interatomic
Figure 1. Schematic representation of the methanol-Na+ geometry considered for the calculation of the potential energy curve.
TABLE 3: Pair-Potential Parameters and Charges for Methanol Molecules atom charges species
charge (e)
C H Om Hm
-0.093 0.100 -0.432 0.225 Lennard-Jones potential
ion pair
ij (eV)
σij (Å)
H-H C-H Om-H C-C Om-C Om-Om
0.00165 0.00338 0.00404 0.00694 0.00828 0.00988
2.450 2.920 2.650 3.475 3.150 2.860
two-body bond potential ion pair
k (eV Å-2)
r0 (Å)
C-H C-Om H-Om
29.56 33.33 46.97
1.105 1.420 0.945
three-body angle potential angle
k (eV rad-2)
θ0 (deg)
C-Om-Hm H-C-Om H-C-H
5.60 5.50 4.40
108.32 106.90 108.38
four-body torsonial potential dihedral angle
k (eV)
R
β
H-C-Om-Hm
0.00762
1.0
3.0
potential, which was subsequently refined by comparison with experimental data. A preliminary DFT calculation based on the PW91 functional29 at the double-numeric-polarized (DNP) level of theory by means of Dmol3 30 was carried out in order to define the optimized geometry of this cluster and thus to provide a suitable starting configuration. It was found that the equilibrium geometry of this cluster is characterized by a Na+-Om distance of 2.4 Å and a C-Om-Na+ angle of 117°. This geometry is in very good agreement with the results of previous experimental and theoretical investigations on the solvation of sodium by methanol.31 The Na+ cation was then displaced from the methanol oxygen at distances of between 1.0 and 8.0 Å, in increments of 0.1 Å, in the plane defined by C-Om-H. The geometries were constrained to maintain the previously defined equilibrium C-Om-Na+ angle. At each increment, a singlepoint energy calculation was performed. These calculations were conducted using Gaussian 9832 at the MP2 level of theory with the 6-31G(d,p) basis set.33 The Z matrix forming the geometry of the methanol cluster was created using optimized bond lengths and angles provided by the preliminary DFT cluster calculation described earlier, which has also been discussed elsewhere.34 An additional single-point energy calculation was performed along the vector at a large distance of 100 Å. The resultant energy value was then used as a “zero” point to calculate interaction energies for each Na+-Om distance to produce the energy profile. The partial charges carried by each atom of the methanol molecule were based on a Mulliken charge distribution analysis performed after a DFT calculation on the cluster formed by the adsorbate molecules and the cation located in site SII. These charges are reported in Table 3. The potential energy curve was then fitted by the combination of a Coulombic
Diffusion of Methanol in Zeolite NaY contribution and a Buckingham potential. In fact, as mentioned above, the MP2 calculations would be expected to overestimate the strength of interaction between methanol and an interstitial cation, because the zeolite framework will screen the charge of the sodium, attenuating its effect compared to that of a gasphase cation. The fitted potential was therefore scaled on the energy axis such that the heat of adsorption of methanol on NaY at the low limit of coverage and 300 K, as calculated using the sorption module of Cerius2,38 was 80 kJ mol-1, as found experimentally by Izmailova et al.9 The values of the A, F, and C parameters are included in Table 2. The heat of adsorption calculation also served to validate the other zeolite-methanol intermolecular terms described above. III. Methanol Intramolecular Terms. Intrasorbate flexibility was described by two-body bond stretches, three-body angle bends, and four-body torsional potentials.
1 1 V(r) ) kr(r - r0)2 + kθ(θ - θ0)2 + kφ[1 + Rcos(βφ)] 2 2 The corresponding parameters were obtained by slightly modifying those given by the cvff force field.24 The validity of the modified potential parameters used for describing both the interactions between methanol-methanol pairs and the guest molecule flexibility was checked by comparing the simulated diffusion coefficient of methanol in the liquid phase at ambient temperature with the value obtained experimentally. All intraand interatomic potential parameters are reported in Table 3. B. Zeolite and Liquid Models. The NaY model had the chemical composition Si136Al56Na56O384. The zeolite framework was built in accordance with Lo¨wenstein’s Al-O-Al avoidance rule.35 The second step consisted of modeling the location of the extraframework cations among the different crystallographic sites. This distribution was defined as follows, based on the structure refined from X-ray diffraction data by Fitch et al.:36 6 cations in sites SI located in the center of the hexagonal prism connecting two sodalite cages, 18 in sites SI′ in the sodalite cage in front of the 6-ring window connected to the hexagonal prism, and 32 in sites SII in 6-ring windows connecting the supercages and sodalite cages. The sites occupied were chosen randomly, with the exception that neighboring SI and SI′ sites were not allowed to be both occupied. This structure was then energy-minimized with the GULP code,37 using the previously described interatomic potentials developed by Ramsahye and Bell26 with the constraint that the cell remained cubic. It has already been noted26 that this minimization results in a very slight change in the lattice parameter compared to experiment. This optimized structure was then loaded with 8, 16, 32, 48, 64, and 96 methanol molecules using the sorption module in the Cerius2 program.38 The loading of 96 methanol molecules per unit cell corresponds to the maximum adsorption capacity measured by microcalorimetry at ambient temperature.39 Prior to the molecular dynamics simulations, all generated structures were optimized using GULP in order to provide lowest-energy starting configurations. The liquid methanol system was simulated using 205 methanol molecules distributed within a cubic simulation box of cell length 24.2068 Å, which was the cell parameter of our model dealuminated zeolite Y structure. This is equivalent to the density of methanol at standard temperature and pressure.40 C. Molecular Dynamics Simulations. The interatomic potentials previously described for modeling the whole system were implemented in the DLPOLY program41 in the NVT ensemble using the Evans isokinetic thermostat.42 We selected the optimized structures obtained by the minimization procedure
J. Phys. Chem. B, Vol. 111, No. 11, 2007 2839 as starting configurations, and the minimized cell dimensions were kept fixed during the molecular dynamics (MD) runs. All components of the system (adsorbate and adsorbent) were then treated as fully flexible during the MD simulations. A time step of 1 fs was selected, with simulations run at loadings of 8, 16, 32, 48, 64, and 96 molecules per unit cell, or an average of 1, 2, 4, 6, 8, and 12 molecules per supercage, respectively. The simulations spanned a range of temperatures between 300 and 700 K, each for 106 steps (i.e., 1 ns), following 50000 steps of equilibration. For the cases of 8 and 16 molecules per unit cell, longer runs of 5 ns were additionally carried out at temperatures between 300 and 450 K. A short-range cutoff of 8.50 Å was used, and electrostatic interactions were evaluated using the Ewald method. The trajectory was recorded every 200 steps during the production stage, and radial distribution functions were recorded every 500 steps. The mean square displacements (MSDs) of the methanol molecules at each loading and at the different temperatures were evaluated by means of the following classical equation:
MSD(t) ) 〈∆rj2(t)〉 )
1
N
1
N
∆rj2(t) ) ∑[rj(t) - rj(0)]2 ∑ N j)1 N j)1
where N corresponds to the number of methanol molecules considered in the computation of the MSD and we used multiple time origins in order to improve the statistics of the calculation. When the MSD corresponded at least to a significant distance (100-150 Å2), the diffusion coefficients evaluated for each loading were obtained by fitting the MSD plots as a function of the time in the region 0-200 ps, assuming the following Einstein relation:
MSD(t) ) A + 6Dt The activation energies corresponding to the self-diffusion processes were then evaluated from the linear least-squares fit to the Arrhenius plots of ln(D) ) f(1000/T). 3. Results and Discussion Figure 2 reports the MSDs for methanol in NaY at the various loadings investigated (8, 16, 32, 48, 64, 96 methanol molecules/ unit cell) as a function of the temperature. At all of these loadings, we observed that the MSDs increase with temperature. Furthermore, for the lower loadings (up to 48 methanol molecules/unit cell), the MSDs appear to be reasonably linear over a broad time domain, indicating unrestricted threedimensional diffusion within the zeolite pore system. The values of the MSDssup to 900 Å2 at 400 ps for higher temperaturess denote unhindered diffusion across multiple supercages. In some cases, it can also be seen that diffusivity does not always increase smoothly with temperature, as, for instance, is evident from the noticeable “gap” between the plots for 500 and 600 K in Figure 2b. Similar behavior can be discerned in the MSD plots for 32 and 48 methanol molecules/unit cell (Figure 2c,d). For higher loadings (64 and 96 methanol molecules/unit cell), the MSD plots (Figure 2e,f) appear less linear, and the methanol trajectories are more restricted (MSDs reach less than 90 Å2 at high temperatures), which suggests displacements limited to the length scale of one supercage, the diameter of which is approximately 12 Å. Similarly, at lower loadings and temperatures, the average displacement of the molecules is less than the dimensions of a supercage. For this reason, we carried out additional longer MD runs of 5 ns at the lower temperatures (300-450 K) for both 8 and 16 molecules/unit cell. The aim was to determine whether long-range motion simply occurs on
2840 J. Phys. Chem. B, Vol. 111, No. 11, 2007
Plant et al.
Figure 2. MSD plots for methanol in NaY at various temperatures and at loadings of (a) 8, (b) 16, (c) 32, (d) 48, (e) 64, and (f) 96 methanol molecules per unit cell.
Figure 3. MSD plots for methanol in NaY at temperatures of between 300 and 450 K and at loadings of (a) 8 and (b) 16 methanol molecules per unit cell.
a slower time scale at lower temperature or whether motion is restricted to a shorter length scale, as seen at higher loadings where the MSDs flatten. The resulting MSD plots are shown in Figure 3. For 8 molecules (Figure 3a), there is still little evidence of long-range displacement at 300 K. At 350 and 400 K, average motion remains within the length scale of a supercage with some evidence of leveling off of the plots, indicative of restricted motion. Similar behavior is observed for 16 molecules/ unit cell except that, at 350 K and above, the magnitude of the MSD suggests motion of methanol across more than one supercage. For the purposes of this discussion, therefore, we define two types of motion, “short-range” and “long-range”. Short-range motion is on the length scale of a single FAU supercage
(intracage), whereas long-range motion is not as restricted and involves diffusion throughout the pore structure combining both intra- and intersupercage motions. The evidence (see also below) suggests that, at lower loadings and temperatures (roughly
