Molecular Dynamics Simulations of Methanol to Olefin Reactions in

Article pubs.acs.org/JPCC

Molecular Dynamics Simulations of Methanol to Olefin Reactions in HZSM-5 Zeolite Using a ReaxFF Force Field Chen Bai, Lianchi Liu, and Huai Sun* School of Chemistry and Chemical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China S Supporting Information *

ABSTRACT: A ReaxFF force field has been extended and used in the present work to simulate methanol to olefin (MTO) reactions in H-ZSM-5 zeolite. By explicitly considering multibody interactions and thermodynamic conditions, the initial reaction network of MTO in acidic zeolite has been obtained. New reaction mechanisms are proposed based on the simulations. For the activation of methanol, a less possible but very important CH3 radical mechanism is identified in addition to the commonly accepted methoxyl mechanism. The commonly accepted chain-growth mechanism, in which ethene interacts with methyloxide, has been observed. However, it is a small contribution to the total production. The more popular route for the chain growth is attributed to the presence of deprotonated Brønsted sites, which are produced via the activation of methanol molecules. Therefore, the hydrocarbon pool is working with the methanol molecules involved. With the hydrocarbon pool, the chain growth is significantly accelerated. Considering the collision probability, the rate-determining step for MTO is not the activation of methanol as suggested by static calculations but the C−C chain formation. The simulation data are consistent with previously reported experimental observations.

1. INTRODUCTION The methanol to olefin (MTO) reaction, catalyzed by ZSM-5 or SAPO-34, is one of the most promising alternative energy technologies. Methanol can be economically obtained from various sources, such as natural gas, biomass, or even organic wastes. Hence, the MTO process has drawn considerable attention from both scientific and industrial communities. However, the reactions are highly complicated, and the reaction mechanisms of MTO are far from being fully understood.1 A previously proposed direct mechanism,2,3 in which the C− C bond formation is only from C1 species, has been challenged by recent experimental and computational works. Lesthaeghe et al.4 have summarized the results of direct mechanisms based on theoretical studies. All direct routes involve very high energy barriers or unstable intermediates. On the other hand, higherorder hydrocarbons are formed with only methanol as the feedstock. 5 Some researchers have suggested that the production may be the result of impurities in the reactant or the catalysts,6 but others have raised different opinions.7 Most researchers nowadays believe that the MTO process likely follows the “hydrocarbon pool” (HP) mechanism.8−10 In this mechanism, the short-chain alkenes or polymethylbenzenes generated during the induction time serve as “seeds” for further chain growth. When the chain lengths reach a certain size depending on the pore sizes and topologies, light olefins such as ethene and propene are produced by cracking the intermediates. Two major HP mechanisms are proposed, i.e., the aromatic pathway or alkene pathway.11−14 However, detailed information on HP formation and functions is still lacking.15,16 © 2012 American Chemical Society

All theoretical studies carried out so far have been based on static energy calculations. Such calculations can be performed only after the mechanisms are proposed. Therefore, these calculations are useful in comparing the proposed mechanisms but cannot be used to predict or discover mechanisms. Chemical reactions are controlled by free energies, and the entropies that are ignored in the static calculations are often important contributions to the free energies. In fact, the entropy contribution may change the order of energy differences and the preferred mechanisms.17−19 To understand how chemical reactions occur in reality, reactions in condensed phases should be studied so that entropy contributions can be included. One of the possible approaches to study chemical reactions in condensed phases involves ab initio molecular dynamics (AIMD).20 Given the extremely large computational requirements, the systems that can be treated using AIMD are small.20−23 The recently developed reactive force field (ReaxFF) provides an alternative approach.24 Molecular dynamics simulations can be applied to simulate chemical reactions in a fairly large scale (up to 106 atoms) because energy hypersurfaces are represented using atom−atom interactions.25 Although ReaxFF has been applied to various studies,19,24−28 its full potential is yet to be explored. One of the major hurdles that must be overcome is the accuracy of the underlying force field parameters. Another major difficulty in using ReaxFF Received: January 8, 2012 Revised: March 1, 2012 Published: March 5, 2012 7029

dx.doi.org/10.1021/jp300221j | J. Phys. Chem. C 2012, 116, 7029−7039

The Journal of Physical Chemistry C

Article

simulations in general is time scale limitation. Unless the reactions are extremely rapid (e.g., combustion or detonation), most chemical reactions take place in a much longer time scale than that which the simulations can reach (usually nanoseconds). A common practice is to raise the temperature of the simulation such that reactions can be observed in the nanosecond time scale. However, given that the entropy contribution is significantly sensitive to the temperature, this approach significantly changes the reaction conditions and undermines the conclusions drawn from the simulations. The time-scale problem is essentially a rare-event issue; i.e., the possibility of crossing the transition state is too low to be seen within the simulation time. For simulations with classical force fields, several algorithms such as accelerated molecular dynamics29 and replica exchange molecular dynamics30 have been proposed to enhance the sampling processes of rare events. However, these methods need to be developed to study reaction mechanisms using reactive force fields. In the present work, the existing ReaxFF24,31 parameters for Si, O, Al, H, and C were extended to enable an accurate representation of acidic zeolites. To address the time-scale issue, a simulated annealing scheme was applied to shorten the reaction time. The simulations were mostly performed at a low (close to experimental) temperature but periodically run at a high temperature for a short period to accelerate the reactions. Although this approach is not thermodynamically rigorous, the ample data obtained in a relatively short simulation period provide sufficient information for studying the reaction mechanisms, which hopefully contributes to the understanding of the mechanisms of MTO chemical reactions.

