Modeling Catalytic Steps on Extra-Framework Metal Centers in

Article pubs.acs.org/JPCC

Modeling Catalytic Steps on Extra-Framework Metal Centers in Zeolites. A Case Study on Ethylene Dimerization Shrabani Dinda,†,§ Agalya Govindasamy,†,§ Alexander Genest,† and Notker Rösch*,†,‡ †

Institute of High Performance Computing, Agency for Science, Technology and Research, 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore ‡ Department Chemie and Catalysis Research Center, Technische Universität München, 85747 Garching, Germany S Supporting Information *

ABSTRACT: In a case study of organometallic catalytic reactions, this work benchmarks density functional theory calculations on zeolite-supported transition metal complexes. Elementary steps of ethylene dimerization and hydrogenation reactions involving the complex [Rh(C2H4)2(H2)]+, supported on faujasite, were examined by comparing explicit QM (quantum mechanics) cluster models as well as QM/MM (molecular mechanics) embedded models to planewave periodic models as reference. Two QM cluster models, 1T and 5T where T refers to tetrahedral units of zeolite, as well as four QM/ MM cluster models were explored. For the MM region, the UFF force field was found preferable to the semiempirical method PM6. The embedded cluster models reproduce barriers of C−C and C−H bond formation with deviations from the reference of at most 10 kJ mol−1. With variations of similar size, the effect of embedding on the energetics of the reactions under study is moderate, likely because of the small nonpolar reactants. For elucidating such catalytic reactions at transition metal species in zeolites, cluster models appear equally well-suited as periodic models but computationally advantageous. rected functionals like B97-D26 and ωB97X-D.31 The ability of such more advanced DFT methods to reproduce experimental adsorption enthalpies in zeolites were probed in a few cases (e.g., for the empirical DFT-D approach32−34 and the nonlocal correlation functional vdW-DF2).32 As small cluster models of zeolites are unable to describe long-range interactions and specifically the diversity of zeolite frameworks, large clusters or periodic (first-principles) calculations may be considered as more appropriate. However, these latter methods entail a notably higher computational effort. Seeking a balance between accuracy and computational cost, embedded cluster methods have been explored20,35−42,43−48 where a partition of the zeolite system, including the reaction site, is treated by a high-level (often DFT) method. Most such strategies employ hybrid quantum mechanical/molecular mechanical (QM/MM) approaches.35,37,41,42,44,46,48 The ONIOM method,20,38,39,45,47 in particular the QM/shell model,36,40,43 are well-known for adequately representing zeolite systems. Besides the computational method (or method combination), also the size of cluster models is an important aspect when addressing catalysis in zeolites. As key characteristic of such models, it is common to state the number of T atoms in the QM partition, where T refers to tetrahedral units of zeolite.

1. INTRODUCTION Zeolites are industrially important catalysts for a wide range of applications due to their microporous architecture and their unique acidic characteristics.1−6 Hydrocarbon transformations are very common among zeolite-supported catalytic reactions. For producing higher alkenes via C−C coupling, significant progress has been made using acidic zeolites as an alternative to conventional acid catalysts. 1,7−10 Recent experimental work11−16 addressing C−C coupling on oxide- or zeolitesupported transition metal complexes motivated the present benchmark examination of computational methods. As a computational case study, our work targets aspects of experimental investigations on the dimerization of ethylene to n-butene over site-isolated Rh(I) complexes, supported on dealuminated zeolite HY.13−15 The Kohn−Sham approach to density functional theory (DFT) provides efficient computational methods for describing such systems and has been widely used in computational catalysis.17−22 However, most standard DFT methods, specifically due to the exchange-correlation approximations commonly used, are not able to deal with dispersion interactions (van der Waals, vdW) 23 which limits the application of such methods to many zeolite-catalyzed reactions.24,25 In recent years, numerous efforts have been described to construct DFT methods that are able to handle vdW interactions (e.g., empirical variants such as the DFT-D scheme),26,27 new functionals like the M06 family28,29 and vdW-DF,30 or self-consistently parametrized dispersion-cor© 2014 American Chemical Society

Received: August 11, 2014 Revised: September 30, 2014 Published: October 22, 2014 25077

dx.doi.org/10.1021/jp508141q | J. Phys. Chem. C 2014, 118, 25077−25088

The Journal of Physical Chemistry C

Article

A few computational studies49−51 have been reported that specifically address C−C coupling reactions in acidic zeolites. Using 4T cluster models of the Brønsted acid site, a B3LYP study demonstrated a concerted mechanism for ethylene dimerization with an activation barrier of 127 kJ mol−1.51 A three-partition ONIOM approach, (MP2:HF:UFF), of a 84T cluster of faujasite revealed a stepwise process, with a comparable barrier of 121 kJ mol−1.50 The effect of the environment on the dimerization of ethylene was discussed using a [B3LYP:MNDO] ONIOM strategy for similarly large cluster models of the zeolites SSZ-13 (74T), ZSM-5 (72T), and ZSM-22 (66T).49 Previous computational studies on zeolite-supported Rh complexes,19,52 related to the experimental work of Gates et al.,19,52 used 1T cluster models to assign vibrational spectra52 and to explore the mechanism of acetylene cyclo-trimerization.19 To the best of our knowledge, there are no computational studies on the mechanism of ethylene dimerization on zeolite-supported transition metal complexes. The present work aims at establishing a reliable, yet economic strategy for describing chemical reactions on zeolite-supported transition metal complexes. For the present case study, we selected C−C coupling and hydrogenation reactions, catalyzed by a [Rh(C2H4)2(H2)]+ complex adsorbed on faujasite, and designed a hierarchical set of models of the active site. Using results of DFT-D calculations on periodic models as a reference, we explored how different levels of cluster embedding affect structures and energetics.

2.2. Method Aspects of Periodic Models Calculations. The DFT calculations on periodic models were carried out with a plane-wave approach, employing the Vienna Ab initio Simulation Package (VASP).57 The projector augmented wave (PAW) method58 was used to describe the electron− ion−core interaction. As an exchange-correlation functional we selected the Perdew−Burke−Ernzerhof (PBE) formulation59 of the generalized gradient approximation (GGA). To describe dispersion interactions, we invoked the semiempirical D2 scheme,26,27 as implemented in VASP. The conjugate-gradient algorithm was used to relax the locations of the atoms within the unit cell until the forces acting on each atom were below 0.02 eV/Å. In a few cases, we selected the quasi-Newton (QN) algorithm to speed up the convergence, once the gradients with the CG method were about 0.05 eV/Å. The kinetic energy cutoff was set to 800 eV for the full optimization of the unit cell and to 400 eV for the relaxation of the ionic positions in the fixed unit cell.32 For the SCF iterations, we selected thresholds of the changes in the total energy of 10−6 eV for minima and 10−7 eV for transition states. The sampling of the Brillouin zone was restricted to the Γ point. A Gaussian smearing (parameter 0.01 eV) was used for the occupation numbers; total energies were extrapolated to zero level broadening.60 Transition state (TS) structures were located with the dimer method,61 starting from structures obtained by running a few steps of the nudge elastic band method.62 All minima and TS structures were verified by a restricted normal-mode analysis using the finite difference method63 implemented in VASP. Only the fragment of the metal complex was included in this analysis while all framework atoms were kept fixed. Including the Al center and the connected −O−Si units in addition to the Rh complex changed the vibrational frequencies of interest at most by 5 cm−1. All intermediates studied were in a singlet state, as confirmed by single-point plane-wave calculations on the triplet states of the reference models. In all cases, the corresponding triplet states are notably less stable, by more than 200 kJ mol−1 (Table S1 of the Supporting Information). 2.3. Method Aspects of Cluster Models Including ONIOM Details. For the cluster calculations, we chose four models: 1T, 5T, 83T, and 240T (Figure 1). We used the first two as pure QM models and the latter two for cluster embedding. The 1T model is composed of the fragment [Al(OH)4]− and the moiety [Rh(C2H4)2(H2)]+ which is anchored on the oxygen centers of two of the hydroxyl groups.19,52 The zeolite fragment of the 5T model consists of an [Al(OSiH3)4]− unit in which the Si atoms at the “edges” are saturated with hydrogen atoms. The 83T cluster, 278 atoms in total, is taken as a segment of a supercage unit with additional sodalite cages and hexagonal prisms. The 240T cluster, 768 atoms in total, comprises a full supercage surrounded by 10 sodalite cages. O centers on the surface of these embedding models are at least 400 pm from the Rh center. Dangling bonds of Si centers at the outer edges of the cluster models 83T and 240T (including the Rh adsorption complex) were saturated with hydrogen atoms. Each Si−H bond was oriented to the nearest (missing) oxygen center of the ideal zeolite structure. Initially the positions of the terminal hydrogen centers were optimized by relaxing only the Si−H distances, while the rest of the zeolite framework was kept fixed at the experimental structure. Subsequently, the corresponding H and Si positions were kept fixed to prevent artifacts due to unrealistic deformations of the cluster models.

