11796
J. Phys. Chem. C 2008, 112, 11796–11812
Physisorption and Chemisorption of Hydrocarbons in H-FAU Using QM-Pot(MP2//B3LYP) Calculations Bart A. De Moor,† Marie-Franc¸oise Reyniers,*,† Marek Sierka,‡ Joachim Sauer,‡ and Guy B. Marin† Laboratorium Voor Chemische Technologie, Ghent UniVersity, Krijgslaan 281 S5, B-9000 Gent, Belgium, and Institut fu¨r Chemie, Humboldt-UniVersita¨t zu Berlin, Unter den Linden 6, 10099 Berlin, Germany ReceiVed: NoVember 22, 2007; ReVised Manuscript ReceiVed: April 25, 2008
The DFT parametrized zeolite force field in the QM-Pot program is extended with carbon-carbon, carbon-hydrogen, and alkoxy bond describing parameters. The extended force field has been combined with B3LYP and with MP2 as the high-level quantum mechanical (QM) method to simulate the physisorption and chemisorption of ethene, isobutene, 1-butene, 1-pentene, and 1-octene in H-FAU (Si/AlF ) 95) and for physisorption of 1-pentene, n-pentane, 1-octene, and n-octane in all silica FAU. The new parametrization predicts more stable chemisorption complexes than physisorbed π complexes, but with smaller chemisorption energies which are more reliable as shown by comparison with experimental results and with accurate hybrid MP2:DFT calculations. An embedded cluster size study shows that, due to the importance of the stabilizing van der Waals part in the MM contribution of the cluster, QM-Pot(MP2//B3LYP) calculations yield more reliable physisorption and chemisorption energies of hydrocarbons in zeolites than QM-Pot(B3LYP). The QM-Pot(MP2//B3LYP) results are in good agreement with available experimental data. In H-FAU, the H+ · · · alkane interaction was found to contribute at most 7 kJ/mol to the total physisorption energy of n-alkanes while the H+ · · · π interaction contributes 20-25 kJ/mol to the total physisorption energy of alkenes. For n-alkene physisorption in H-FAU, a linear increase of both the physisorption and chemisorption energies of 8.7 kJ/mol per C-atom is found. The protonation energy of n-alkenes in H-FAU was found to be independent of the C-number and amounts to -50 kJ/mol for the formation of secondary alkoxides. The formation of tertiary alkoxides in H-FAU suffers slightly from steric constraints imposed by the zeolite framework. 1. Introduction Zeolites constitute a unique class of oxide materials build from corner sharing SiO4 and AlO4 tetrahedra. Due to their microporous structure and acidity, zeolites are widely used for hydrocarbon conversion reactions such as cracking, isomerization, and alkylation of aromatics.1,2 Knowledge of the molecular details of the various steps in zeolite-catalyzed reactions is rather limited, as even a simple catalytic reaction in zeolites consists of several elementary steps. However, the ever increasing need to improve the efficiency of chemical processes by optimizing the use of existing catalysts and/or the development of new catalysts requires insight into the molecular details of all the reaction steps that occur on the catalyst surface. To gain insight into the molecular behavior of reactants in zeolites, molecular simulations have been proven to be a very useful tool. Molecular simulation of catalytic reactions in zeolites is still a challenge to computational chemistry. Besides the structural complexity of the zeolite, atomistic modeling of reactions in zeolites requires an accurate description not only of the bond breaking/forming process but also of the nonbonded interactions (e.g., hydrogen bonding and van der Waals interactions) between the hydrocarbon and the zeolite. The computational simplicity of empirical force fields allows to study up to 106 atoms and long time scales (up to nanoseconds). State-of-the-art force fields are extremely well suited to describe the relaxation of the zeolite * Corresponding author. Tel.: +32 9 264 5677. Fax: +32 9 264 5824. E-mail:
[email protected]. † Ghent University. ‡ Humboldt-Universita ¨ t zu Berlin.
frame and of the nonbonded interactions between the hydrocarbon molecules and the zeolite atoms.3 They fail however to describe chemical reactions. Quantum chemical (QM) calculations are required to describe the processes of bond breaking and bond formation that occur at the active site. Traditionally, two different structural models have been used to represent the zeolite in QM calculations. In the cluster approach, the local environment of the active site is modeled by a small set of atoms representing a fragment of the zeolite. The major drawback of this approach is that small clusters used to model the active site are not typical of any particular zeolite framework and are therefore not able to explain different catalytic behavior exhibited by structurally different zeolites. In periodic density functional theory (DFT) calculations, where the zeolite is represented as an infinite perfect crystal, the effects of the zeolite framework can be accounted for in the calculations. However, periodic DFT calculations are not only computationally very demanding but also fail to yield reliable energies for the van der Waals interactions that dominate the adsorption-desorption steps, especially for long chain hydrocarbons.4–6 To overcome this shortcoming of periodic DFT calculations, Hafner et al.5 and van Santen et al.6 proposed to correct the periodic DFTresults for van der Waals interactions using an add-on LennardJones potential acting between the hydrocarbon atoms and the zeolite oxygen, silicon, and aluminum atoms. While the latter approach uses fixed DFT structures, the DFT plus damped dispersion approach recently parametrized for many atoms and a variety of functionals offers the possibility to optimize structures with inclusion of dispersion.7,8 Tuma and Sauer9,10
10.1021/jp711109m CCC: $40.75 2008 American Chemical Society Published on Web 07/11/2008
Physisorption/Chemisorption of Hydrocarbons in H-FAU
Figure 1. Schematic representation of the physisorption and chemisorption complex on a 3T embedded cluster.
suggested using hybrid MP2:DFT calculations, but these are computationally extremely demanding. Also, the use of a hybrid DFT/force field approach has been suggested.11,12 In this approach, the DFT part describing the bond rearrangement at the active site is intentionally limited to a small part of the zeolite, typically a 3T cluster, while the van der Waals interactions with the extended zeolite frame are described using a computationally much less expensive force field. In this work, the QMPOT program developed by Sauer and co-workers13,14 is used to study physisorption and chemisorption of hydrocarbons in H-FAU. Physisorption and chemisorption of hydrocarbons are the first elementary steps in most acid catalyzed hydrocarbon conversion processes and are extensively investigated both experimentally15–22 and theoretically.10–12,23–41 Physisorption of alkanes in zeolites is characterized by interactions between the hydrocarbon atoms and the highly polarizable oxygen atoms of the zeolite wall resulting from dispersive van der Waals forces and can amount to about 10-15 kJ/mol per carbon atom.42 A larger hydrocarbon molecule can, depending on the dimensions and the topological characteristics of the pores and cages, be stabilized by many zeolite oxygen atoms, and hence the total interaction energy can become important. This type of interactions gives rise to the substantial heat of physisorption of hydrocarbons in zeolites. Physisorption of alkenes in zeolites is characterized by the additional formation of a π-complex resulting from the interaction between the double bond of the alkene and the acid proton of the zeolite as depicted in Figure 1a. Protonation of physisorbed alkenes yields a chemisorbed intermediate, at 0 K found to be an alkoxide by the formation of a carbonsoxygen σ-bond as shown in Figure 1b. In the literature, the debate on whether the most stable form of the protonated intermediate at experimentally relevant temperatures is an alkoxide or a carbenium ion is still going on.10,41,43 The QM-Pot program uses the quantum mechanicssinteratomic potential functions approach and the periodic molecular mechanics (MM) calculations are done using a DFT-parametrized core-shell model for the zeolite.44 The core-shell representation for the oxygen atoms accounts for polarization. From the study of zeolite deprotonation energies,45–47 zeolite proton mobilities,48 and ammonia adsorption49,50 in zeolites, Sauer and co-workers showed that the QM-Pot approach yields very accurate results, independent of the size of the embedded cluster. Generally, the variations with the size of the embedded cluster are smaller than 7 kJ/mol. A limited number of studies on hydrocarbon conversion reactions using the QM-Pot program have been reported.12,26,27 In these studies, DFT/B3LYP calculations were performed for a small cluster to ensure that the van der Waals interaction between the hydrocarbon and the zeolite was described by the force field, while the zeolite force field was extended with Lennard-Jones terms taken from literature51 and Coulomb charges were derived from fits to the QM
J. Phys. Chem. C, Vol. 112, No. 31, 2008 11797 electrostatic potential (ESP). Lennard-Jones interactions were considered to act between the hydrocarbon atoms and the zeolite oxygen, silicon, and aluminum atoms. A good agreement with periodic DFT calculations has been observed for reactions in which dispersion energies are less important. In the study of physisorption and chemisorption of C3-C5 alkenes in H-FER, the chemisorbed σ complexes have been predicted to be much more stable than the physisorbed π complexes.12 However, the force field used in these studies lacked an adequate parametrization for the critical zeolite-alkoxide structures. In the first part of this paper, a new parametrization, in particular, for the zeolite-alkoxide bond, is presented. The new parametrization still predicts more stable chemisorption complexes than physisorbed π-complexes, but with smaller chemisorption energies which are more reliable as shown by comparison with experimental results and with accurate hybrid MP2: DFT calculations. It is demonstrated that inclusion of the zeolite-alkoxide bond potentials for the MM calculations is necessary to obtain reliable chemisorption energies for hydrocarbons in zeolites. It is also found that hybrid QM-Pot(MP2// B3LYP) calculations yield more reliable results than QMPot(B3LYP) calculations. In the second part of this paper, a QM-Pot(MP2//B3LYP) study of physisorption and chemisorption of various alkanes and alkenes in H-FAU (Si/AlF ) 95) is presented. The contribution of the H+ · · · hydrocarbon interaction to the total physisorption energy of hydrocarbons is evaluated. The relative stabilities of π and σ complexes and of primary, secondary, and tertiary alkoxides are discussed. The calculated data are compared with available experimental data. 2. H-FAU Structural Model and Computational Procedures 2.1. H-FAU Structure. The starting geometry of the faujasite structure has been taken from the database of the Structure Commission of the International Zeolite Association (IZA).52 The large pore H-FAU framework is composed of sodalite cages, which are connected by double 6-membered rings, giving rise to the characteristic supercages of faujasite, interconnected through 12-membered ring pores. The noncubic FAU unit cell (144 atoms) is doubled in the [100] direction to ensure sufficient separation between the periodic images of the hydrocarbon. The numbering of the T sites and the oxygen atoms corresponds to Hriljac et al.53 One silicon atom was replaced by an aluminum atom and the resulting negative charge is compensated for by adding a hydrogen atom at the O1 position, i.e., at the most likely and preferred neighboring oxygen atom of the aluminum,44 yielding a unit cell composition of HAlSi95O192 and Si/ AlF ) 95. A constant pressure optimization of the unit cell is performed using GULP applying the shell-model ion-pair potential zeolite force field. The resulting unit cell parameters are a ) 3497.1 pm, b ) 1743.0 pm, c ) 1749.3 pm, R ) 59.995°, β ) 59.870°, and γ ) 59.898°, in good agreement with the results obtained by Clark et al.26 using the same optimization procedure. The resulting optimized unit cell has been used in the QM-Pot calculations, not allowing any change in the unit cell parameters when optimizing hydrocarbon species inside the zeolite structure. Optimization of the noncubic FAU unit cell (144 atoms) with VASP, not doubled in the [100] direction, leads to the unit cell parameters a ) 1740.1 pm, b ) 1736.2 pm, c ) 1742.4 pm, R ) 59.869°, β ) 59.818°, and γ ) 59.879°. 2.2. Computational Procedures. 2.2.1. QM-Pot Calculations. In this study, the hydrocarbon species inside the zeolite pores have been optimized using a combined quantum
11798 J. Phys. Chem. C, Vol. 112, No. 31, 2008
De Moor et al.
∆EQM-Pot ) ∆ELR + ∆E(C)QM
Figure 2. Energy levels corresponding to the different interactions between the gas-phase hydrocarbon and the zeolite.
mechanicssinteratomic potential functions approach, using periodic boundary conditions and treating the entire zeolite structure. Calculations are performed with the QMPOT program,13,14 which couples TURBOMOLE54,55 for the QM calculations and GULP56,57 for the force field calculations. The QM-Pot energy (EQM-Pot) of the periodic system S is obtained by the following subtraction scheme:
EQM-Pot ) E(S)Pot + (E(C)QM - E(C)Pot)
(1)
The QM-Pot energy of the periodic system S is obtained from the interatomic potential functions energy of the total system S, E(S)Pot, by adding a correction that is the difference between the QM energy of the cluster, E(C)QM, and the interatomic potential energy of the cluster C, i.e., the active site and the hydrocarbon, E(C)Pot. From this definition of the QM-Pot energy, consistent expressions are derived for combined gradients and second derivatives.14 The physisorption energy, ∆Ephys, is referred to the gas-phase hydrocarbon and is defined as
∆Ephys ) Eπ - Ezeolite - Ehydrocarbon
(2)
The chemisorption energy of an alkene, ∆Echem, is also referred to the gas-phase hydrocarbon and is defined as
∆Echem ) Eσ - Ezeolite - Ehydrocarbon
(3)
The protonation energy of an alkene, ∆Eprot, is defined as the energy difference between the chemisorbed σ complex and the physisorbed π complex, i.e., the reaction energy for the alkoxide formation from physisorbed alkenes
∆Eprot ) Eσ - Eπ
(4)
These definitions are illustrated in Figure 2. The long-range MM contribution, ∆ELR, to the QM-Pot physisorption, chemisorption, and protonation energies is defined as
∆ELR ) ∆E(S)Pot - ∆E(C)Pot
(5)
The stabilizing van der Waals part in the long-range MM contribution is referred to as EvdW,LR. The stabilizing van der Waals part in the cluster-hydrocarbon MM contribution is referred to as EvdW,C. The sum of both yields the van der Waals part in the total system MM contribution is Evdw,S. Using these definitions, the QM-Pot sorption energies can be rewritten as the sum of a low-level long-range MM contribution, ∆ELR, and a high-level QM contribution, ∆E(C)QM
(6)
The embedded clusters used in this study are presented in Figure 3. All optimizations of the physisorption and chemisorption intermediates have been performed using DFT for the high-level calculations. Preliminary calculations on ethene physisorption and chemisorption in H-FAU have shown that no significant structural differences were obtained when using 3T, 8T, and 14T embedded clusters during optimization. Therefore, all optimizations of the physisorption and chemisorption intermediates have been performed on a 3T cluster embedded in H-FAU. The B3LYP functional58,59 has been used in combination with the DZP basis set on the H, C, Si, and Al atoms and the TZP basis set on the O atoms, referred to as the T(O)DZP basis set.60,61 The hydrogen link atoms lay on the broken O-Si bond, and the (Si)O-H and (Al)O-H termination distances are kept fixed at 96.66 and 96.28 pm, respectively. DFT calculations on the larger embedded clusterss8T and 14Tsare single point energy calculations on the optimized structure obtained on the 3T embedded cluster. It is known however that long-range dispersion interactions, which play a major role in the stabilization of hydrocarbons in zeolites, are problematic for DFT methods.4 As the size of the embedded cluster increases, it can be expected that, along with the van der Waals contributions being accounted for by means of the Lennard-Jones potential in the force field, van der Waals contributions in the embedded cluster itself will also become important. Therefore, single point MP2 energy calculations have been performed for the different embedded clusters. The resolution-of-the-identity (RI) approximation62 as implemented in the TURBOMOLE package has been used to obtain the MP2 energies. MP2 calculations have been performed employing Ahlrichs’ optimized triple-ζ valence basis set combined with a single set of polarization functions with exponents 0.8, 0.8, 1.2, 0.3, and 0.35 for the hydrogen, carbon, oxygen, aluminum, and silicon atoms, respectively. This basis set is referred to as the TZVP basis set.61 The total QM-Pot energies obtained by single point energy calculations of the embedded cluster at B3LYP/ T(O)DZP and MP2/TZVP levels of theory are referred to as QM-Pot(B3LYP) and QM-Pot(MP2//B3LYP), respectively. Force Fields Used in QM-Pot Calculations. In this section, the inclusion of carbon-carbon, carbon-hydrogen, and alkoxide bond related potentials in the DFT-parametrized zeolite force field described by Sierka and Sauer44 is presented. The hydrocarbon and alkoxide related bonds and angles are described using Morse and Three-Body potentials
EMorse ) De[(1 - exp(-a(r - r0)))2 - 1]
(7)
1 EThree-body ) k(θ - θ0)2 2
(8)
The Morse potential is able to describe attraction as well as repulsion between atoms using three parameterssDe, a, and r0sfor each bond type taken into account. For the Three-Body potentials, the equilibrium angle θ0 is kept fixed at a value of 109.47° or 120.00°, depending on the hybridization of the atom involved, and hence, only the estimation of the force constant k is necessary. To allow good flexibility of the force field, three different types of carbon atoms have been considered: sp2, sp3, and aromatic carbon atoms, respectively, denoted as Csp2, Csp3, and Carom. The alkoxide oxygen atom is treated analogously as the oxygen of a hydroxyl bridging group of the zeolite, and one type of hydrogen atom is considered in addition to the two types (i.e., the acid proton and the terminating hydrogen atom) already defined in the zeolite force field. In total, 5 Morse
Physisorption/Chemisorption of Hydrocarbons in H-FAU
J. Phys. Chem. C, Vol. 112, No. 31, 2008 11799
Figure 3. The 3T, 8T, and 14T embedded cluster models (without the H-link atoms), used in the cluster size study of the physisorption and chemisorption of hydrocarbons in H-FAU. (silicon, blue; aluminum, pink; oxygen, red; hydrogen, white).
potentials and 9 Three-Body potentials have been included in the force field requiring the estimation of 24 parameters. The fitting facility within the GULP program has been used to estimate the parameters of the interatomic potentials by fitting to ab initio data (the “observables”). The optimum parameters are obtained by minimizing the sum of squares, F
F)
∑
(fcalc - fobs)2
(9)
all observables
where fcalc and fobs are the calculated (by the force field) and observed (ab initio) quantities. GULP uses a Newton-Raphson functional minimization approach. A simultaneous fittings optimization procedure has been used in which not only the parameters of the potentials are adjusted, but also the shell positions of the oxygen atoms and the atomic positions of the hydrocarbon atoms are adjusted. Here, the observables are the ab initio gradients, i.e., the first derivative of the energy with respect to the atomic coordinates. Fitting to ab initio energies was not performed, as it was not the aim of this work to develop a stand-alone force field and inclusion of the energies as observables in the fitting procedure would significantly increase the complexity of the parameter estimation. All calculations have been performed with TURBOMOLE at the B3LYP/T(O)DZP level on a database consisting of 129 structures. Three types of structures were considered in the database (A) gas-phase hydrocarbon structures of methane, ethene, propene, isobutene, 2-pentene, 3-methylpentane, 2-propylbenzene, and 3-methyl-1,3-hexadiene (95 structures). (B) physisorption complexes of isobutene, 1-pentene, 1-octene, and 2-octene on a 14T cluster with SiO-H termination (19 structures). (C) chemisorption complexes of ethoxy, 1-butoxy, 1-pentoxy, 2-pentoxy, and t-butoxy on a 14T cluster with SiO-H termination (15 structures). The database contains optimized as well as distorted structures, ensuring a wide applicability of the obtained interatomic potentials. The majority of the database structures are gas-phase hydrocarbons. There are several reasons for this: (1) It is computationally less demanding to do calculations on smaller structures. (2) All potentials except the C-O related ones can be obtained from the gas-phase hydrocarbons. (3) The zeolite related potentials are retained, i.e., all the zeolite related parameters are kept fixed during the fitting procedure. Next to these carbon-carbon, carbon-hydrogen, and alkoxide bond related potentials, other forces that act on the hydrocarbon atoms result from Coulomb interactions (between the hydrocarbon atoms and between the hydrocarbon atoms and the zeolite atoms) and Lennard-Jones interactions (hydrocarbon-zeolite). Since these potentials may have an influence
on the estimated parameters, it is important to include sufficient physisorption (B) and chemisorption structures (C). The Lennard-Jones terms are considered to act between the C/H hydrocarbon atoms and the Si/Al/O zeolite atoms and were taken from the consistent valence force field (CVFF).51 These parameters have been shown to yield reliable adsorption energies.12 The 6-12 Lennard-Jones potential is applied with a cutoff radius of 40 Å. To avoid interference with the C-O bond describing Morse potential, the lower limit for the C-O Lennard-Jones potential is set at 2 Å. A complete list of all potentials added to the zeolite force field together with their cutoffs can be found as Supporting Information. The choice of the atomic point charges used to determine the Coulomb interactions requires special attention. The charges on the zeolite atoms are determined by the core-shell parameters, whereas for the charges on the hydrocarbon fragment, no unambiguous charge scheme has been proposed up to now. In previous QM-Pot studies of hydrocarbon reactions in zeolites, charges based on an ESP fit were used.12,26,27 It should however be mentioned that use of the core-shell zeolite force field in the MM calculations imposes the condition that the total sum of charges on the hydrocarbon part equals 0 in the case of physisorption and +1 in the case of chemisorption. Especially for the chemisorption complexes, a quite large adjustment of the charges derived from an ESP fit is therefore necessary, since they usually add up to values around only +0.4. The necessity to do ESP fits for all the different hydrocarbon complexes and the a posteriori correction of the ESP-derived charges required to be compatible with the zeolite force field complicates their use in QMPot calculations. Therefore, in this work we propose the use of an a priori fixed charge distribution scheme for the hydrocarbon atoms. In the case of physisorption, all hydrocarbon hydrogen atoms are assigned the charge +0.1 and each carbon atom is assigned the charge -0.1 multiplied by its number of hydrogen atoms. This scheme ensures that the formal charge on the hydrocarbon is zero. For the chemisorption complexes on the other hand, the total charge on the alkoxide fragment should amount to +1.0 to obey the unit cell charge neutrality condition, and a zero charge is assigned to all hydrocarbon atoms except for the alkoxide carbon atom, which is assigned a charge of +1.0. Results obtained from QM-Pot calculations using the reparameterized force field for the MM calculations were compared with those obtained using the force field described by Nieminen et al.12 in the MM calculations. These authors also used the DFT-parametrized shell model for the zeolite as developed by Sierka and Sauer44 with the same Lennard-Jones potentials as those used in the extended force field and with 40 Å cutoffs for
11800 J. Phys. Chem. C, Vol. 112, No. 31, 2008
De Moor et al.
TABLE 1: Point Charges on the Hydrocarbon Atoms Used in the MM Calculations QM-Pot point charge model atom
ESP charges
ESP deriveda
extended FFb
-0.32 0.16 0.00
-0.20 0.10 0.00
-0.30 0.60 0.20 0.05 +1.00
0.00 1.00 0.00 0.00 +1.00
Ethene C (dCH2) H (dCH2) total
-0.27/-0.32 0.14-0.23 0.09
C (-CH3) C (O-C) H (-C-O) H (-CH3) total
-0.40 0.30 0.02-0.04 0.11-0.15 0.36
C (dCH2) C (dC