2. MODELS AND METHODS The zeolite was represented by a three-dimensional periodic 2 × 1 × 1 super cell (Figure 1) of the ZMS-5 unit cell. The unit cell parameters and atomic coordinates were taken from experimental data.32 The simulation box had 192 tetrahedral (192T) sites. It includes intersections of straight and zigzag channels. The silicon atoms at the T12 and T9 sites33 were replaced by aluminum, resulting in a Si/Al ratio of 176/16. Three types of models with different chemical components were constructed for the simulations: (A) Acidic zeolite with 15 methanol molecules. (B) Methoxylated zeolite with 10 ethene molecules. (C) Acidic zeolite with 15 methanol and 6 ethene molecules. Model A was designed to investigate the activation of methanol; model B was used to study the initial C−C bond formation; and model C simulated a system with ethene as the HP element. The concentrations of the reactants were three times to six times higher than the normal concentrations in the industrial process. NVT molecular dynamics simulations were performed using the ReaxFF package in the LAMMPS software suite (12 June 2011 version).34 The time step of integration was 0.25 fs, and the total simulation time is 1 ns. The simulation and data collection temperature was 700 K, which was close to the experimental temperature. To speed up reactions, a high temperature (1600 K) was applied for a short period in a fixed interval of simulation time. The temperature control was as follows: running at 700 K for 100 ps, heating up to 1600 K in 10 ps, running at 1600 K for 20 ps, cooling to 700 K in 10 ps, and repeating the above steps. The 20 ps simulation time at 1600 K was estimated such that about 1% to 2% of the slowest

Figure 1. Plane (a) and vertical (b) views of a typical system, a 192T model with a Si/Al ratio of 176/16.

reactions occur based on estimated activation energies, which are discussed in the following section of the current paper. The temperatures were controlled using a Berendsen thermostat with a damping constant of 0.1 ps. To maintain the framework stability, atoms of the zeolite framework were fixed, except for those in the active site fragments (−SiOAl(OH)2O(H)Si−). The cutoff distance for energy evaluation was 10 Å. The QEq method was used to calculate atomic charges at every step.35 The trajectories collected at 700 K were scanned to collect the molecular species based on the bond order information. Unit reactions were identified by tracking the species in the trajectories from the selected intermediates or products. To enhance the statistical reliability of the collected data, three independent runs were performed for models A and B, whereas 15 independent runs were performed for model C. The initial force field parameters for Si, O, C, H, and Al were obtained from the literature.24,31 However, the original parameters for aluminum were derived for aluminum oxide, in which the interaction with silicon is missing. New parameters of aluminum were derived by fitting quantum mechanics density functional theory (DFT) data calculated on model molecules representing an acidic zeolite. The DFT calculations were carried out using the B3LYP functional and the 631g(d,p) basis set. Parameterization was done using the successive one-parameter search technique as described by 7030



Article

Figure 2. Comparisons of energy differences between the force field and QM values of the typical bonds and angle in clusters. (a) Dissociation of the Al−O bond in the Al(OH)3 cluster, (b) energy of H2Al−O−H as a function of the Al−O−H valence angle, (c) energy of Al(OH)3 as a function of the O−Al−O valence angle, and (d) energy of H3SiO(H)AlH3 as a function of the Si−O−Al valence angle.

van Duin et al.36 Quantum DFT calculations were carried out using the Gaussian 03 package.37

deviation between the force field and DFT results is 4.07 kcal/ mol. The maximum difference between the calculated energy barriers of the force field and the DFT results is 4.5 kcal/mol. The order of barrier heights is correctly described by the force field. The force field against the QM/MM results26 calculated using a 128T zeolite cluster model was further examined. The comparison results are given in Table 1. The standard deviation between the force field and the QM/MM results is 4.4 kcal/ mol, consistent with the quality of comparison using the cluster model. The actual energy values are different because of the environmental changes, and the species are stabilized by the framework as indicated by previous calculations.39−41 The order of energy barrier heights for the small and large cluster calculations is the same. 3.2. Model A. The population of the major species obtained from the 1 ns simulation averaged over three independent runs of Model A is shown in Figure 4. Data are collected only for the trajectories at a low temperature (700 K), and the temperature curve is plotted in the figure. In addition to methanol, the major species identified in the trajectories are water, hydrocarbons, or oxides of hydrocarbons with variable numbers (1 to 4) of carbon atoms, CxHy(O), and H2 molecules. Methanol is rapidly decomposed, and 14 methanol molecules are consumed at approximately 250 ps. Water molecules are immediately generated as methanol is decomposed. The population is