2. COMPUTATIONAL MODELS AND METHODS 2.1. Zeolite Models. Faujasite has a highly symmetric cubic structure that belongs to the space group Fd3m ̅ . Its unit cell contains 192 T atoms. There is a single crystallographically distinct T site surrounded by four symmetry-distinct oxygen centers (Figure S1 of the Supporting Information).53 The faujasite structure is built from supercages, sodalite cages, and hexagonal prisms which are made up of 4-, 6-, and 12-member rings (MR). To reduce the computational cost of the periodic calculations, caused by the rather large unit cell, we invoked a strategy54,55 (Figure S1b of the Supporting Information) that employs a primitive cell of 48 T atoms with a volume that is only one-fourth of the original unit cell. The translational vectors of this smaller primitive cell, of rhombohedral shape, accordingly are reduced by a factor 1/√2, and the angles between the vectors are reduced to 60°. The crystal structure of all-silica faujasite was taken from the IZA Web site.56 The primitive cell Si48O96 of all-silica faujasite was optimized (Method Aspects of Periodic Models Calculations), yielding the lattice parameters a = b = c = 17.32 Å. To introduce acidic centers, we replaced one of the symmetry-equivalent Si atoms by an Al atom, yielding a framework with a ratio Si/Al = 47. The net negative charge of the lattice due to this substitution was compensated by a proton on a bridging oxygen center of the type Al−O−Si. The optimized reduced unit cell of the Al substituted faujasite model HAlSi47O96 is triclinic, and the lattice parameters are a = 17.35 Å, b = 17.34 Å, c = 17.32 Å; α = 60.2°, β = 60.0°, γ = 60.0°. For modeling the zeolite-supported cationic complex [Rh(C2H4)2]+, formally Rh(I), we kept this shape of the unit cell and replaced the acidic proton in the lattice by the [Rh(C2H4)2]+ fragment (see Method Aspects of Periodic Models Calculations and Faujasite-Supported Complex [Rh(C2H4)2]+). 25078



Article

Vibrational frequencies were calculated at the same level of theory as the corresponding geometry optimization. The structures reported are either local minima with no imaginary frequency or transition states with exactly one imaginary frequency. In some cases, we determined the internal reaction coordinate (IRC) to verify that the obtained TS structure lies on the path that connects the intended initial and final state structures. To evaluate calculated structures, we used rootmean-square deviations (RMSD) of the distances X-Rh to all other atoms X of the moiety [Al(O)2Rh(C2H4)2(H2)] (including the two O centers of the zeolite framework bound to Rh) with respect to the optimized structure of the corresponding periodic model.

3. RESULTS AND DISCUSSION For the present case study, we selected five typical elementary reactions that may occur during ethylene dimerization catalyzed by a [Rh(C2H4)2(H2)]+ complex adsorbed on faujasite (Scheme 1): C−C coupling as direct dimerization of ethylene, ethylene insertion into a Rh-ethyl bond, and ethylene hydrogenation to ethyl. The results for C−H bond formation Scheme 1. Schematic Representation of Five Reactions under Studya

Figure 1. Sketches of various cluster models used to study C−C and C−H formation reactions on the complex [Rh(C2H4)2]+, anchored at an Al center in the supercage of faujasite. QM partitions of the various models are given in ball-and-stick fashion; bonds in the embedding (low-level) partition are rendered as simple sticks. For instance, model 11/83T comprises a QM partition of 11 T atoms out of a zeolite cluster model of 83T atoms in total. Si−H bonds at the cluster edges are shown explicitly. The color coding of the atoms is given at the bottom of the figure.

All cluster calculations were carried out with Gaussian09.64 Stationary points of the models 1T and 5T were located using the PBE exchange correlation functional62,59 and 6-31G(d,p) basis sets65 for H, C, O, Al, and Si, as well as the Stuttgart− Dresden effective core potential MWB2866 for Rh. We adopted an ONIOM scheme45 with two partitions of the models 83T and 240T. For each of them, we selected high-level (QM) partitions of two sizes, 5T and 11T (Figure 1), to check the effect of this model aspect. The rest of the zeolite clusters were assigned as low-level partition. The resulting ONIOM models are labeled as 5/83T, 11/83T, 5/240T, and 11/240T. The high-level partition of the ONIOM calculation was treated with the same level of theory used for the models 1T and 5T. To describe the environment of the high-level partitions, we used the semiempirical PM6 method67 or the MM universal force field (UFF)68 approach. The combined theoretical method ONIOM(PBE:UFF) was applied to all embedding cluster models. The 5/83T models optimized in this way were also treated by single-point QM calculations for the full 83T cluster. We chose the DFT functionals PBE, PBE-D2,27 B97D,26 and ωB97X-D.31 In addition, ONIOM(PBE:PM6) calculations were carried out on the 5/83T model to check the effect of the method chosen for the low-level partition.

a

Reaction 1→2 represents a concerted C−C coupling and H2 activation; reaction 3→4 is a concerted C−C coupling and H−H bond formation, reactions 5→6 and 5′→6′ are C−C coupling reactions involving ethylene insertion into the Rh−ethyl bond, forming a butyl ligand, and 3→5 is a C−H bond formation reaction. The ethyl group of 6 shows an agostic interaction of ethyl with the Rh center, whereas complex 6′ does not exhibit such an interaction. 25079



Article

Can one suggest an economic cluster model that yields sufficient accuracy for studying this type of catalysis problem? 3.1. Faujasite-Supported Complex [Rh(C2H4)2]+. Experimental studies as well as previous theoretical investigations on small cluster models suggest that zeolite-anchored Rh(I) complexes bind to two O centers of the acidic Al site.13,15,19,52,69 A QM/MM study of [Rh(CO)2]+ binding on faujasite, addressing the confining effect of the zeolite cage, showed that this metal complex prefers to bind to two framework oxygen atoms of the 12-MR over other bidentate sites on 4- and 6-MR.70 In the present work on [Rh(C2H4)2]+, we chose the same bidentate binding site of the 12-MR of faujasite. In the periodic reference model, the [Rh(C2H4)2]+ moiety binds to the framework in an asymmetric fashion, Rh−O = 218 and 225 pm (Table S3 of the Supporting Information). A similar asymmetric binding mode was determined with all cluster models, except 1T and 5T, emphasizing the need for an adequate description of the environment of the active site (Table S3 of the Supporting Information). The Rh−O distances in both bare cluster models are notably shorter than in the reference: 211 pm in 1T and 215 pm in 5T. The Rh−O distances in the embedded cluster models falls into the ranges 220−223 and 226−229 pm, thus are slightly longer than in the reference structure. The remaining structural parameters, Rh− C and C−C, of [Rh(C2H4)2]+ deviate in all cluster models at most by 1−2 pm from the corresponding values of the periodic model (Table S3 of the Supporting Information). We also compared characteristic vibrational frequencies of faujasite-supported [Rh(C2H4)2]+ complexes to experimental data52 (Tables S4 and S5 of the Supporting Information). After suitable scaling with respect to the experimental spectrum of free ethylene (Table S4 of the Supporting Information), one finds that the C−H stretching frequencies of the complex, calculated with the periodic model, fit the experimental spectrum of the complex within 8 cm−1, but the largest deviation is 24 cm−1 (Table S5 of the Supporting Information). More revealing for the metal−ligand bonding in [Rh(C2H4)2]+ is the C−C stretching mode, measured at 1623 cm−1 in free ethylene. The present as well as earlier DFT calculations52 show that, in the adsorption complex, the locally totally symmetric CH2 scissoring mode of ethylene, measured at 1342 cm−1 in the free molecule, couples to the C−C mode. In the zeolite embedded complex [Rh(C2H4)2]+, only the two modes (one for each ligand) with the higher frequency were

are given in Section S3 and Table S2 of the Supporting Information. Figure 2 shows the TS structures of the four C−C coupling reactions, as obtained with periodic DFT-D calculations. We

Figure 2. Optimized structures of four TSs for the C−C coupling reactions under study. Selected distances in picometers. Color coding: Rh, blue; O, red; C, gray; and H, white.

will base the evaluation of methods and models on calculated reaction energies Er (Table 1) and activation barriers Ea (Table 2) as well as on the corresponding deviations ΔEr and ΔEa from the results of the corresponding periodic DFT-D calculations. Table S3 of the Supporting Information provides information on structures of intermediates and TSs. With such reactions in mind, we examined geometric and energetic aspects of cluster models, by addressing the following questions using periodic models as reference. (1) Does a QM/ MM approach provide an adequate description of the active site and its reactivity? (2) What level of embedding is advisable? (3)

Table 1. Reaction Energies Er (kJ mol−1) of Four C−C Coupling Reactions (Scheme 1) on the Complex [Rh(C2H4)2]+, Anchored at an Al Center in the Supercage of Faujasite, from Various Modeling Strategies,a Calculated at the PBE Levelb models 1T 5T 5/83T, UFF 5/83T, PM6 11/83T, UFF 5/240T, UFF 11/240T, UFF Periodic, PBE-D

