Transition State Modeling for Catalysis - American Chemical Society

Chapter 28

Acidic Catalysis by Zeolites: Ab Initio Modeling of Transition Structures Joachim Sauer, Marek Sierka, and Frank Haase

Downloaded by STANFORD UNIV GREEN LIBR on August 7, 2012 | http://pubs.acs.org Publication Date: April 8, 1999 | doi: 10.1021/bk-1999-0721.ch028

Humboldt-Universität, Institut für Chemie, Arbeitsgruppe Quantenchemie, Jägerstraße 10-11, D-10117 Berlin, Germany

Two different methods for treating transition states in large catalytic systems are presented: (i) a combined quantum mechanics - empirical valence bond (QM– EVB) method for localizing transition structures in embedded clusters and (ii) the constraint Car-Parrinello molecular dynamics for a fully periodic ab initio treatment. The power of the methods is demonstrated for two different elemen tary reactions in zeolite catalysts, proton jump between different proton localiza tion sites and C-C formation between two adsorbed methanol molecules yieldig ethanol and water.

Zeolites are microporous materials which are extensively used in the chemical industry as heterogeneous catalysts. They contain Bronsted acidic OH groups, SiO(H)Al(OSi) , which represent the active sites and are capable of activating chemical bonds. While the different zeolites have this type of active site in common, they differ by the structure of the (alumo)silicate framework and, hence, by the shape and size of the micropores. Theoretical studies of the catalytic activity involve different stages of increasing complexity. (1) The Bronsted acidity is characterized by calculating energies ofdeprotonation. This has been done for different framework structures (1, 2). This is a direct, but hypothetical process and does not involve a transition structure. (2) The catalytic activity may be related to the proton mobility within the cavity of the unloaded zeolite. Introducing an Al atom into the three-dimensional cornersharing net of Si0 tetrahedra creates a negative charge which is compensated by the acidic proton. The proton can be attached to any of the four oxygen atoms of the tetrahedron. Depending on the particular framework considered, the different proton positions have different ener gies. It is known that, e.g., in H-faujasite protons occupy preferentially the 03 and 01 sites. In this study, we are interested in the transition structure and the energy barrier for proton jumps between framework oxygen atoms belonging to the same A10 tetrahedron. 3

4

4

358

© 1999 American Chemical Society

In Transition State Modeling for Catalysis; Truhlar, D., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1999.


359 (3) The proton jumpfromthe Brensted site to an adsorbed molecule which is activated for further reaction by protonation. We have shown for methanol by different methods that for a loading of one molecule per Brensted site the protonated methanol in the zeolite cavity is a low lying transition structure or a local minimum about 10-20 kJ/mol above the neutral methanol hydrogen-bonded to the zeolite surface (3, 4). (4) The ultimate interest is in elementary steps of the catalytic conversion of the feedstock molecules activated by protonation. We are particularly interested in methanol which is the feedstock of the two industrially important processes, methanol to olefin (MTO) and methanol to gasoline (MTG). For the methanol conversion different mechanisms have been proposed. They are based on transition structures and intermediates localized by density functional calculations which represent the catalysts by a small cluster model (5). A particularly interesting step in this mechanism is the formation of thefirstC-C bond which we study here. Localizing transition structures is routine now for gas phase reactions, but there are substantial difficulties for reactions involving solid catalysts. A realistic approach which aims at reactivity differences between different zeolites with differentframeworkstruc tures must include periodic boundary conditions and take the structure of the micropores into account. This means that the number of degrees offreedomis very large which makes the localization of transition structures very demanding. We present two strategies to overcome such difficulties. Thefirstis the use of a combined quantum mechanics - inter atomic potential approach (QM-Pot; hybrid or embedded cluster method) for getting in formation on the potential energy surface. To allow bond breaking and making in the interatomic potential part we adopt Warshel's empirical valence bond (EVB) method (6). Transition structures are localized by the a modification of the trust region method (7-9). For so many degrees offreedom(up to 435 in this study) this can be completed only be cause an initial transition structure and an initial Hessian are calculated with the EVB po tential functions only. The final transition structure search is made using combined QME VB gradients which are also used to update the Hessian. We demonstrate the success of the method by localizing transition structures for a well defined elementary step, proton jumps between differentframeworkoxygen atoms in two different zeolite structures (stage 2 above). Significantly different barriers are predicted for different crystallographic sites in a given zeolite and for zeolites with differentframeworkstructure. The second strategy aims at exploring larger regions of a potential energy surface with arichstructure of local minima and saddle points. It is the method of constraint mo lecular dynamics (10) which we use in combination with a full ab initio calculation of the zeolite applying periodic boundary conditions (Car-Parrinello MD) (11). Such simula tions are computationally very demanding and presently limited to catalysts with modest unit cell sizes. We study the C-C bond formation between two methanol molecules in cha bazite (stage 4 above). Combined QM-EVB Studies of Proton Jumps in Zeolites Methods. The embedding scheme used (1, 12) partitions the entire system (S) into two parts: the inner part (I) containing site in question, and the outer part (O). Saturating the dangling bonds of the inner part with link atoms (hydrogen atoms) yields the cluster (C). The energy of the total system is obtained approximatelyfromthe QM energy of the clus-