3. RESULTS AND DISCUSSION 3.1. Parameterization. The model molecules used for parametrization are given in Figure 2. The figure also shows comparisons between the force field and DFT-calculated energy results for Al−O bond stretching, Al−O−H angle, O−Al−O angle, and Si−O−Al angle bendings. A standard deviation of 6.1 kcal/mol is obtained for these comparisons. For the bond stretching energies, the DFT calculations are conducted with singlet and triplet electronic states. The singlet state is lower in energy at a short separation between Al and O atoms, but the triplet energy becomes lower at a larger separation. The force field results well agree with the DFT data in different regions. The Al−O−H, O−Al−O, and Si−O−Al angles are sampled using H2AlOH, Al(OH)3, and H3SiO(H)AlH3 molecules, respectively. The comparisons given in Figure 2b−d show that the force field represents the energy curves well near the minimum, but the deviations become larger with increased angles. Chemical species previously proposed in reaction mechanisms38 were included in the training set. Figure 3 shows a comparison of the energy diagrams calculated using the force field and using the DFT. The calculations were performed on 3T clusters (HO)3SiOAl(OH)2O(H)Si(OH)3. The standard 7031



Article

Figure 3. Comparison of energy differences between the FF and QM results of the three reaction pathways in the 3T model. (a) First step of the mechanism: decomposition of methanol and formation of surface methoxyl group. (b) Second step of the mechanism: carbon chain growth and formation of propene molecules. (c) Concerted mechanism: carbon chain growth while water molecules are dissociated.

from the beginning; the population monotonically increases and approaches a plateau near the end. The formation of hydrogen molecules is associated with the production of heavier hydrocarbons. On the basis of Figure 4, the reactions can be roughly described as follows. Methanol decomposition produces water, and C1 species as commonly indicated in the literature.7,38,42 The C1 species are very active, and they quickly react with other molecules or species in the system. The C2 species may be directly formed from the C1 species, although the probability is low. As the carbon chain increases, hydrogen

stable as methanol disappears at approximately 250 ps. This finding indicates that the generation of water is directly associated with the decomposition of methanol. The C1 hydrocarbons immediately appear as the reactions start. However, the population is flattened and eventually slightly decreases at ca. 500 ps. Higher-order hydrocarbons and their oxides (C2 to C4) are generated in later stages. Specifically, two C2H4 molecules are observed in the simulations, indicating the possibility of a direct mechanism of C−C bond formation, but the probability is low. The hydrogen molecules are produced 7032



Article

Table 1. Energy Differences between the Force Field and QM/MM Results of Two Selected Pathways of Propene Formation from Ethenea 1st step of SW mechanism

a

2nd step of SW mechanism

concerted mechanism

structure

FF

QM/MM

structure

FF

QM/MM

structure

FF

QM/MM

CH4O + H−Z Ads1_SW1 Ads2_SW1 TS_SW1 Int_SW1 H2O + CH3−Z Ea1

0 −28.5 −25.6 12.1 −6.2 8.5 37.7

0 −26.3 −24.2 17 −7.8 2.9 41.2

C2H4 + CH3−Z Ads_SW2 TS1_SW2 Int_SW2 TS2_SW2 Prod_SW2 propene + H−Z Ea2 Ea3

0 −12.3 15.6 −34.3 −12.5 −35.8 −12.7 27.9 21.8

0 −10.5 10.9 −36.8 −10.3 −32.5 −13 21.4 26.5

CH4O + C2H4 + H−Z Ads_Con TS1_Con Int_Con H2O + C3H7−Z TS2_Con Prod_Con H2O + C3H6 + H−Z Ea4 Ea5

0 −42.4 −10.6 −49 −37.4 −15.6 −27.2 −16 31.8 21.8

0 −36.6 −7.6 −47.5 −33 −7.6 −33.7 −15.1 29 25.4

Both are in the 128T cluster model.