1→2

3→4

5→6

MAD

5′→6′

Er

ΔEr

Er

ΔEr

Er

ΔEr

Er

ΔEr

14c −1 −7 −27 −3 6 4 14

0c −15 −21 −41 −17 −8 −10 −

7 −2 11 −6 4 16 10 10

−3 −12 1 −16 −6 6 0 −

11 13 19 44 24 18 19 31

−20 −18 −12 13 −7 −13 −12 −

5 5 9 46 6 7 8 14

−9 −9 −5 32 −8 −7 −6 −

8.0 13.5 9.8 25.5 9.5 8.5 7.0

The embedding environment is treated at the UFF or PM6 level. bAlso shown are the differences ΔEr of the activation energies relative to the corresponding values from the periodic models and their mean absolute deviation (MAD). cReaction energy estimated by constraining the initial state to the structure optimized for complex 2.

a

25080



Article

Table 2. Activation Barrier Ea (kJ mol−1) of Four C−C Coupling Reactions (Scheme 1) on the Complex [Rh(C2H4)2]+, Anchored at an Al Center in the Supercage of Faujasite, from Various Modeling Strategies,a Calculated at the PBE Levelb models

1→2

1T 5T 5/83T, UFF 5/83T, PM6 11/83T, UFF 5/240T, UFF 11/240T, UFF Periodic, PBE-D

3→4

5→6

MAD

5′→6′

Ea

ΔEa

Ea

ΔEa

Ea

ΔEa

Ea

ΔEa

108 96 102 100 103 106 98 93

15 3 9 7 10 13 5 −

132 130 129 116 129 132 132 142

−10 −12 −13 −26 −13 −10 −10 −

136 137 146 165 143 141 137 145

−9 −8 1 20 −2 −4 −8 −

94 94 93 132 92 92 91 90

4 4 3 42 2 2 1 −

9.5 6.8 6.5 23.8 6.8 7.3 6.0

The embedding environment is treated at the UFF or PM6 level. bAlso shown are the differences ΔEa of the activation energies relative to the corresponding values from the periodic models and their mean absolute deviation (MAD).

a

Table 3. Root Mean Square Deviation (RMSD, pm) of the Distances X-Rh to All Other Atoms X of the Moietya [Al(O)2Rh(C2H4)2(H2)] with Respect the Optimized Structure of the Corresponding Periodic Model model 1T 5T 5/83T PM6 11/83T 5/240T 11/240T a

TS

TS

TS

TS

1

1−2

2

3

3−4

4

5

5−6

6

5′

5′-6′

6′

6 5 8 10 8 6 6

6 6 4 16 8 12 6

16 22 18 42 13 20 11

3 2 4 6 4 4 5

19 19 19 18 19 19 19

3 3 5 8 4 6 6

5 6 5 9 4 6 6

4 4 4 32 4 6 6

5 5 4 6 5 5 6

5 3 4 5 4 5 5

4 3 8 9 4 9 6

5 4 7 23 7 8 9

Including the two O centers of the zeolite framework bound to Rh.

3.2. Low-Level Embedding of the QM Partition. To examine the effect of the low-level method in ONIOM calculations, used for treating the environment of the highlevel partition, we compare results of (PBE:UFF) calculations on the embedded cluster model 5/83T to results of (PBE:PM6). For all four C−C coupling reactions (Scheme 1), the models treated with the (PBE:UFF) combination reproduce the energetics of the periodic reference model better than models treated with the (PBE:PM6) approach (Tables 1 and 2). In all cases, but one, the deviations ΔEr of the reaction energies obtained with the PM6 method are largest by absolute value (Table 1); the same finding holds for the deviations ΔEa of the activation energies (Table 2). As a global measure for the impact of the description of the environment on the overall structure, we compared RMSD results from both types of ONIOM approaches for the fragment [Al(O)2Rh(C4H10)] in the QM partition of the embedded cluster model 5/83T to the corresponding PBE-D results of the periodic reference (Table 3). The RMSD value, 15 pm, of the ONIOM(PBE:PM6) method, averaged over all reactions, is somewhat larger than the corresponding result ONIOM(PBE:UFF), 10 pm. In summary, embedding in a partition modeled by the PM6 method does not seem to be a good choice for describing reactions in a zeolite framework, especially when the active site comprises a transition metal center. Therefore, in the following, we will discuss only ONIOM results obtained with the combination (PBE:UFF). This conclusion agrees with results of earlier studies where UFF as low-level method had been found to work surprisingly well for studying C−C coupling and adsorption in acidic zeolites.50,71,72 3.3. Selected Elementary Steps of C−C Coupling. Now we turn to the four C−C coupling reactions (Scheme 1). In

identified by infrared (IR) spectroscopy as a broader peak, at 1440 cm−1.52 Thus, in the complex [Rh(C2H4)2]+, one notes a significant shift to lower frequencies due to back-donation from the Rh center into the ethylene π* orbitals. Our calculations on the periodic reference model give the corresponding scaled frequencies at 1475 and 1454 cm−1; the average of these two frequencies is 25 cm−1 higher than the experimental value. The modes with the lower frequencies, calculated at 1195 and 1181 cm−1, are not accessible via IR measurements.52 Cluster models reproduce these frequencies of the periodic reference system rather well (Table S5 of the Supporting Information). In contrast, the average of the two higher frequencies, calculated with the embedded cluster models, is about 20−35 cm−1 larger than in the reference system and their splitting a bit smaller than in the reference. All cluster models, even those without embedding, reproduce the two lower frequencies within just a few wavenumbers. The cluster models 1T and 5T reproduce the calculated C− H stretching frequencies of the periodic model within 8 cm−1 with an average absolute deviation of ∼3 cm−1 (scaled values, Table S5 of the Supporting Information). The same models reproduce the four C−C stretching frequencies slightly worse, with an average absolute deviation of ∼7 cm−1 from the result of periodic models. Cluster embedding notably shifts the C−H stretching frequencies to higher values, by 10−20 cm−1, increasing the deviations from the periodic model by 9−36 cm−1 (Table S5 of the Supporting Information). The higher two of the C−C stretching frequencies are even more strongly overestimated, with embedding by up to 42 cm−1 [model 11/ 83T, ONIOM(PBE:UFF)], while the lower two frequencies of all embedded systems remain within 3 cm−1 at the values of the periodic model, resulting in average deviations of 11−20 cm−1. 25081



Article