360 ter and the difference between the energies of the total periodic system and the cluster, calculated by the interatomic potential functions. Hence, when calculating the forces on an atom a of the I region all three energies make contributions: F„ , = F (C) + F (S)-F (C)


e

0M

Pot

(1)

Pot

For an atom of the O region, the force is obtainedfromthe interatomic potential function alone. A link atom is not moved according to the force acting on it. It is instead kept fixed on the bond which it terminates (12). Our subtraction scheme (equation 1) requires an analytical interatomic potential func tion both for the outer and the inner region. For treating reactions including transition structure (TS) searches this becomes a serious problem. Interatomic potential functions generally are not capable of describing bond breaking and bond making processes. The empirical valence bond (EVB) approach of Warshel (6) is a solution which we adopt. The single minimum potential energy surfaces (PES) of the reactant and product states are coupled within an empirical 2x2 valence bond problem to yield the adiabatic PES which describes both minima connected by the transition structure. The crucial part of the EVB model is the exchange element V which we calculate following Chang and Miller (13): l2

2

T

T

[Vi (r)] = A exp (B Ar + A r C Ar) 2

(2)

where Ar = r - r , r is the current structure, r is a reference structure in terms of internal coordinates, and A, B and C are a constant scalar, vector and matrix, respectively, which are derived using quantum mechanically calculated energies, gradients and Hessian for the TS. Differently from the original proposal we define V in terms of internal coordi nates of the atoms with the largest displacement along the reaction path only, and not of all the atoms. This allows us to derive values of A, B, C and r using quantum mechanical calculations on small cluster models for the TS. These values are assumed to be transfer able and used to describe the reaction in the periodic lattice. Our modification of the meth od relies on the observation that the reaction coordinate involves only negligible motions of atoms sufficiently distantfromthe reaction site. We use a modification of the trust region method (7-9) with BFGS (7) and GMSP (14) update methods of the Hessian for localizing minima and transition structures, re spectively. First, the optimum cell parameters for the periodic zeolite lattice are deter mined using potential functions alone (constant pressure optimization). Then a constant volume optimization of all degrees offreedomis performed within the combined QMEVB scheme. The initial Hessian matrix calculated at the EVB level is updated in every iteration. We stress that localizing the TS is achieved without explicit calculation of the Hessian for the QM part, which in many cases would be computationally extremely ex pensive. Convergence in the gradient norm of0.0001 E a,/ is usually achieved within 12 cycles for the minimum and within 20 cycles for the TS search (H-faujasite, 435 degrees offreedom).For structures with a smaller QM region (3 or 4 tetrahedra) the exact com bined QM-EVB Hessian matrices are evaluated at the end of the search to proof that a sta tionary point of correct order is found. For the QM part we apply the density functional (DFT) method and adopt the B3LYP (15) functional. The basis set is triple-zeta on oxygen and double-zeta on all other 0

0

1 2

0

1

h


361


atoms. Polarization functions are added to all atoms. The interatomic potential functions have the functional form of a shell-model ion-pair potential which has the advantage of including the mutual polarization of the inner and outer regions. They are parametrized on B3LYP calculations using the same basis set as for the QM cluster (16). The present implementation^/ 7) of the QM-EVB scheme uses the GULP program(18) for evaluating the interatomic potential function contributions and the TURBOMOLE (19, 20) program for the quantum mechanical calculations. Problem and Models. Experimetally, the jump rate and the barrier are not easily acces sible. The zeolite has to be completelyfreeof any adsorbed molecule and it is not easy to completely dehydrate and deammoniate the samples. NMR experiments have been per formed to estimate the activation energy for proton jumps in zeolites, but the results are scanty and inconsistent (21-24). For example for faujasite type zeolites the barriers of 2040 kJ/mol (Si/Al = 1.2 and 2.6) (21) and 61 kJ/mol (Si/Al = 3) (24) have been reported. Early quantum mechanical calculations by Sauer et al. (25) used finite cluster models and arrived at an estimate of 52± 10 kJ/mol for the proton jump barrier. This applies to an "av erage" zeolitic Bronsted site, but cannot explain differences between different zeolites.