adsorbed on a Brønsted acid site (denoted by H−Z in the figure) to form a van der Waals (vdW) complex, CH3OH2Z. From the complex, two major routes are identified. First is the commonly referred route, in which one H2O molecule is released and the methoxyl group CH3−Z on the surface is formed. From the CH3−Z groups, various reactions are observed. The most common ones are interactions with other species to form new C−C bonds (chain growth). The second route, which has not yet been reported in the literature to the best of our knowledge, involves methanol taking the hydrogen atom from the Brønsted acid site as well as forming a protonated methanol CH3OH2 and a deprotonated Brønsted (basic) site. This reaction occurs only with a collision on the adsorbed CH3OH2Z by another molecule in the system, as shown by an image in Figure 6. This is a dynamic phenomenon that has not been perceived in static calculations. The free CH3OH2 is highly unstable. It quickly dissociates to a water molecule and a free CH3 radical. The CH3 radical only exists in the high-temperature pulse, and its average lifetime is 14 ps. The existence of CH3 and the deprotonated Brønsted (basic) site is important because they are very active, and many species are generated due to the existence of CH3. Several types of reactions are observed for the CH3 radical. It may take a hydrogen atom from the acidic site to form CH4. Most popularly, it interacts with another species to form new C−C bonds (chain growth). The two observed direct formations of C2 are both from the CH3 radical. As the snapshots in Figure 7 (Ia−Ic) show, a free CH3 group attacks a methanol adsorbed on the Brønsted site. Subsequently, the C−O bond breaks, and the C−C bond forms, producing ethane and water molecules. The C−C bond may be formed between the radical and methoxyl group as shown in Figures 7IIa−IIc. The surface methoxyl is attacked by a free CH3 group to form a carbene-like intermediate. The intermediate returns a hydrogen to the framework and produces a C2H5 radical. The C2H5 radical can easily release a hydrogen atom to form C2H4. The existence of radicals has been detected via electron spin-resonance (ESR) experiments conducted by Kolboe et al.,53 who have recorded the ESR signals of a ZSM-5 catalyst during the induction period. Former experimental works have suggested that surface methoxyl groups may account for the first C2 formation. IR and NMR investigations have indicated that the broken C−H bond of the surface methoxyl group is responsible for the formation of the first hydrocarbons during the MTO process.43,44 Methoxyl groups may play the role of precursors of surface-

Figure 4. Product distributions as a function of simulation time in model A (methanol/acidic zeolite); the numbers of species are the average of three simulations. The lower figure zooms in the species with a small number. High-temperature regions are not shown.

molecules are released, and the molecular weight increases. Subsequently, unsaturated carbohydrates are formed. An analysis of the simulated trajectories reveals the mechanisms in detail. The observed major reaction pathways are illustrated in Figure 5. First, the methanol molecule is

Figure 5. Mechanisms for the activation of methanol molecules observed during molecular dynamics simulations. Species in red are only observed in high-temperature regions. 7033



Article

Figure 6. Formation of free CH3OH2 groups by collision. (a) A methanol absorbed on the acid site is attacked by a free methanol molecule, (b) collision between the two molecules, and (c) release of a CH3OH2 group and formation of a deprotonated site.

Figure 7. Two pathways of formation of the first C2 species in model A (only methanol as reactant). I(a) A methanol absorbed on the acid site is attacked by a free CH3 group; I(b) the C−O bond of methanol stretches and breaks, resulting in the formation of a water molecule; and I(c) A C−C bond is formed and a water molecule is absorbed at the deprotonated site. II(a) A surface methoxyl group is attacked by a free CH3; II(b) a C−C bond is formed, and one C−H bond from the methoxyl group stretches and breaks; and II(c) an acid site is restored, and a C2H5 species is produced.

stabilized carbene or ylide.1,45 Using 13C MAS NMR-UV/vis spectroscopy analysis, Wang et al.7 have reported that methoxyl groups work as methylation agents of cyclohexane via a “carbene-like species”. Very recently, Yamazaki et al.46 have presented evidence of a “carbene-like intermediate” for a carbon chain growth reaction between ethene and surface methoxyl groups.46 These experimental data are consistent with the C2 formation mechanisms observed in our simulations. 3.3. Model B. The simulation performed on Model B was designed to study how the carbon chain grows by reacting with the methoxyl groups, which are known to be methylation

agents in the reactions.42,43 Figure 8 shows the averaged population of major species as a function of time obtained in the three parallel 1 ns simulations, along with the temperature curve. Compared with methanol decomposition, the rate of ethene reaction is relatively low. At the end of the simulation, only five ethene molecules are consumed. Associated with the disappearance of ethene, hydrocarbon species with C3 and C4 appear. The first C3 is immediately generated as the first ethene disappears, indicating that C3 directly originates from C2. The first C4 species appears at ca. 300 ps, and its appearance is associated with the disappearance of C3, indicating that C4 is 7034



Article

Figure 10. Snapshot of an absorbed ethene molecule on a deprotonated site to form a C2H4···Z structure.

Figure 8. Product distributions as a function of the simulation time in Model B (ethene/methoxylate zeolite). The numbers of species are the average of three simulations. High-temperature regions are not shown.