each case, we first discuss the structural parameters and the chemistry based on the periodic DFT-D reference models. Subsequently, we explore the efficiency of the cluster strategies by comparing in detail geometric and energy characteristics of various cluster models. The results of the C−H bond formation are provided in Section S3 of the Supporting Information. To probe temperature and entropy effects on the energetics, we calculated the free energy of reaction 5′→6′ (to be discussed later) and the corresponding free energy of activation as an example, comparing the results for the 5/83T embedding model to those of the periodic reference, at standard conditions, 298.15 K. These results agree within 0.1 kJ mol−1 with those at the experimental temperature, 303 K. Reaction and activation free energies, calculated both for the periodic reference and the cluster model, differed only by 1−4 kJ mol−1 from the corresponding energy results. Hence, we decided to restrict the discussion to energy values. Reaction 1→2. With an H2 molecule coadsorbed at the active site, the dimerization of ethylene can proceed either before or after H2 activation via oxidative addition of H2 to the Rh center. The reaction 1→2 involves oxidative C−C coupling with simultaneous activation of an η2-coordinated H2 molecule (Scheme 1). These are two oxidative additions where a Rh(I) complex with 16e− is directly oxidized to a 16e− Rh(V) species, concomitantly reducing chemisorbed H2 to dihydride and the two ethylene moieties to metallocyclopentane (MCp), respectively. In species 1, H2 binds to the metal center in a side-on fashion: Rh−H = 190 and 184 pm; H−H = 82 pm (Table S3 of the Supporting Information). For the resulting MCp-dihydride complex 2, we calculated Rh−H = 154 pm and H−H = 155 pm. In the imaginary vibrational mode (361i cm−1) of the corresponding transition state 1−2, the centers C2 and C3 move toward each other; simultaneously, the centers H1 and H2 move moderately away from each other. The structures 1−2 shows evidence of C−C bond formation with the distance C2−C3 reduced from 286 to 206 pm and the H− H distance enlarged from 82 to 90 pm (Figure 2). We calculated a reaction energy Er = 14 kJ mol−1 and an activation energy Ea = 93 kJ mol−1 (Tables 1 and 2). The QM clusters 1T and 5T reproduce key structural parameters of complex 1 and TS 1−2 in adequate fashion (Table S3 of the Supporting Information). The embedded cluster models show similar trends for the structural parameters of the active site, except for an increased elongation in Rh−H1 bonds of complex 1. Noticeable differences in the 5/240T model in the TS 1−2 are the elongated Rh−C2/C3 distances, by 15−18 pm (Table S3 of the Supporting Information). The carbon centers C2 and C3 are those involved in the C−C coupling. We quantified the structural change at the active site by RMSD values of atomic positions (relative to Rh) in the fragment [Al(O)2Rh(C2H4)2(H2)] (Method Aspects of Cluster Models Including ONIOM Details, Table 3). In general, low RMSD values correlate with low values of ΔEr and ΔEa. The RMSD values for this reaction scatter in the range of 4−22 pm, where structure 1 and TS 1−2 yield smaller values, while complex 2 gives higher values in all cases. The scattering of the RMSD values is due to the fact that the MCp ring of 2 in the periodic model is more puckered than in the cluster models. In the 1T model, the distance H1−H2 in 2 is constrained to 155 pm, to ensure close structural similarity to the periodic model. Without this constraint, the optimization yields a butyl complex that is energetically more favorable by 39 kJ mol−1, with an agostic interaction of the type Rh−H2−C4 (Figure S2

of the Supporting Information). In contrast to the 1T model, the 5T model retains the MCp configuration, although the RMSD is high, 22 pm. This deviation is mostly due to shifts of the hydrogen atoms attached to C3. Cluster embedding does not improve the RMSD value in any of the models. With the smaller embedding model 5/83T, the distances Rh−C3 and Rh−H3 are calculated noticeably longer than with the periodic reference, Rh−C3/H3 = 277/294 pm versus 248/238 pm, respectively. These structural findings may reflect that the agostic interaction is calculated too weak in the cluster model (Figure S2 of the Supporting Information). The high RMSD value, 18 pm in the product complex 2 for the 5/83T model, also reflects these changes. To check whether this discrepancy is due to differences in treating dispersion interaction, we reoptimized the 5/83T structure at the PBE-D level, starting from the configuration of the periodic model. However, the resulting optimized structure is similar to the one previously obtained without any D correction. Comparing structures with the 5/83T model as reference, the RMSD values of the other embedded-cluster models are much smaller, ranging from 2 to 8 pm only (Table S6 of the Supporting Information). Thus, the structure of the periodic model differs slightly while all cluster models yield similar structures, irrespective of the level of embedding. Recall that periodic slab model calculations predict the reaction 1→2 as slightly endothermic. All cluster models, except 1T, show negative reaction energies ΔEr, from −8 to −21 kJ mol−1 (Table 1). The higher ΔEr values, −8 and −10 kJ mol−1 are calculated for the models 5/240T and 11/240T, while the lowest deviations ΔEr are determined for the models 5/83T and 11/83T at −21 kJ mol−1 and −17 kJ mol−1, respectively. Thus, the larger embedding cluster with 240 T atoms leads to reaction energies somewhat closer to the periodic reference. Although the 1T and 5T models are certainly too small for describing many aspects of the zeolite systems, the calculated activation barriers are quite close to the reference values of the periodic model, with ΔEa = 15 and 3 kJ mol−1, respectively (Table 2). Most embedded cluster models yield very satisfactory activation energies, ΔEa = 5−10 kJ mol−1. A slightly larger deviation, ΔEa = 13 kJ mol−1, is obtained with the model 5/240T. The productlike TS structure, according to the Rh−C2/C3 distances that are elongated by 15−18 pm with respect to the periodic reference (Table S3 of the Supporting Information), occurs only for this model. Therefore, we interpret these findings as an indication for a somewhat different TS isomer with a slightly larger discrepancy in ΔEa. Reaction 3→4. In the reaction 3→4, the C−C coupling is probed after H2 activation by the metal center (Scheme 1). The transformation 1→3 of a η2-H2 complex to a dihydride complex via oxidative addition is quite common and straightforward in organometallic catalysis.73,74 In the reaction 3→4, the dihydride 3, with Rh−H = 154 and 158 pm and H−H = 203 pm (Table S3 of the Supporting Information), is oxidized to dihydrogen, and the ethylene moieties are reduced to MCp. Thus, the metal center retains its +3 oxidation state. Just as TS 1−2, also TS 3− 4 is an “early” transition state regarding the hydrogen configuration (Figure 2). In the periodic model, going from the initial state 3 to TS 3−4 the Rh−H distances remain essentially unchanged and the H−H separation is hardly reduced, from 203 to 196 pm, while the distance C2−C3 shortens notably, from 308 to 177 pm (Table S3 of the Supporting Information). The resulting complex 4 comprises a 25082



Article

MCp moiety and an η2-coordinated dihydrogen, with C2−C3 = 153 pm and H−H = 93 pm. A normal-mode analysis of TS 3− 4, with one imaginary frequency at 290i cm−1, confirms this transformation. The reaction 3→4 is endothermic by 10 kJ mol−1 (Table 1). The activation energy of 3→4, Ea = 142 kJ mol−1, is notably larger than the barrier of reaction 1→2, Ea = 93 kJ mol−1 (Table 2). Now we compare structures of the cluster models to those of the periodic models. The RMSD of the structures of the active site is uniformly quite high, 19 pm for all cluster models (Table 3). The largest individual deviations, 31−45 pm, are calculated for the distances Rh−HC2 and Rh−HC3 (Figure S3 of the Supporting Information). As before, the periodic model differs notably from all cluster models. With the use of the 5/83T structure as reference, the RMSD values of TS 3 → 4 for other cluster models, embedded or not, drop to 1−4 pm (Table S6 of the Supporting Information). The dihedral angle C1−C2−C3− C4 is 46° in the periodic model and about −47° (!) in all cluster models (Table S3 of the Supporting Information). One notes a short contact in the periodic model, between a hydrogen at C2 of one ethylene and a framework oxygen, HC2− O = 234 pm, hinting at a van der Waals interaction between HC2 and the zeolite framework. In the embedded cluster models, this particular oxygen center is directly at the boundary between the QM and MM partitions, potentially entailing an underrepresentation of the interaction with the zeolite wall. This may explain why the optimization of embedded models yields a different arrangement of the four carbon centers. To corroborate this assumed underrepresentation of the interaction with the zeolite wall, we look at the favorable arrangement of these four carbon atoms in the models without zeolite framework. Indeed, with the bare models 1T and 5T, we calculated essentially the same arrangement between the two ethylene moieties as in the embedded models, as exemplified by the dihedral angle C1−C2−C3−C4 (Table S3 of the Supporting Information). Therefore, we conclude that the embedded cluster models do not provide the same interaction pattern of the ethylene moieties with the zeolite walls as the periodic model. We also probed the stability of the arrangement of the two ethylene units. The conformation of the two ethylene moieties determined with the periodic model but used with the 1T model, differs only by 4 kJ mol−1 from the genuine 1T conformation. Indeed, both arrangements may easily be transformed into each other, with an approximate barrier of 10 kJ mol−1 (or 6 kJ mol−1 in reverse direction). Also, the distance H1−H2 is different in the two types of models: 196 pm in the periodic reference and ∼187 pm in the 1T cluster model. A detailed analysis of the product structure 4 did not reveal notable RMSD values, which are at most 6 pm in all cluster models (Table 3). The reaction energies for the C−C coupling 3→4 are relatively well-reproduced by the various cluster models, with a largest (absolute) deviation ΔEr occurring in the 5T model of −12 kJ mol−1. The other deviations ΔEr are in the range from −6 to 6 kJ mol−1. Similarly, one does not calculate any significant embedding effect on the activation barriers TS 3−4 (Table 2); they all are by ∼12 kJ mol−1 lower than in the periodic model. This systematic deviation, ΔEa being negative for all cluster models, correlates with the special structure of TS 3−4, calculated by the periodic model (see above). Reactions 5→6 and 5′→6′. Next, we discuss the C−C coupling between a Rh-coordinated ethylene moiety and an ethylene ligand. In these cases, the hydrogenation 3→5 of