Figure 1. Fragments of the faujasite (a) and chabasite (b)frameworksand topologies around an Al site (c and d, respectively) which are identical for both frameworks, except the numbering.



362 The combined QM-EVB method presented above properly accounts for the structure of the periodic lattice and the long range crystal interactions at modest computational ex pense. Specifically, we tackle the problem of differences between the proton jump barri ers between different crystallographic positions within the H-faujasite lattice and between two different zeolite frameworks, H-faujasite and H-chabasite. In the unit cells of faujasite, Figure la, and chabasite, Figure lb, all the tetrahedral atoms are crystallographically equivalent but there are four inequivalent oxygen posi tions, labelled 01 - 04. Bothframeworksconsists of four- and six-membered silicate rings and include the double six-membered ring (hexagonal prism) as secondary building unit. Theframeworktopology around the Al atom is similar for bothframeworks(Figure lc). We perform QM-EVB calculations using different sizes of the QM region (cluster model). Thefirstconsists of three T 0 tetrahedra (3T model, Figure 2a). For H-faujasite we define a large QM region consist of 3 four-membered silicate rings (4T model, Figure 2b). For comparison, we also examine the finite cluster model used by Sauer et. al. (25) (Figure 2c). 4

3

Results and Discussion. Table I shows bond angles and distances which characterize the proton jump transition structure. For 3T cluster models embedded into the two framework structures the exact combined QM-EVB Hessian marices have been calculated and we confirmed that they have one negative eigenvalue (imaginaryfrequenciesof about 1300 cm" ). For the larger models the Hessian initialized by the EVB method, but updated using gradientsfromthe combined QM-EVB method, always retained one negative eigenvalue 1


363 during the whole optimization process. In all the structures the proton, the Al atom and the two bridging oxygen atoms form a plane within 5 degree. The proton is located be tween the two oxygen atoms with O-H distances 120.4 - 127.7 pm. The formation of TS involves large distortion of the O-Al-0 tetrahedral angle by about 30 degrees. This dis tortion is larger for thefreecluster model due to the lack offrameworkconstraints. The structures obtained for the 01-04 jump in H-faujasite using the 3T and 4T cluster models are virtually identical.


3

Table I. Bond Distances (pm) and Angles (Degree) of Proton Jump Transition Structures. zeolite

jump sites QM cluster

free cluster H-chabasite H-faujasite

r(0'-H)

r(Al-H)

O-Al-O'

122.5

122.5

194.9

74.8

123.3 127.7 127.0

126.3 125.3 125.1

194.3 194.4 194.2

76.9 78.0 77.6

187.2

77.7

Ol-H-02 Ol-H-04 Ol-H-04

4T

3

Ol-H-02

4T

3

126.7

120.4

4T

3

123.6

121.6

185.7

77.0

4T

3

123.4

121.2

184.0

77.7

123.7

124.5

194.3

76.3

Ol-H-03 03-H-04 03-H-02

3T 3T

r(0-H)

4 X

3

Table II shows the reaction energies and activation barriers calculated for the proton jumps. The total QM-EVB energies can be decomposed into the quantum mechanical contribution, AE . , and the long range contribution AE . both evaluated for the structure obtainedfromQM-EVB optimizations. For almost the same cluster size the quantum mechanical contribution of the embedded cluster is larger than the barrier for thefreecluster which is due to the constraint imposed by the framework on the struc ture adjustment for the transition structure. Both clusters consists of three tetrahedra, but differ in the termination, OH for the embedded cluster and H for the free space cluster. Moreover, the long range interaction also makes a positive and significant contribution to the energy barrier. Since both effects neglected in thefreespace clusters work in the same direction, the results for embedded cluster calculations are substantially larger than the free cluster result. This is different from deprotonation and ammonia adsorption energies for which the two effects act in opposite direction (1, 2). Increasing the quantum mechanical part from the 3T to the 4T model reduces the long range contribution to 1 /3 while the total barrier changesfrom107 to 101 kJ/mol only. This is the expected behaviour as in the limit of a very large model the long range correc tion should vanish. For H-chabasite we could afford fully periodic single point calcula tions using the Dsolid code (26). The QM-EVB and the periodic QM results agree within 1 kJ/mol. (#)

QM//QM

EVB

(#)

LR//QM

EVB

3


364 #

Table II: Reaction Energies, AE, and Activation Barriers, A E , as well as Quan tum Mechanical and Long Range Contributions to the Barrier, QM//QM-EVB and LR//QM-EVB, for Proton Jumps in Zeolites, kJ/mol. positions QM cluster

zeolite


free cluster

QM

LR

AE

AE

0.0

42.0

42.0

-

51.2

24.5

#

H-chabasite

Ol-H-02

3T

21.7

75.7

H-faujasite

Ol-H-04

3T

22.2

107.0

66.0

41.0

Ol-H-04

4T

3

13.8

100.9

87.0

13.9

Ol-H-02

4T

3

19.6

98.2

77.8

20.4

Ol-H-03

4T

3

5.8

105.6

89.5

16.1

03-H-04

4T

3

8.0

109.6

94.0

15.6

03-H-02

4T

3

13.8

66.9

59.9

7.0

The calculations reveal that proton jump barriers show significant differences be tween different zeolites, i.e. between H-chabasite and H-faujasite, and there are also large differences between the different sites and paths in a given zeolite. In H-faujasite the low est barrier found is 67 and the highest one 110 kJ/mol (4T model, Table 2). It has been suggested (24) that the proton jump barrier is proportional to the acidity differences be tween the two oxygen sites involved in the jump. Our results do not support this assump tion. For example, for the 01-04 and 02-03 jumps in H-faujasite the energy barriers differs by more than 30 kJ/mol while the relative energies of the two minima are the same (13.8 kJ/mol). These relative energies correspond directly to differences in the deproto nation energies, the usual measure of acidity. 3

Constraint Car-Parrinello Molecular Dynamics for C-C Bond Formation between Two Methanol Molecules in H-Chabasite Simulations. Constant temperature MD simulations at 500 K are performed. First, unconstraint simulations are made for the initial state of two adsorbed methanol molecules and the product state, adsorbed ethanol and water molecules. Because the likely reaction path is known to involve a rather high activation barrier the method of the 'Blue Moon ensemble is adopted (27). It requires the choice of a holonomic constraint which is intro duced into the Car-Parrinello Lagrangian. The C-C distance is the natural choice of a reaction coordinate for C-C bond forma tion. By performing several MD simulations for different values of the constraint the re action coordinate is scannedfromthe average C-C distance in the reactant to the average C-C distance in the product state. Hence, this method does not localize a single transition structure, but samples the phase space in its vicinity. Free energy differences are then computed by integrating the statistically averaged constrained forces along the chosen re action coordinate. The electronic problem is solved by density functional theory employ ing the Perdew-Burke-Ernzerhof functional. We use Vanderbilt ultrasoft pseudopotentials and expand the Kohn-Sham orbitals in a plane wave basis set including all plane waves with kinetic energies up to 25 Ry. 1



365 Results. We have carried out constraint MD simulations for six values of the C-C dis tance between 150.8 pm, the equilibrium distance of the ethanol molecule, and 359.8 pm, the C-C distance in the adsorbed methanol dimer. The following observations are made: (1) The zeolitic proton is immediately transferred to the methanol dimer and virtually no fluctuations are seen which return the proton. (2) At a C-C distance near the transition state (210 - 230 pm) the C-0 bond of the protonated methanol molecule becomes significantly elongated. (3) The transition state structures comprise a carbonium ion and a weakly interacting wa ter molecule. Configurations are seen in which the proton of the carbonium ion is shared between the two carbon atoms similar to a non-classical bridged ethyl cation.

Figure 3. Constraint MD simulation of two methanol molecules in H-chabasite. Snapshot showing a configuration close to the transition state. (4) The water molecule is split offfromthe transition state complex. (5) Eventually the proton of the carbonium ion is transfered back to the zeolite framework restoring the acid site. Later on the ethanol molecule forms a hydrogen-bonded surface complex with its OH group pointing to the zeolitic proton. Figure 4 shows the statistically averaged constrained forces along the reaction coor dinate. Integration from the initial state to the transition state yields afreeenergy of acti vation of about 200 kJ/mol. To estimate the error connected with the relatively small number of integration points we made an additional MD simulation at r(CC)=4.1 bohr close to the transition state. Thefreeenergy of activation increased by 20 kJ/mol. Hence, ourfinalestimate of the free energy barrier is 220±20 kJ/mol. We also computed an en ergy profile of the reaction by constrained optimization (zero Kelvin) of representative configurations taken from the 500K trajectories. According to this profile the activation energy is 180 kJ/mol. Blaszkowski et al. (5) calculated activation energies for the same reaction, but re placed the periodic zeolite structure by a small cluster model, HOHAl(OH) OH. It com prises just the acidic proton and the A10 -tetrahedron terminated by hydrogen atoms. The activation energy calculated was 310 kJ/mol, significantly larger than our value. This is not surprising because the transition structure is an ionic complex which in our periodic scheme can be stabilized by the zeoliticframework.However, the transition structure found in the cluster calculations is similar to that found in our CPMD simulations. 2

4


366

transition state

0.09

J3 o

0.07

i

0.05


1

0.03 0.01

2MeOH

EtOH+H?0

O

o

-0.01

3.5

4.5

5.5

6.5

Figure 4. Statistically averaged constrained forces as a function of the C-C distance (cubic spline of 6 points).

Acknowledgement This work has been supported by the "Max-Planck-Gesellschaft" and the "Fonds der Chemischen Industrie". Computer time on a CRAY-T3E has been provided by the "Hochstleistungsrechenzentrum Stuttgart (HLRS)". Literature Cited 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.

Eichler, U.; Brändle, M.; Sauer, J. J. Phys. Chem. B 1997, 101, 10035. Brändle, M.; Sauer, J. J. Am. Chem. Soc. 1998, 120, 1556. Haase, F.; Sauer, J. J. Am. Chem. Soc. 1995, 117, 3780. Haase, F.; Sauer, J.; Hutter, J. Chem. Phys. Lett. 1997, 266, 397. Blaszkowski, S. R.; van Santen, R. A. J. Am. Chem. Soc. 1997, 119, 5020. Aqvist, J.; Warshel, A. Chem. Rev. 1993, 93, 2523. Fletcher, R. Practical Methods of Optimization; John Wiley&Sons: New York, NY, 1987. Helgaker, T. Chem. Phys. Lett. 1991, 182, 503. Culot, P.; Dive, G.; Nguyen, V. H.; Ghuysen, J. M. Theoret. Chim. Acta 1992, 82, 189. Ryckaert, J. P.; Cicotti, G. J. Chem. Phys. 1983, 78, 7368. Galli, G.; Parrinello, M. In Computer Simulations in Material Science; Meyer, M.; Pontikis, V., Eds.; Kluwer: Dordrecht, 1991; pp 283-304.



367 12. Eichler, U.; Kölmel, C.; Sauer, J. J. Comput. Chem. 1997, 18, 463. 13. Chang, Y.-T.; Miller, W. H. J. Phys. Chem. 1990, 94, 5884. 14. Bofill, J. M. Chem. Phys. Lett. 1996, 260, 359. 15. Becke, A. D. J. Chem. Phys. 1993, 98, 5648. 16. Sierka, M.; Sauer, J. Faraday Discuss. 1997, 106, 41. 17. Sierka, M. QMPOT; Humboldt-University: Berlin, 1998. 18. Gale, J. D. J. Chem. Soc., Faraday Trans. 1997, 93, 629. 19. Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Chem. Phys. Lett. 1989, 162, 165. 20. Treutler, O.; Ahlrichs, R. J. Chem. Phys. 1995, 102, 346. 21. Freude, D.; Oehme, W.; Schmiedel, H.; Staudte, B. J. Catal. 1974, 32, 137. 22. Freude, D. Chem. Phys. Lett. 1995, 235, 69. 23. Baba, T.; Inoue, Y.; Shoji, H.; Uematsu, T.; Ono, Y. Microporous Mater. 1995, 3, 647. 24. Sarv, P.; Tuherm, T.; Lippmaa, E.; Keskinen, K.; Root, A. J. Phys. Chem. 1995, 99, 13763. 25. Sauer, J.; Kölmel, C. M.; Hill, J.-R.; Ahlrichs, R. Chem. Phys. Lett. 1989, 164, 193. 26. DSOLID; Molecular Simulations Inc.: San Diego, 1994. 27. Carter, E. A.; Ciccotti, G.; Hynes, J. T.; Kapral, R. Chem. Phys. Lett. 1989, 156, 472. 28. Munson, E. J.; Kheir, A. A.; Lazo, N. D.; Haw, J. F. J. Phys. Chem. 1992, 96, 7740.


Transition State Modeling for Catalysis - American Chemical Society

Recommend Documents