produced from the C3 species. The chart also includes C1 and C2 species, which are tracked as different species from the methoxyl C1 and feeding ethene C2. They are intermediates that exist in the early stage of the simulations. The fact that C1, C2, C3, and C4 species appear in the order as given indicates chain growth as the simulation proceeds. The observed major reaction channels are illustrated in Figure 9. The most common one is the same as that reported in the literature.38 In this route, an ethene molecule collides with a methoxyl group to form a CH2CH2CH3 carbenium. The carbenium undergoes intramolecular proton transfer to form a more stable carbenium, CH3CHCH3, which loses a hydrogen atom to produce propene. Another less popular pathway is observed, as illustrated in the figure. The ethene molecule is captured by a deprotonated Brønsted site to form an adsorbed complex C2H4···Z. Figure 10 is a snapshot of the structure. As shown in the figure, one hydrogen atom of ethene is within a short distance of ca. 1.85 Å to the Brønsted oxygen, and the distance between the hydrogen and carbon atom stretches to 1.13 Å. The adsorbed structure easily loses a hydrogen atom to the framework and forms a C2H3 species. The C2H3 species reacts with other species in the system to form Cn+2 species. Considering that the C2H4···Z route has not been reported in the literature, DFT calculations were carried out to validate this mechanism. The energy changes in the hydrogen transfer from C2H4 to the deprotonated Brønsted site were calculated using a 3T cluster model and the B3LYP/6-31 g(d,p) method. Two electronic models, the open shell for radicals and closed shell for ionic species, were evaluated. The closed shell results indicate a very high energy barrier height (120.4 kcal/mol), whereas the open shell results are much more modest. Figure 11 illustrates the energy diagram obtained using the open shell model. The adsorption energy is a little less than 5 kcal/mol; the energy barrier height is only 5.8 kcal/mol; and the overall process is slightly endothermic. A previous report has stated4,47

Figure 11. QM validation of the hydrogen transfer process from ethene to a deprotonated site, a 3T cluster at the B3LYP/6-31 g(d,p) level. The data are calculated using the open shell (radical) model.

that the deprotonated Brønsted oxygen lacks the strong basic character to break the C−H bond of methanol, but the unsaturated hydrocarbons such as ethene have not been taken into consideration. Our data show that the C2H4 molecule can be easily converted into the C2H3 radical by the deprotonated Brønsted site. 3.4. Comparison of the Initial Reaction Rates. With the same simulation protocol, the decomposition of methanol (Model A) is significantly faster than that of the ethane (Model B). The kinetic characteristics of both reactions were derived based on simulations performed at different temperatures from 1500 to 2000 K at a 100 K interval. The rate constants at different temperatures were obtained by fitting the logarithm of number densities to straight lines. As shown in Figure 12, the logarithms of the reaction constants are plotted versus the reciprocal of temperature (1/T). The effective activation

Figure 9. Mechanisms of ethene chain growth. Species in red are only observed in high-temperature regions. 7035



Article

Figure 13. Representation of collisions between reactants: (a) collisions between a methanol molecule and a H site in both good and bad directions and (b) collision between an ethene molecule and a CH3 site in both good and bad directions.

vdW surfaces, the differences were quantified by defining the effective collision probability (ECP) between atoms A and B ECP = Figure 12. Arrhenius extrapolation for (a) methanol activation and (b) ethene activation. The temperature range is from 1500 to 2000 K.

where ESA is the effective accessible surface area and TASA is the total accessible surface area. The accessible surface area for a molecule is calculated as the van der Waals surface area of the molecule accessible by a probing atom (H, O, or C). With the vdW radius of each element and molecular equilibrium structures, the ECP values were calculated to be 9.1% for methanol and Brønsted hydrogen but only 1.8% for ethene and the methoxyl group. 3.5. Model C. Model C that contains both methanol and ethene molecules mimics a condition wherein a HP is involved in the reactions. The populations of important species identified in the simulations are shown in Figure 14. Both methanol and ethene molecules start to react from the beginning. Although the decay rate is higher for methanol than for ethene, the rate for methanol consumption is slower than that found in Model A, but the rate of consumption of ethene is faster than that in Model B. Generally, the coexistence of methanol and ethene enables rapid chain growth. As shown in the chart, 5 of 6 ethene molecules are consumed at 250 ps. Along with the consumption of ethene, the populations of C1, C2, and C3 appear in the order of the number of carbon. The patterns of these populations show that they are intermediates that increase to a maximum and then decrease. This phenomenon is not clear from the data obtained for Model B (Figure 8), presumably due to the slower rate of chain growth C2 in Model B than that in Model C. The water amount steadily increases but reaches a plateau near the end, coincident with the disappearance of methanol. This finding is consistent

energies Ea and frequency factors A are obtained from the fitted lines. The data are summarized in Table 2. The apparent Table 2. Comparison of Apparent Activation Energies (kcal/ mol) Obtained by the Extrapolation of MD Simulations and Static DFT Calculations methanol activation method MD DFT

Ea 40.7 40.6

ethene activation

A 3.6 × 10

14

Ea

A

33.6 26.4

2.2 × 1013

⎛ ESA ⎞ ⎛ ESA ⎞ ⎜ ⎟ ⎟ ×⎜ ⎝ TASA ⎠ A ⎝ TASA ⎠B

activation energies calculated from the simulation are 40.7 kcal/ mol for methanol and 33.6 kcal/mol for ethene. The results are in reasonable agreement with previous DFT calculation results.38,48−51 However, the difference between the frequency factors is significant. As shown in the table, the frequency factors are 3.6 × 1014 for methanol decomposition, 1 order of magnitude higher than the value of 2.2 × 1013 for ethene decomposition. Intuitively, methanol activation only takes place when the lonepair electron of oxygen atoms contacts the Brønsted hydrogen, as illustrated in Figure 13a. Carbon activation only occurs when the electrons of the π-orbital in ethene collide with the center of the CH3 group of the methoxyl group (Figure 13b). Assuming that the effective collision is proportional to exposed 7036