ethylene (Table S2 of the Supporting Information) precedes the C−C coupling step (Scheme 1). In both reactions 5→6 and 5′→6′, an ethylene moiety is inserted into the Rh-ethyl bond to form the butyl ligand while maintaining the +3 oxidation state of the metal center. Each of the two pairs, 5 and 5′ as well as 6 and 6′, comprises rotational isomers. Complex 5 is 39 kJ mol−1 more stable than 5′ due to an agostic interaction between the CH3 group of the ethyl ligand and the metal center. Complex 5 exhibits a pseudo-octahedral geometry where the agostic H1 center occupies one of the axial positions trans to η2-ethylene. In contrast, complex 5′ is a penta-coordinated trigonal-bipyramidal species. Key structural characteristics of the metal-ethyl moiety of 5 are Rh−C2−C1 = 82°, Rh−C2−C1−H1 = 4°, Rh−C1 = 236 pm, and Rh−H1 = 191 pm (Table S3 of the Supporting Information). In contrast, the ethyl moiety of 5′ points away from the metal center, changing these key characteristics to Rh−C2−C1 = 126°, Rh− C2−C1−H1 = 42°, Rh−C1 = 316, and Rh−H1 = 331 pm. Both resulting complexes 6 and 6′ feature square-pyramidal geometries with an agostic C−H bond of the butyl ligand at one of the square-planar corners (Table S3 of the Supporting Information). In complex 6, the C−H moiety in the δ position at the butyl ligand participates in an agostic interaction with the Rh center (Rh−C1 = 238 pm, Rh−H = 172 pm; Table S3 of the Supporting Information), whereas in 6′ the agostic interaction is due to the C−H moiety in γ-position of the butyl ligand, but showing very similar structural characteristics (Rh−C2 = 235 pm, Rh−H = 173 pm; Table S3 of the Supporting Information). As structure 6 is 21 kJ mol−1 lower in energy than 6′, the products of these two reactions 6 and 6′ are more similar in energy than the reactants 5 and 5′. Important structural parameters of the TSs TS 5−6 and 5′− 6′ are given in Figure 2. The higher activation barrier of TS 5− 6 likely reflects the stability of 5 due to the agostic interaction. Both reactions are endothermic, with reaction energies of 31 and 14 kJ mol−1, respectively (Table 1). The corresponding barriers were calculated at 145 kJ mol−1 (5 → 6) and 90 kJ mol−1 (5′ → 6′, Table 2). The RMSD of the active site structures of reactant 5, TS 5− 6, and product 6 are below 7 pm in all cluster models (Table 3). A detailed examination of the active site in complex 6 reveals that the hydrocarbon moiety (CH3)(CH2)3 of the fragment [Al(O)2Rh(C4H10)] is differently oriented relative to the Al(O)2Rh group of the pure cluster models 1T and 5T on the one hand and the embedded cluster models and the periodic model on the other hand. This can be demonstrated by inspecting the positions of the hydrocarbon moiety and the Al center relative to the plane O−Rh−O (Figure S5 of the Supporting Information). In the pure cluster models, both groups are located on the same side of the plane, whereas they are in opposite half-spaces in all models that account for the zeolite framework. Starting the optimization of the 5/83T model from the 1T configuration yields the originally calculated configuration with the moieties in question in opposite halfspaces. Hence, the pure cluster models 1T and 5T are unable to reproduce the reference structure. As in reaction 3→4, the structure variation with the sophistication of the models hints at an interaction of the metal fragment with the zeolite wall. In product 6, the shortest contacts in the periodic model are about 250 pm. No such deviating structures were determined for the reaction 5′→6′. All stationary states are similar among all types of models, periodic, pure clusters, and embedded clusters. 25083



Article

Table 4. Reaction Energies Er of Four C−C Coupling Reactions (Scheme 1) at a Zeolite-Supported Rh Center from SinglePoint Calculations on the Full 83T Clustera Using Various Exchange-Correlation Approximations with and without Dispersion Correctionb functional

1−2 Er

ΔEr

Er

ΔEr

Er

ΔEr

Er

ΔEr

5/83T,b UFF PBE PBE-D B97D ωB97x-D Periodic, PBE-D

−7 1 3 22 20 14

−21 −13 −11 8 6 −

11 17 19 29 10 10

1 7 9 19 0 −

19 36 32 33 12 31

−12 5 1 2 −19 −

9 13 6 16 −20 14

−5 −1 −8 2 −34 −

3−4

5−6

MAD

5′−6′

9.8 6.5 7.3 7.8 14.8

a Structure of the model 5/83T, determined at the ONIOM(PBE:UFF) level. bAlso shown are the differences ΔEr of the reaction energies relative to the values from the corresponding periodic models and the corresponding mean absolute deviation (MAD). All energies are in kilojoule per inverse mol.

Table 5. Activation Energies Ea (kJ mol−1) of Four C−C Coupling Reactions (Scheme 1) at a Zeolite-Supported Rh Center from Single-Point Calculations on the Full 83T Clustera Using Various Exchange-Correlation Approximations with and without Dispersion Correctionb functional 5/83T,b UFF PBE PBE-D B97D ωB97x-D periodic, PBE-D

1−2

3−4

5−6

MAD

5′-6′

Ea

ΔEa

Ea

ΔEa

Ea

ΔEa

Ea

ΔEa

102 85 84 102 95 93

9 −8 −9 9 2 −

129 128 133 149 153 142

−13 −14 −9 7 11 −

146 148 148 155 129 145

1 3 3 10 −16 −

93 96 91 98 70 90

3 6 1 8 −20 −

6.5 7.8 5.5 8.5 12.3

Structure of the model 5/83T, determined at the ONIOM(PBE:UFF) level. bAlso shown are the differences ΔEa of the activation energies relative to the values from the corresponding periodic models and the corresponding mean absolute deviation (MAD). All energies in kilojoule per inverse mol. a

various models via the MAD values of the reaction energies Er (Table 1), one notices, not surprisingly, that the embedded cluster model 11/240T performs best on average with MAD = 7.0 kJ mol−1. This is a very satisfactory result, given the notably larger computational effort required by the periodic reference. The model 5/240T performs almost as well, with MAD = 8.5 kJ mol−1. Both smaller embedding clusters (i.e., the models x/ 83T) are not much worse, with MAD ≈ 9.5 kJ mol−1. The good performance of the bare 1T model, by this criterion, seems fortuitous as a comparison with the result for the bare 5T model suggests MAD = 13.5 kJ mol−1. This latter assessment is confirmed by the overall accuracy of the activation energies, calculated with the 1T model (see below; Table 2). Finally, the very large MAD value, 25.5 kJ mol−1, for the embedding at the PM6 level (5/83/T) very clearly demonstrates that embedding by an environment at the UFF level is much superior (Table 1). The analogous analysis of the MAD values of the activation energies Ea (Table 2) yields a very similar picture. In fact, the overall quality of the calculated reaction barriers is even higher because the pertinent energies are notably larger than the reaction energies. The model 11/240T again shows the best overall performance, with MAD = 6.0 kJ mol−1. However, the other UFF embedded cluster models, even the 5T bare cluster, all fare very well, with MAD values ranging from 6.5 to 7.3 kJ mol−1. The 1T model now shows its limits, with MAD = 9.5 kJ mol−1; this value still seems reasonable if only a rough approximation is required. Just as before for reaction energies, the cluster 5/83T, embedded in a PM6 environment, is also not useful for calculating activation energies; MAD = 23.8 kJ mol−1 (Table 2).

The low RMSD, less than 9 pm, for all models corroborates this statement (Table 3). Regarding the energetics of reaction 5→6, it is not unexpected that all cluster models predict uniformly negative deviations ΔEr of −20 to −12 kJ mol−1 from the periodic model. The 11/83T model is a slight exception, not in trend, but in size, with ΔEr = −7 kJ mol−1. The deviations ΔEa of the corresponding barriers span a similarly narrow interval, from 1 to −9 kJ mol−1. For reaction 5′→6′, the cluster models describe even better both the reaction energy and the activation barrier. The ΔEr values range from −9 to −5 kJ mol−1, and the ΔEa values fall in the range of 1−4 kJ mol−1. Finally, we address some complications encountered when applying the cluster models to the reaction 5→6. With all cluster models, bare and embedded, we determined the same two isomers for the reactant 5 where the alternative structure is by 2−4 kJ mol−1 higher in energy, depending on the cluster model used, than the “standard” structure 5 presented above. The two isomers are the result of alternative relative orientations of the CH2 groups at the centers C2 and C3 (i.e., the centers which are to be bound in that reaction). Similarly, an alternative rotamer structure was also found for the TS 5−6. In contrast, that alternative TS structure was even lower in energy, by 4−12 kJ mol−1, compared to the corresponding “standard” structure discussed above. We were unable to locate these alternative rotamer structures with the periodic reference model. From this point of view, the cluster models, bare and embedded, added variants to reaction 5→6 beyond the reference model. 3.4. Overall Comparison of the Various Cluster Models. When one analyzes the overall performance of the 25084



Article

are reproduced well, within 10 kJ mol−1 (Tables 4 and 5). It is well-known that proper accounting of dispersion interactions is indispensable for describing the adsorption enthalpies in zeolites, in particular in H-ZSM-5.32−34 Speybroeck et al.34 showed that experimental and theoretical adsorption enthalpies of alcohols and nitriles, obtained by PBE, deviate on average by about 35 kJ mol−1 (mean absolute deviation − MAD); inclusion of dispersion interactions (PBE-D2) notably reduced the corresponding MAD value to ∼7 kJ mol−1.34 Unlike adsorption, all bond-breaking and bond-making reactions studied in the present work take place at the extra framework metal center Rh. Therefore, the effect of including dispersion corrections is rather small, despite a few short contacts between the hydrocarbons and the zeolite wall. The MAD values of PBE and PBE-D results differ by 0.8 kJ mol−1 for reaction energies (Table 4) and by 2.3 kJ mol−1 for activation energies (Table 5).