Article

Figure 14. Product distributions as a function of the simulation time of Model C (methanol/ethene/acidic zeolite). The numbers of species are the average of 15 simulations. The lower figure zooms in on the species with a small number. High-temperature regions are not shown. Figure 15. Mechanism collection from a complex model. The numbers in brackets indicate the number of times the reactions occur in 15 parallel complex system simulations. Species in red only exist in high-temperature regions.

with the mechanism of water generation as discussed above. The population of H2 molecules steadily increases, accompanied with the generation of large unsaturated hydrocarbon molecules, which will be further scrutinized later in this paper. 3.6. Reaction Network. On the basis of a careful analysis of the trajectories obtained from 15 independent simulations performed on Model C, a reaction network is drawn, as shown in Figure 15. The network consists of two cycles that are interconnected in a sense that the intermediates generated in one cycle may be used in another. The lower part is the methanol cycle in which the methanol is consumed to generate various intermediates. The upper cycle is the chain-growth cycle in which higher-order hydrocarbon chains are generated. Between the two cycles is H−Z, which represents the acidic zeolite. To its left is one of the feeding molecules, ethene; to its right is another, methanol. The major species observed in the simulations are denoted as nodes, and the reaction paths identified in the simulations are represented by lines that connect the nodes. The fraction numbers along the paths indicate the relative probabilities of different reactions from a species (node). The denominator indicates the number of species marked at the previous node that reacted in the 15 simulations. The numerator represents how many of them take the indicated path. All of the reactions start from methanol being adsorbed on a Brønsted acid site. For a total of 225 methanol molecules (15 simulations each having 15 methanol molecules), 210 CH3OH2−Z complexes were found. The majority (183/210) of them went through the common route by releasing a H2O molecule and forming a methoxyl group CH3−Z on the surface. The CH3−Z groups are fairly stable and serve as intermediates in various reactions in the network. The most common one

(165/186) interacts with another hydrocarbon to form a new C−C bond. A much less popular (27/210) route, as discussed above (Figure 6), is via the vdW complex CH3OH2−Z, forming an active species CH3OH2, which releases a water molecule and becomes a free CH3 radical. The CH3 radicals are quickly associated with nearby species. Among the CH3 radicals identified in the simulations, 9/27 form CH4 by combining a hydrogen atom from the acidic site of zeolite, leaving the deprotonated acidic center (Z), which also plays a very important role in the reaction network; 3/27 combine with the deprotonated acidic center (Z) to form a methoxyl group; 9/27 form new C−C bonds with another hydrocarbon species; and 6/27 interact with a methoxyl group to form the C2 species, C2H5. C2H5 reacts with a deprotonated acidic site (Z) and returns a hydrogen atom to the framework to form a C2H4 molecule, and the catalyst is recovered. Interestingly, the last step generates an ethene molecule. However, only six new C2H4 molecules are generated in 15 simulations. This result means that the chance of directly observing ethene formation from C1 species in 1 ns is only 6/15. A total of 96 ethene molecules including those generated from the methanol cycle are in the simulations, and 67 of them are reacted. The reactions for ethene can be classified as three types depending on the coreactant: with the methoxyl (C1) group, with higher-order hydrocarbons (C2+), and with water. Surprisingly, only a small portion (5/67) of the ethene molecules follows the well-known route by reacting with the methoxyl group to form a C3 species. The majority (62/67) of 7037



Article

Figure 16. Hydrogen formation process identified during molecular dynamics simulations, starting from (a) an ethene molecule attacking a Brønsted acid site, (b) hydrogen molecule formation after the breakage of C−H and O−H bonds, and (c) bonding of a C2H3 group to the framework to form a C2H3−Z complex.

the ethene molecules interacts with the deprotonated Brønsted site to form an intermediate, C2H4···Z, and then becomes a C2H3 radical, as discussed in the Model B simulations. Among the C2H3 radicals, 9/62 interact with C1 species (methoxyl group) to form C3 species, and 36/62 interact with large carbon species (>C2) to form C4+ species. The Model A and B simulations do not reveal that 16/62 of them interact with water, forming C2H3OH2. Ten of the 16 C2H3OH2 molecules release a hydrogen atom to the framework and become the aldehyde C2H4O. 3.7. Hydrogen and Coke Formation. A significant amount of hydrogen molecules are produced during the late stage of our simulations. The observed mechanisms for generating hydrogen are similar. The hydrogen atom on the Brønsted site is attacked by a hydrocarbon CmHn to form a transition state in which the H2 molecule forms. The carbon atom of the hydrocarbon then bonds to the zeolite framework to form an alkyl oxide CmHn‑1−Z. The alkyl oxide may undergo further reactions with other species. This mechanism can be seen by a snapshot taken from the simulation, as illustrated in Figure 16. An ethene molecule collides with one Brønsted acid site; the C−H bond of ethene and the O−H bond of zeolite are broken; and two hydrogen atoms form a new H2 molecule, leaving ethyl oxide (C2H3−Z). The dehydrogenation rate depends on the size of the hydrocarbon. Figure 14 clearly shows that hydrogen molecules are produced more rapidly during the late stage of simulations when higher-order hydrocarbons are already formed. This finding stems from the stabilization of the transition state of the dehydrogenation process by the longer carbon chain. The identified process is similar to the dehydrogenation mechanism of propane obtained via the transition pathway sampling method.52 Accompanying the dehydrogenation is the formation of hydrocarbons with high degrees of unsaturation. Figure 17 is animage of a large hydrocarbon molecule formed near the end of the simulations. The molecule has a structure close to that of benzene, a typical compound of coke molecules.

Figure 17. Types of high C/H ratio species identified during the late stage of molecular dynamics simulations (>500 ps).

been analyzed, and new reaction mechanisms are identified. (2) By considering collision probabilities, the rate-determining step for MTO is not the activation of methanol, as suggested by static QM calculations, but the C−C chain formation between ethene and methyloxide. Due to the higher collision probability, methanol decomposition is faster than the carbon chain growth. However, the activation energy of methanol decomposition is higher than that of the carbon chain growth. (3) For the activation of methanol, a less possible but very important CH3 radical mechanism is identified in addition to the commonly accepted methoxyl mechanism. The commonly referred chain-growth mechanism in which ethene interacts with methyloxide to form a new C−C bond has been observed. However, this is a small contribution to the total production. The most popular route for ethene activation is due to the deprotonated Brønsted sites that are produced by the activation of methanol molecules. (4) Therefore, the HP is effective when methanol molecules are involved. Without methanol, the chain growth is a slow process. With methanol, which undergoes decomposition and produces deprotonated acidic sites, the hydrogen pool starts to work. Consequently, chain growth is significantly accelerated.

4. CONCLUSIONS (1) The ReaxFF force field has been extended for the simulation of catalytic reactions in acidic zeolite. To produce sufficient data for the analysis of the reaction mechanism, a simulated annealing technique has been used in the simulations. Using three simulation models, the initial reaction network of MTO in acidic zeolite has

■

ASSOCIATED CONTENT

S Supporting Information *

ReaxFF force field parameters used in this work. This material is available free of charge via the Internet at http://pubs.acs.org. 7038



■

Article

(31) Zhang, Q.; Qi, Y.; Hector, L. G. Jr.; Ç ağın, T.; Goddard, W. A. III Phys. Rev. B 2005, 72, 045406(12p). (32) Olson, D. H.; Kokotailo, G. T.; Lawton, S. L.; Meier, W. M. J. Phys. Chem. 1981, 85, 2238−2243. (33) Alvarado-Swaisgood, A. E.; Barr, M. K.; Hay, P. J.; Redondo, A. J. Phys. Chem. 1991, 95, 10031−10036. (34) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (35) Rappe, A. K.; Goddard, W. A. J. Phys. Chem. 1991, 95, 3358− 3363. (36) van Duin, A. C. T. B., J. M. A.; van de Graaf, B J. Chem. Soc. 1994, 90, 2881−2895. (37) 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 C.02; Gaussian, Inc.: Wallingford, CT, 2004. (38) Maihom, T.; Boekfa, B.; Sirijaraensre, J.; Nanok, T.; Probst, M.; Limtrakul, J. J. Phys. Chem. C 2009, 113, 6654−6662. (39) Brändle, M.; Sauer., J. J. Am. Chem. Soc. 1998, 120, 1556−1570. (40) Clark, L. A. S., M.; Sauer, J. J. Am. Chem. Soc. 2003, 125, 2136− 2141. (41) Derouane, E. G. C., C. D. Microporous Mesoporous Mater. 2000, 425, 35−36. (42) Cui, Z.-M.; Liu, Q.; Bain, S.-W.; Ma, Z.; Song, W.-G. J. Phys. Chem. C 2008, 112, 2685−2688. (43) Wang, W.; Buchholz, A.; Seiler, M.; Hunger, M. J. Am. Chem. Soc. 2003, 125, 15260−15267. (44) Jiang, Y.; Hunger, M.; Wang, W. J. Am. Chem. Soc. 2006, 128, 11679−11692. (45) Hutchings, G. J.; Hunter, R. Catal. Today 1990, 6, 279−306. (46) Yamazaki, H.; Shima, H.; Imai, H.; Yokoi, T.; Tatsumi, T.; Kondo, J. N. Angew. Chem. 2011, 123, 1893−1896. (47) Lesthaeghe, D.; Van Speybroeck, V.; Marin, G. B.; Waroquier, M. Angew. Chem. 2006, 118, 1746−1751. (48) Andzelm, J.; Govind, N.; Fitzgerald, G.; Maiti, A. Int. J. Quantum Chem. 2003, 91, 467−473. (49) Vos, A. M.; Nulens, K. H. L.; De Proft, F.; Schoonheydt, R. A.; Geerlings, P. J. Phys. Chem. B 2002, 106, 2026−2034. (50) Svelle, S.; Rønning, P. O.; Kolboe, S. J. Catal. 2004, 224, 115− 123. (51) Svelle, S.; Rønning, P. O.; Olsbye, U.; Kolboe, S. J. Catal. 2005, 234, 385−400. (52) Bučko, T. J. Chem. Phys. 2009, 131, 214508. (53) A. Holmen, K.-J. J. a. S. K. Nat. Gas Convers. 1991, 413−420.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was partially funded by the National Science Foundation of China (No. 21073119 and No. 21173146) and the National Basic Research Program of China (No. 2007CB209701).