3.5. Comparison of Various Exchange-Correlation Functionals. To evaluate the performance of other DFT methods regarding the energetics, we carried out single-point DFT calculations on the 83T cluster model, taking the structure of the 5/83T model optimized with the ONIOM(PBE:UFF) method. The strategy of carrying out full QM single-point calculations on a QM/MM structure has been shown before to be valuable.34,47 The reaction energies and activation barriers obtained with the functionals PBE, PBE-D, B97D, and ωB97x-D are collected in Tables 4 and 5, respectively. These functionals were selected as follows. PBE serves as standard GGA, PBE-D is its dispersion-corrected extension and was used by us as reference (see Section 2), ωB97x-D was shown before to yield good results for adsorption enthalpies in zeolites32 and B97D as another functional that comprises a dispersion-correction term. Recall that in the present study we systematically compared PBE-D2 energies of the periodic reference with ONIOM(PBE:UFF) energies of the embedded cluster models. We started the present set of single-point calculations by obtaining PBE results for our test set. At this all-QM level for the 83T cluster, the mean absolute deviation (MAD) of the differences ΔEr to the periodic PBE-D reference are reduced from 9.8 to 6.5 kJ mol−1 (Table 4). Inclusion of the dispersion interaction, PBE-D2, yields only small changes further with respect to the periodic PBE-D reference, increasing the MAD by 0.8 to 7.3 kJ mol−1. From these quality criteria, one concludes that already the more economic embedded cluster calculations at the ONIOM(PBE:UFF) level overall show very satisfactory agreement with the periodic reference. The reaction energies obtained with the B97D functional are quantified on average by MAD = 7.8 kJ mol−1 (i.e., they are of comparable quality as the PBE-D results). In contrast, the functional variant ωB97x-D yields quite substantial deviations ΔEr from the periodic PBE-D reference, with MAD = 14.8 kJ mol−1; this deterioration is due to the reactions 5→6 and 5′→6′ (Table 4). Inspection of Table 5 shows that the corresponding singlepoint results for the activation energies are overall of similar quality. The MAD values of the deviations ΔEa from the reference are mostly smaller than those of the ΔEr values, although this does not seem to be a systematic trend. The allQM results at the PBE level for the 83T cluster do not agree quite well with the periodic PBE-D reference as the results at the ONIOM(PBE:UFF) level. The corresponding MAD values of the differences ΔEa are 6.5 and 7.8 kJ mol−1, respectively (Table 5). For reaction barriers, inclusion of the dispersion interaction improves the agreement with the reference, in contrast to the reaction energies. The MAD value of ΔEa, 5.5 kJ mol−1, is 2.3 kJ mol−1 smaller than the MAD value of the PBED2 results. Overall, we reach the same conclusion as for the reaction energies. The ONIOM cluster calculations, with dispersion interactions reflected by the UFF embedding, agree satisfactorily with the periodic reference. The activation energies obtained with the B97D functional are quantified on average by MAD = 8.5 kJ mol−1 (i.e., they are of similar quality as the PBE-D2 results). The functional variant ωB97x-D again yields substantial deviations ΔEa from the periodic PBE-D reference, with MAD = 12.3 kJ mol−1, once again due to the reactions 5→6 and 5′→6′ (Table 5). Overall, we obtained improved energies, calculated for the 83T cluster model, when the PBE energy was corrected for dispersion interaction (i.e., at the PBE-D level). After applying the D-correction, almost all energies, Er and Ea, of the test set

4. SUMMARY AND OUTLOOK Several variants of a catalytic elementary step, namely C−C coupling on a Rh−diethylene complex inside faujasite, were studied by cluster model strategies. We used two small pure clusters, 1T and 5T, and four cluster models, 5/83T, 11/83T, 5/240T, and 11/240T, where small QM partitions (5T, 11T) were embedded in a force-field environment, described at the UFF level. We validated these cluster model results by comparison with reference results, obtained with a threedimensional periodic model. A detailed analysis revealed that the semiempirical PM6 method for embedding of the QM partition, instead of the UFF, is notably less satisfactory for structures and energetics. The small cluster 1T is unable to reproduce the structure of complex 2. For lack of confinement effects, neither simple model, 1T and 5T, offers an adequate representation of complex 6. Recall that both pure cluster models are unable to reproduce the asymmetric arrangement between Rh and the two oxygen atoms of the zeolite. In other cases, lack of the framework gave rise to various isomers of the hydrocarbon species at the metal center (e.g., for TS 5−6). Nevertheless, these very simple models reproduce the energetics of the reactions studied within ∼10 kJ mol−1 and therefore can be employed when one needs to quickly assess reaction profiles, especially for complex systems like metal-containing zeolites. Among the embedded cluster models, the smaller model 5/83T suffices for adequately describing the reaction mechanism. Larger embedded cluster models do not yield any significant improvement in structure or energetics, neither by extending the QM partition, 5T→11T, nor the MM partition, 83T→ 240T. The largest change in energetics, obtained due to these extensions of either partition, was 8 kJ mol−1, excluding one exceptional case. Throughout our analysis, we only witnessed exceptional energetics when the periodic reference structure differs notably and rather uniformly from the structures provided by all cluster models. In the C−C coupling steps 1→2 and 3→4, where a MCp ring is produced, conformational fluxionality either in the TS or the product may trigger structural differences between the periodic reference and the embedded cluster models. In addition, some rotational flexibility during C−C bond formation in 5→6 or 5′→6 was attributed to a weaker agostic interaction in the cluster models than in the periodic model. Our results do not admit a definitive comment regarding the importance of functional corrections for dispersion interactions. 25085



Article

(3) Kubička, D.; Kubičková, I.; Č ejka, J. Application of Molecular Sieves in Transformations of Biomass and Biomass-Derived Feedstocks. Catal. Rev.: Sci. Eng. 2013, 55, 1−78. (4) Serrano, D.; Pizarro, P. Synthesis Strategies in the Search for Hierarchical Zeolites. Chem. Soc. Rev. 2013, 4004−4035. (5) Van Santen, R.; Kramer, G. J. Reactivity Theory of Zeolitic Brønsted Acidic Sites. Chem. Rev. 1995, 95, 637−660. (6) Yu, M.; Noble, R. D.; Falconer, J. L. Zeolite Membranes: Microstructure Characterization and Permeation Mechanisms. Acc. Chem. Res. 2011, 44, 1196−1206. (7) Lukyanov, D. B.; Beckett, S. J.; Vazhnova, T. Toward Better Understanding of the Catalytic Action of Acidic Zeolites: Investigation in Methane and Ethane Activation and Transformation. ACS Catal. 2012, 2, 2596−2601. (8) Mlinar, A. N.; Zimmerman, P. M.; Celik, F. E.; Head-Gordon, M.; Bell, A. T. Effects of Brønsted-Acid Site Proximity on the Oligomerization of Propene in H-MFI. J. Catal. 2012, 288, 65−73. (9) Olsbye, U.; Svelle, S.; Bjørgen, M.; Beato, P.; Janssens, T. V.; Joensen, F.; Bordiga, S.; Lillerud, K. P. Conversion of Methanol to Hydrocarbons: How Zeolite Cavity and Pore Size Controls Product Selectivity. Angew. Chem., Int. Ed. 2012, 51, 5810−5831. (10) Van der Mynsbrugge, J.; Visur, M.; Olsbye, U.; Beato, P.; Bjørgen, M.; Van Speybroeck, V.; Svelle, S. Methylation of Benzene by Methanol: Single-Site Kinetics over H-ZSM-5 and H-beta Zeolite Catalysts. J. Catal. 2012, 292, 201−212. (11) Bayram, E.; Lu, J.; Aydin, C.; Uzun, A.; Browning, N. D.; Gates, B. C.; Finke, R. G. Mononuclear Zeolite-Supported Iridium: Kinetic, Spectroscopic, Electron Microscopic, and Size-Selective Poisoning Evidence for an Atomically Dispersed True Catalyst at 22° C. ACS Catal. 2012, 2, 1947−1957. (12) Lu, J.; Martinez-Macias, C.; Aydin, C.; Browning, N.; Gates, B. C. Zeolite-Supported Bimetallic Catalyst: Controlling Selectivity of Rhodium Complexes by Nearby Iridium Complexes. Catal. Sci. Technol. 2013, 3, 2199−2203. (13) Serna, P.; Gates, B. C. A Bifunctional Mechanism for Ethene Dimerization: Catalysis by Rhodium Complexes on Zeolite HY in the Absence of Halides. Angew. Chem. 2011, 123, 5642−5645. (14) Serna, P.; Gates, B. C. Zeolite-and MgO-Supported Rhodium Complexes and Rhodium Clusters: Tuning Catalytic Properties to Control Carbon−Carbon Vs. Carbon−Hydrogen Bond Formation Reactions of Ethene in the Presence of H2. J. Catal. 2013, 308, 201− 212. (15) Serna, P.; Gates, B. C. Zeolite-Supported Rhodium Complexes and Clusters: Switching Catalytic Selectivity by Controlling Structures of Essentially Molecular Species. J. Am. Chem. Soc. 2011, 133, 4714− 4717. (16) Yardimci, D.; Serna, P.; Gates, B. C. Surface-Mediated Synthesis of Dimeric Rhodium Catalysts on MgO: Tracking Changes in the Nuclearity and Ligand Environment of the Catalytically Active Sites by X-Ray Absorption and Infrared Spectroscopies. Chem.−Eur. J. 2013, 19, 1235−1245. (17) De Moor, B. A.; Reyniers, M.-F.; Sierka, M.; Sauer, J.; Marin, G. B. Physisorption and Chemisorption of Hydrocarbons in H-FAU Using QM-Pot (MP2//B3LYP) Calculations. J. Phys. Chem. C 2008, 112, 11796−11812. (18) Hemelsoet, K.; Van der Mynsbrugge, J.; De Wispelaere, K.; Waroquier, M.; Van Speybroeck, V. Unraveling the Reaction Mechanisms Governing Methanol-to-Olefins Catalysis by Theory and Experiment. ChemPhysChem 2013, 14, 1526−1545. (19) Kletnieks, P. W.; Liang, A. J.; Craciun, R.; Ehresmann, J. O.; Marcus, D. M.; Bhirud, V. A.; Klaric, M. M.; Hayman, M. J.; Guenther, D. R.; Bagatchenko, O. P.; Dixon, D. A.; Gates, B. C.; Haw, J. F. Molecular Heterogeneous Catalysis: A Single-Site Zeolite-Supported Rhodium Complex for Acetylene Cyclotrimerization. Chem.−Eur. J. 2007, 13, 7294−7304. (20) Maihom, T.; Boekfa, B.; Sirijaraensre, J.; Nanok, T.; Probst, M.; Limtrakul, J. Reaction Mechanisms of the Methylation of Ethene with Methanol and Dimethyl Ether over H-ZSM-5: An ONIOM Study. J. Phys. Chem. C 2009, 113, 6654−6662.

For the present examples, such corrections are very small, at most 4 kJ/mol. Comparison of the cluster models with periodic PBE-D results confirmed that the embedding ONIOM scheme provides an appealing alternative for studying catalytic reactions on zeolite-supported metal complexes. The embedded cluster strategy is quite economic and can be applied without encountering technical difficulties. Further advantages of the embedded cluster strategy are a wider range of analysis tools in many quantum chemistry codes and, last but not least, the costeffective usage of hybrid functionals, exclusively applied to the active site. Another aspect of embedding models seems worth noting. The number of isomers and pathways determined with cluster models was in some cases notably larger than encountered with the periodic model. We even obtained a different isomer for one transition state (TS 3−4), where the interaction between a hydrogen atom of a substrate and a zeolite wall was not correctly represented in the embedded cluster models because this interaction occurred exactly at the boundary between highand low-level partitions. This finding reflects a potentially weak point of an embedding strategy that may not be easy to anticipate.

■

ASSOCIATED CONTENT

S Supporting Information *

Structure of faujasite and designators of unique oxygen center; energies of singlet and triplet states of all intermediates (1 to 6′), data for and discussion of the C−H bond formation reaction 3→5; structural parameters of reactants, transition structures, and products of all elementary steps under study: 1→2, 3→4, 5→6, and 5′→6′; pertinent vibrational frequencies; RMSD values of stationary structures obtained with various cluster models with respect to the model 5/83T; total energies of reactants, transition structures, and products, Cartesian coordinates of all pertinent structures are included as XMol XYZ files. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions §

S.D. and A.G. contributed equally.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Georgi N. Vayssilov for valuable discussions and suggestions. This work was partially supported by the Institute of Chemical and Engineering Sciences, Agency for Science, Technology and Research (A*STAR), Singapore. We acknowledge generous computing resources provided by the A*STAR Computational Resource Centre.

■

REFERENCES

(1) Corma, A. Inorganic Solid Acids and Their Use in Acid-Catalyzed Hydrocarbon Reactions. Chem. Rev. 1995, 95, 559−614. (2) Degnan, T. F. The Implications of the Fundamentals of Shape Selectivity for the Development of Catalysts for the Petroleum and Petrochemical Industries. J. Catal. 2003, 216, 32−46. 25086



Article

(21) Nasluzov, V. A.; Ivanova, E. A.; Shor, A. M.; Vayssilov, G. N.; Birkenheuer, U.; Rösch, N. Elastic Polarizable Environment Cluster Embedding Approach for Covalent Oxides and Zeolites Based on a Density Functional Method. J. Phys. Chem. B 2003, 107, 2228−2241. (22) Nie, X.; Janik, M. J.; Guo, X.; Song, C. Shape-Selective Methylation of 2-Methylnaphthalene with Methanol over H-ZSM-5 Zeolite: A Computational Study. J. Phys. Chem. C 2012, 116, 4071− 4082. (23) Grimme, S. Accurate Description of Van Der Waals Complexes by Density Functional Theory Including Empirical Corrections. J. Comput. Chem. 2004, 25, 1463−1473. (24) Svelle, S.; Tuma, C.; Rozanska, X.; Kerber, T.; Sauer, J. Quantum Chemical Modeling of Zeolite-Catalyzed Methylation Reactions: Toward Chemical Accuracy for Barriers. J. Am. Chem. Soc. 2009, 131, 816−825. (25) Van Speybroeck, V.; Van der Mynsbrugge, J.; Vandichel, M.; Hemelsoet, K.; Lesthaeghe, D.; Ghysels, A.; Marin, G. B.; Waroquier, M. First Principle Kinetic Studies of Zeolite-Catalyzed Methylation Reactions. J. Am. Chem. Soc. 2011, 133, 888−899. (26) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (27) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A Consistent and Accurate Ab initio Parametrization of Density Functional Dispersion Correction (DFT-D) for the 94 Elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (28) Zhao, Y.; Truhlar, D. G. Density Functionals with Broad Applicability in Chemistry. Acc. Chem. Res. 2008, 41, 157−167. (29) Zhao, Y.; Truhlar, D. G. The M06 Suite of Density Functionals for Main Group Thermochemistry, Thermochemical Kinetics, Noncovalent Interactions, Excited States, and Transition Elements: Two New Functionals and Systematic Testing of Four M06-Class Functionals and 12 Other Functionals. Theor. Chem. Acc. 2008, 120, 215−241. (30) Lee, K.; Murray, É. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy Van Der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (31) Chai, J.-D.; Head-Gordon, M. Long-Range Corrected Hybrid Density Functionals with Damped Atom−Atom Dispersion Corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615−6620. (32) Chiu, C.-c.; Vayssilov, N. G.; Genest, A.; Borgna, A.; Rösch, N. Predicting Adsorption Enthalpies on Silicalite and HZSM-5. A Benchmark Study on DFT Strategies. J. Comput. Chem. 2014, 35, 809−819. (33) Hansen, N.; Kerber, T.; Sauer, J.; Bell, A. T.; Keil, F. J. Quantum Chemical Modeling of Benzene Ethylation over H-ZSM-5 Approaching Chemical Accuracy: A Hybrid MP2:DFT Study. J. Am. Chem. Soc. 2010, 132, 11525−11538. (34) Van der Mynsbrugge, J.; Hemelsoet, K.; Vandichel, M.; Waroquier, M.; Van Speybroeck, V. Efficient Approach for the Computational Study of Alcohol and Nitrile Adsorption in H-ZSM-5. J. Phys. Chem. C 2012, 116, 5499−5508. (35) Bakowies, D.; Thiel, W. Hybrid Models for Combined Quantum Mechanical and Molecular Mechanical Approaches. J. Phys. Chem. 1996, 100, 10580−10594. (36) Boscoboinik, J. A.; Yu, X.; Yang, B.; Fischer, F. D.; Włodarczyk, R.; Sierka, M.; Shaikhutdinov, S.; Sauer, J.; Freund, H.-J. Modeling Zeolites with Metal-Supported Two-Dimensional Aluminosilicate Films. Angew. Chem., Int. Ed. 2012, 51, 6005−6008. (37) Cheng, L.; Curtiss, L. A.; Assary, R. S.; Greeley, J.; Kerber, T.; Sauer, J. Adsorption and Diffusion of Fructose in Zeolite HZSM-5: Selection of Models and Methods for Computational Studies. J. Phys. Chem. C 2011, 115, 21785−21790. (38) Chu, Y.; Han, B.; Zheng, A.; Yi, X.; Deng, F. Pore Selectivity for Olefin Protonation Reactions Confined inside Mordenite Zeolite: A Theoretical Calculation Study. J. Phys. Chem. C 2013, 117, 2194− 2202.

(39) Fu, H.; Xie, S.; Fu, A.; Ye, T. Theoretical Study of the CarbonylEne Reaction between Formaldehyde and Propylene on the Mgy Zeolite. Comp. Theor. Chem. 2012, 982, 51−57. (40) Götz, K.; Meier, F.; Gatti, C.; Burow, A. M.; Sierka, M.; Sauer, J.; Kaupp, M. Modeling Environmental Effects on Charge Density Distributions in Polar Organometallics: Validation of Embedded Cluster Models for the Methyl Lithium Crystal. J. Comput. Chem. 2010, 31, 2568−2576. (41) Lin, H.; Truhlar, D. G. QM/MM: What Have We Learned, Where Are We, and Where Do We Go from Here? Theor. Chem. Acc. 2007, 117, 185−199. (42) Lomratsiri, J.; Probst, M.; Limtrakul, J. Structure and Adsorption of a Basic Probe Molecule on H-ZSM-5 Nanostructured Zeolite: An Embedded ONIOM Study. J. Mol. Graphics Modell. 2006, 25, 219− 225. (43) Sauer, J.; Sierka, M. Combining Quantum Mechanics and Interatomic Potential Functions in Ab initio Studies of Extended Systems. J. Comput. Chem. 2000, 21, 1470−1493. (44) Shor, A. M.; Shor, E. A. I.; Laletina, S.; Nasluzov, V. A.; Vayssilov, G. N.; Rösch, N. Effect of the Size of the Quantum Region in a Hybrid Embedded-Cluster Scheme for Zeolite Systems. Chem. Phys. 2009, 363, 33−41. (45) Svensson, M.; Humbel, S.; Froese, R. D.; Matsubara, T.; Sieber, S.; Morokuma, K. ONIOM: A Multilayered Integrated MO + MM Method for Geometry Optimizations and Single Point Energy Predictions. A Test for Diels-Alder Reactions and Pt(P(t-Bu)3)2 + H2 Oxidative Addition. J. Phys. Chem. 1996, 100, 19357−19363. (46) Vayssilov, G. N.; Petrova, G. P.; Shor, E. A. I.; Nasluzov, V. A.; Shor, A. M.; Petkov, P. S.; Rösch, N. Reverse Hydrogen Spillover on and Hydrogenation of Supported Metal Clusters: Insights from Computational Model Studies. Phys. Chem. Chem. Phys. 2012, 14, 5879−5890. (47) Yadnum, S.; Choomwattana, S.; Khongpracha, P.; Sirijaraensre, J.; Limtrakul, J. Comparison of Cu-ZSM-5 Zeolites and Cu-MOF-505 Metal-Organic Frameworks as Heterogeneous Catalysts for the Mukaiyama Aldol Reaction: A DFT Mechanistic Study. ChemPhysChem 2013, 14, 923−928. (48) Zimmerman, P. M.; Head-Gordon, M.; Bell, A. T. Selection and Validation of Charge and Lennard-Jones Parameters for QM/MM Simulations of Hydrocarbon Interactions with Zeolites. J. Chem. Theory Comput. 2011, 7, 1695−1703. (49) Chu, Y.; Han, B.; Zheng, A.; Deng, F. Influence of Acid Strength and Confinement Effect on the Ethylene Dimerization Reaction over Solid Acid Catalysts: A Theoretical Calculation Study. J. Phys. Chem. C 2012, 116, 12687−12695. (50) Namuangruk, S.; Pantu, P.; Limtrakul, J. Investigation of Ethylene Dimerization over Faujasite Zeolite by the ONIOM Method. ChemPhysChem 2005, 6, 1333−1339. (51) Svelle, S.; Kolboe, S.; Swang, O. Theoretical Investigation of the Dimerization of Linear Alkenes Catalyzed by Acidic Zeolites. J. Phys. Chem. B 2004, 108, 2953−2962. (52) Liang, A. J.; Craciun, R.; Chen, M.; Kelly, T. G.; Kletnieks, P. W.; Haw, J. F.; Dixon, D. A.; Gates, B. C. Zeolite-Supported Organorhodium Fragments: Essentially Molecular Surface Chemistry Elucidated with Spectroscopy and Theory. J. Am. Chem. Soc. 2009, 131, 8460−8473. (53) Olson, D.; Dempsey, E. The Crystal Structure of the Zeolite Hydrogen Faujasite. J. Catal. 1969, 13, 221−231. (54) Kramer, G. J.; Van Santen, R. Theoretical Determination of Proton Affinity Differences in Zeolites. J. Am. Chem. Soc. 1993, 115, 2887−2897. (55) Sastre, G.; Katada, N.; Suzuki, K.; Niwa, M. Computational Study of Brønsted Acidity of Faujasite. Effect of the Al Content on the Infrared OH Stretching Frequencies. J. Phys. Chem. C 2008, 112, 19293−19301. (56) Baerlocher, C.; McCusker, L. B.; Olson, D. H. Atlas of Zeolite Framework Types; 7th ed.; Elsevier: Amsterdam, 2007; http://www. iza-structure.org/databases. 25087



Article

(57) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169−11186. (58) Blöchl, P. E. Projector Augmented-Wave Method. Phys. Rev. B 1994, 50, 17953−17979. (59) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (60) Kresse, G.; Hafner, J. Ab initio Molecular-Dynamics Simulation of the Liquid-Metal−Amorphous-Semiconductor Transition in Germanium. Phys. Rev. B 1994, 49, 14251−14269. (61) Henkelman, G.; Jónsson, H. A Dimer Method for Finding Saddle Points on High Dimensional Potential Surfaces Using Only First Derivatives. J. Chem. Phys. 1999, 111, 7010−7022. (62) Jønsson, H.; Mills, G.; Jacobsen, K. Classical and Quantum Dynamics in Condensed Phase Simulations; Berne, B. J., Ciccotti, G., Coker, D. F., Eds. World Scientific: Singapore, 1998; p 385. (63) Kresse, G.; Furthmüller, J.; Hafner, J. Ab initio Force Constant Approach to Phonon Dispersion Relations of Diamond and Graphite. Europhys. Lett. 1995, 32, 729−734. (64) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. et al. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (65) Binkley, J. S.; Pople, J. A.; Hehre, W. J. Self-Consistent Molecular Orbital Methods. 21. Small Split-Valence Basis Sets for First-Row Elements. J. Am. Chem. Soc. 1980, 102, 939−947. (66) Bergner, A.; Häußermann, U.; Dolg, M.; Stoll, H.; Preuß, H. Energy-Adjusted Ab initio Pseudopotentials for the Second and Third Row Transition Elements. Theor. Chim. Acta 1990, 77, 123−141. (67) Stewart, J. J. Optimization of Parameters for Semiempirical Methods V: Modification of NDDO Approximations and Application to 70 Elements. J. Mol. Model. 2007, 13, 1173−1213. (68) Rappé, A. K.; Casewit, C. J.; Colwell, K.; Goddard, W., III; Skiff, W. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (69) Goellner, J. F.; Gates, B. C.; Vayssilov, G. N.; Rösch, N. Structure and Bonding of a Site-Isolated Transition Metal Complex: Rhodium Dicarbonyl in Highly Dealuminated Zeolite Y. J. Am. Chem. Soc. 2000, 122, 8056−8066. (70) Laletina, S. S.; Shor, A. M.; Shor, E. A.; Nasluzov, V. A.; Rubaylo, A. I. The Study of Rhodium (I) Dicarbonyl Coordination in Faujasite Cavities by DFT-Method. J. Sib. Fed. Univ., Chem. 2008, 2, 190−199. (71) Boekfa, B.; Choomwattana, S.; Khongpracha, P.; Limtrakul, J. Effects of the Zeolite Framework on the Adsorptions and HydrogenExchange Reactions of Unsaturated Aliphatic, Aromatic, and Heterocyclic Compounds in ZSM-5 Zeolite: A Combination of Perturbation Theory (MP2) and a Newly Developed Density Functional Theory (M06-2X) in ONIOM Scheme. Langmuir 2009, 25, 12990−12999. (72) Namuangruk, S.; Tantanak, D.; Limtrakul, J. Application of ONIOM Calculations in the Study of the Effect of the Zeolite Framework on the Adsorption of Alkenes to ZSM-5. J. Mol. Catal. A: Chem. 2006, 256, 113−121. (73) Kubas, G. J. Metal−Dihydrogen and σ-Bond Coordination: The Consummate Extension of the Dewar−Chatt−Duncanson Model for Metal−Olefin π Bonding. J. Organomet. Chem. 2001, 635, 37−68. (74) Macgregor, S.; Whittlesey, M.; Perutz, R. A Theoretical Study of [M(Ph3)4] (M = Ru or Fe), Models for the Highly Reactive d8 Intermediates [M(dmpe)2] (dmpe = Me2PCH2CH2PMe2). Zero Activation Energies for Addition of CO and Oxidative Addition of H2. J. Chem. Soc., Dalton Trans. 1998, 291−300.

25088


Modeling Catalytic Steps on Extra-Framework Metal Centers in

Recommend Documents