■

REFERENCES

(1) Stocker. Microporous Mesoporous Mater. 1999, 29, 3−48. (2) Blaszkowski, S. R.; van Santen, R. A. J. Am. Chem. Soc. 1997, 119, 5020−5027. (3) Tajima, N.; Tsuneda, T.; Toyama, F.; Hirao, K. J. Am. Chem. Soc. 1998, 120, 8222−8229. (4) Lesthaeghe, D.; Van Speybroeck, V.; Marin, G. B.; Waroquier, M. Ind. Eng. Chem. Res. 2007, 46, 8832−8838. (5) Mokrani, T.; Scurrell, M. Catal. Rev.: Sci. Eng. 2009, 51, 1−145. (6) Song, W.; Marcus, D. M.; Fu, H.; Ehresmann, J. O.; Haw, J. F. J. Am. Chem. Soc. 2002, 124, 3844−3845. (7) Wang, W.; Hunger, M. Acc. Chem. Res. 2008, 41, 895−904. (8) Dahl, I. M.; Kolboe, S. J. Catal. 1994, 149, 458−464. (9) Dessau, R. M. J. Catal. 1986, 99, 111−116. (10) Haw, J. F.; Marcus, D. M. Top. Catal. 2005, 34, 41−48. (11) Haw, J. F.; Song, W.; Marcus, D. M.; Nicholas, J. B. Acc. Chem. Res. 2003, 36, 317−326. (12) Svelle, S.; Joensen, F.; Nerlov, J.; Olsbye, U.; Lillerud, K.-P.; Kolboe, S.; Bjørgen, M. J. Am. Chem. Soc. 2006, 128, 14770−14771. (13) Bjørgen, M.; Svelle, S.; Joensen, F.; Nerlov, J.; Kolboe, S.; Bonino, F.; Palumbo, L.; Bordiga, S.; Olsbye, U. J. Catal. 2007, 249, 195−207. (14) Bjørgen, M.; Joensen, F.; Lillerud, K.-P.; Olsbye, U.; Svelle, S. Catal. Today 2009, 142, 90−97. (15) Zhu, Q.; Kondo, J. N.; Setoyama, T.; Yamaguchi, M.; Domen, K.; Tatsumi, T. Chem. Commun. 2008, 5164−5166. (16) Vandichel, M.; Lesthaeghe, D.; Mynsbrugge, J. V. d.; Waroquier, M.; Van Speybroeck, V. J. Catal. 2010, 271, 67−68. (17) Schmidt, K.-H.; Jurado, B. Phys. Rev. Lett. 2010, 104, 212501(4p). (18) Ivanova, E.; Pang, J.; Jowitt, T. A.; Yan, G.; Warwicker, J.; Sutcliffe, M. J.; Lu, H. Proteins: Struct., Funct., Bioinf. 2011, 80, 602− 615. (19) Liu, L.; Bai, C.; Sun, H.; Goddard, W. A. J. Phys. Chem. A 2011, 115, 4941−4950. (20) Tilocca, A. J. Mater. Chem. 2010, 20, 6848−6858. (21) Chen, H.; Yan, T.; Voth, G. A. J. Phys. Chem. A 2009, 113, 4507−4517. (22) Ohta, Y.; Okamoto, Y.; Page, A. J.; Irle, S.; Morokuma, K. ACS Nano 2009, 3, 3413−3420. (23) Christie, J. K.; Tilocca, A. Chem. Mater. 2010, 22, 3725−3734. (24) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A. J. Phys. Chem. A 2003, 107, 3803−3811. (25) Nomura, K.-i.; Kalia, R. K.; Nakano, A.; Vashishta, P. Comput. Phys. Commun. 2008, 178, 73−87. (26) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. Angew. Chem. 2009, 121, 7766−7770. (27) Mayernick, A. D.; Batzill, M.; van Duin, A. C. T.; Janik, M. J. Surf. Sci. 2010, 604, 1438−1444. (28) Neyts, E. C.; Shibuta, Y.; van Duin, A. C. T.; Bogaerts, A. ACS Nano 2010, 4, 6665−6672. (29) Voter, A. F. Phys. Rev. Lett. 1997, 78, 3908−3911. (30) Sugita., Y; Okamoto., Y. a. Chem. Phys. Lett. 1999, 314, 141− 151. 7039


Molecular Dynamics Simulations of Methanol to Olefin Reactions in

Recommend Documents