- Email: [email protected]
planewave-DFT scheme for large chemical systems: proton jumps in zeolites
Chemical Physics Letters 387 (2004) 388–394 www.elsevier.com/locate/cplett
A hybrid MP2/planewave-DFT scheme for large chemical systems: proton jumps in zeolites Christian Tuma, Joachim Sauer
*
Humboldt-Universit€at zu Berlin, Institut f€ur Chemie, Unter den Linden 6, Berlin D-10099, Germany Received 20 January 2004; in final form 17 February 2004 Published online: 10 March 2004
Abstract We present an embedding scheme to introduce local corrections at post Hartree–Fock level to density functional theory (DFT) calculations. As a first application we study proton jump reactions in the zeolite HSSZ-13 and show that energy barriers and rate constants are significantly changed by second-order Møller–Plesset perturbation theory (MP2) corrections to plane wave based DFT calculations. Electronic energy barriers increase from 68 to 81 kJ/mol (dry zeolite), and from 22 to 30 kJ/mol (hydrated zeolite). The predicted heats of adsorption of one water molecule onto the Brønsted acidic sites O1 and O2 are 73 and 69 kJ/mol, respectively. Ó 2004 Elsevier B.V. All rights reserved.
1. Introduction In electronic structure studies of materials and molecular systems the use of density functional theory (DFT) has met with great success. The main advantage of DFT is its computational efficiency and scaling behaviour which, especially in combination with a plane wave basis set, allows the application to extended molecular systems or to solids. For chemical problems (in particular involving transition states) hybrid density functionals are often the best choice among available functionals, see e.g., [1–3]. Unfortunately, the calculation of the exact exchange part is computationally too expensive with plane waves (e.g., [4]), and one is limited to various gradient-corrected functionals. Although the plane wave basis set quality can be improved systematically by increasing the kinetic energy cutoff, DFT itself does not offer methodological hierarchies of increased accuracy as post Hartree–Fock (HF) methods [5,6] do. The motivation of our work is the desire to combine the computational efficiency of plane wave based DFT *
Corresponding author. Fax: +493020937136. E-mail address: [email protected] (J. Sauer). URL: http://www.chemie.hu-berlin.de/ag_sauer.
0009-2614/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.02.056
with the accuracy and systematic convergence of post HF quantum mechanical (QM) methods. We present a mechanical embedding technique for chemical systems with an important fragment, e.g., a reaction site. We use DFT in a plane wave approach as a ‘low’-level method and add local corrections at post HF level using secondorder Møller–Plesset perturbation theory (MP2). This hybrid QM/QM approach extends the pool of existing high-level/low-level combinations [7–9]. As a first application we study a proton jump reaction in the zeolite HSSZ-13 (the acidic form of the zeolite chabasite), either loaded with a single water molecule per Brønsted site or being completely dry. This type of reaction has already been studied using different approaches [10–17]. Without inclusion of exact exchange, DFT calculations underestimate reaction barriers for proton jump reactions. For example, using a 3T cluster model, Auerbach et al. [12] report 72.8, 90.2 and 99.1 kJ/ mol for the proton jump barrier in the zeolite H–Y when BLYP, B3LYP, and MP2, respectively, are used. Ryder et al. [15] investigated a proton jump in the zeolite HZSM-5 employing a 5T cluster model and found energy barriers of 85.8, 108.4 and 135.1 kJ/mol (BPW91, B3LYP, BH&HLYP, respectively). For water assisted proton jumps the cluster model study by Krossner and Sauer [10] yielded 8.1 and 16.4 kJ/mol (BP86, MP2,
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
respectively) for the energy barrier. These examples illustrate the necessity to introduce at least local improvements over GGA type density functionals in studies of proton jump barriers.
2. Method We follow the mechanical embedding scheme described in [13,18]. The total system S (e.g., unit cell or large molecule) is divided into an inner part containing the active site and an outer part. When bonds are cut on definition of the inner part, link atoms L are added to saturate the dangling bonds of the inner part. The inner part and the link atoms together form the cluster model C. The total system S and the cluster model C are calculated at ‘low’-level, in this study DFT with a plane wave basis set and atomic pseudopotentials. The cluster model C is also calculated at ‘high’-level, here MP2 with standard gaussian basis sets. The new energy hypersurface is defined by a subtraction scheme as EðS; LÞ ¼ EðCÞhigh þ ½EðSÞlow EðCÞlow ;
ð1Þ
or equivalently EðS; LÞ ¼ EðSÞlow þ ½EðCÞhigh EðCÞlow :
ð2Þ
Our embedding scheme imposes constraints on the positions of the link atoms L. This implies additional contributions to forces on atoms of broken bonds which ensures translational and rotational invariance of the derivatives of EðS; LÞ. The present embedding scheme differs from the IMOMO method of Morokuma et al. [7] and other hybrid approaches [8,9] by details of the implementation of constraints. The term in squared brackets in Eq. (1) corresponds to the ‘low’-level long-range correction to the ‘high’level cluster result. It accounts for both long-range electronic effects and structural constraints imposed by the zeolite framework on the reaction site. Small cluster models without embedding miss such corrections and therefore neglect specific characteristics of the framework and its potential impact on the problem studied. Another view is given in Eq. (2). The term in squared brackets can be considered as the local ‘high’-level correction (here MP2) to the full periodic ‘low’-level (in this study DFT) result. Our mechanical embedding scheme is implemented in the QMPOT code [13]. For the present QM/QM coupling QMPOT has been extended with an additional interface to the CPMD code [19]. Along with its interface to the Turbomole programme package [20–22] we use QMPOT for combined MP2/DFT calculations.
389
3. Computational details 3.1. DFT The GGA by Perdew, Burke, and Ernzerhof (PBE) [23] is employed in connection with a computational efficient Pade form for exchange and correlation [24]. For the silicon, aluminium, oxygen, and hydrogen atoms of the zeolite system norm-conserving pseudopotentials are used (Troullier and Martins [25]) with core radii of 1.4, 1.6, 1.15, and 0.7 a0 , respectively. The kinetic energy cutoff of the plane wave basis set is 70 Ry. This value is commonly used in DFT studies of aqueous systems with this type of pseudopotential. A series of extensive test calculations showed that the presence of oxygen requires a higher cutoff of about 100 Ry for well converged DFT results. For this reason we also perform calculations at 100 Ry cutoff. All DFT plane wave calculations use the CPMD code [19]. The unit cell parameters of HSSZ-13 (periodic boundary conditions, stoichiometry AlHSi11 O24 per unit cell) have been obtained from ion-pair shell-model potential calculations with the hydrogen atom at its energetically favoured O1 position: a ¼ 945:3 pm, b ¼ 938:7 pm, c ¼ 944:1 pm, a ¼ 94:0°, b ¼ 95:4°, and c ¼ 94:8°, see [13] for details. In our calculations these cell parameters are kept constant. The cluster model (described in Section 3.2) is treated as an isolated system, i.e., without periodic boundary conditions. To enclose the electrons of this molecular system virtually completely, a computational box bigger than the unit cell of HSSZ-13 is required. This means that the basis sets for the periodic system and the cluster model include a different number of plane waves. However, test calculations on isolated molecules show that the quality of the plane wave basis is determined by its kinetic energy cutoff and not by the number of plane waves. For the cluster model convergence in the total energy is achieved for computational boxes bigger than Under this condition total energy fluctuations are 12 A. about 1.5 and 0.2 mEh at 70 and 100 Ry cutoff, respectively, which can be neglected in the context of our embedding scheme. We use a computational box of 14 A for the water molefor the cluster model, and of 10 A cule. 3.2. MP2 The cluster model for the ‘high’-level calculations consists of three tetrahedral Si[O1=2 ]4 units (3T). It is highlighted using a ball-and-stick presentation in Fig. 1. Eight hydrogen link atoms HL have been added to saturate dangling bonds of the inner part. The HL –OT bond lengths have been fixed at 95.3 pm (T ¼ Al) and 95.4 pm (T ¼ Si), mean values from bond length optimisations on a 2T cluster model. Ahlrichs’ fully
390
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
Fig. 1. (a–f) Proton jump reaction between oxygen sites O1 and O2 in the zeolite HSSZ-13 without (a–c) and in the presence of a water molecule (d–f). Structures for initial and final states (a,c,d,f) and transition state structures (b,e) are shown, the ball-and-stick parts represent the 3T cluster model for the MP2/DFT mechanical embedding calculations (see text for details).
optimised triple-f basis set is used for oxygen and double-f basis sets for the remaining elements [26]. Polarisation functions with exponents 1.2 (O), 0.8 (H), 0.35 (Si), and 0.3 (Al) have been added. In the MP2 correlation scheme all electrons are included. This basis set can be regarded as a small one for MP2 calculations. Higher angular momentum functions are necessary to achieve reasonable MP2 results, and diffuse functions are needed to reduce the basis set superposition error (BSSE). Therefore, we made a systematic MP2 study applying Dunning’s correlation consistent (aug-)ccpVXZ basis sets [27–29] to a 2T + H2 O cluster model. Single point MP2 calculations for the water binding energy and corresponding counterpoise corrections to
account for the BSSE [30] are made using a RI (resolution of the identity) implementation [22]. Electrons in the 1s orbitals of Al, Si, and O are not explicitly correlated (‘frozen core’). The results (Table 1) show that the basis set for the oxygen atoms is decisive for the BSSE, i.e., at least triple-f quality and one set of diffuse functions are necessary. Table 1 also includes the MP2 correlation energy parts of the binding energies corrected cpc for the BSSE, DEcorr: . The results for the large basis sets show that most of the BSSE (up to 90%) arises from the cpc correlation energy. Taking the DEcorr: values for different cpc basis sets, DEcorr: ðX Þ, X ¼ 2; 3; 4 (‘D;T;Q’), we carried out 2-point extrapolations to the complete basis set 1 limit, DEcorr: [31]:
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
391
Table 1 Single point RI-MP2 calculations for water adsorption on a 2T cluster model using different basis sets n
H, Al, Si
O
DE
DEcpc
DEcorr:
cpc DEcorr:
f
ðnÞa
1 DEcorr:
1 2 3 4 5 6 7 8 9 10 11 12
cc-pVDZ
cc-pVDZ cc-pVTZ cc-pVQZ aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVQZ cc-pVTZ cc-pVQZ aug-cc-pVTZ aug-cc-pVQZ cc-pVQZ aug-cc-pVQZ
88.0 72.6 67.6 64.4 66.3 67.1 74.5 68.7 65.3 65.5 69.7 65.9
53.1 57.5 59.5 54.6 58.9 61.3 58.5 60.6 59.9 61.6 61.3 61.8
21.6 20.5 20.1 17.7 20.0 21.3 22.9 21.1 20.1 20.7 22.5 21.2
5.9 12.6 15.5 11.4 15.3 17.1 13.7 16.0 15.8 17.2 16.7 17.5
0.45 0.52 0.57 0.64 0.64 0.72 0.58 0.63 0.80 0.90 0.69 0.90
(1) (2) (4) (5) (7) (9)
17.0 18.5 17.6 18.6 18.8 18.7
cc-pVTZ
cc-pVQZ
Water binding energies and corresponding BSSE-corrected values, DE and DEcpc , as well as their RI-MP2 correlation energy parts, DEcorr: and cpc 1 , are given. The fraction of the RI-MP2 correlation based BSSE on the full BSSE, f , is shown. DEcorr: is the 2-point extrapolated value for the DEcorr: complete basis set limit of the RI-MP2 correlation part of the water binding energy (all energies in kJ/mol). a Line n from which the second point is taken for extrapolation.
3
cpc cpc X 3 DEcorr: ðX Þ ðX 1Þ DEcorr: ðX 1Þ
X 3 ðX 1Þ
3
1 ¼ DEcorr:
ð3Þ
1 Table 1 shows different values for DEcorr: depending on the extrapolation points chosen. The most reliable extrapolation here (X ¼ 3 and 4, diffuse basis functions 1 on oxygen) predicts 18.7 kJ/mol for DEcorr: . When using cc-pVTZ (H, Al, Si) and aug-cc-pVQZ (O) basis sets 1 92% (17.2 kJ/mol, Table 1) of DEcorr: are recovered. Therefore we apply these basis sets to our 3T cluster model in the MP2 part of a second series of embedding calculations. To cope with the large number of basis functions (up to 1136) we use the ‘frozen core’ RIapproximation described before. The (RI)-MP2 calculations employ the Turbomole programme package [20–22].
3.3. MP2/DFT Using the QMPOT programme geometry optimisations are carried out at the combined MP2/DFT level. Two different basis set combinations (individually described above) are used: Ahlrichs’ basis sets for MP2 and 70 Ry plane wave cutoff for DFT (labelled ‘A’), as well as Dunning’s basis sets for MP2 and 100 Ry plane wave cutoff for DFT (labelled ‘B’). For comparison DFT-only calculations with CPMD are made. The convergence criterion during the geometry optimisation is a maximum cartesian component of the gradient of 1 104 Eh =a0 . Finally, numerical force constant calculations are made (cartesian displacements 0.01 a0 ) to characterise the stationary points found (DFT-only calculations with 70 and 100 Ry cutoff, and MP2/DFT with basis set A). Frequencies are calculated after projecting out the translational degrees of freedom.
4. Results and discussion For the proton jump from oxygen sites O1 to O2 both minimum structures with the proton located either on O1 or O2 and the corresponding transition state structures are located. This is done using the new combined MP2/DFT scheme and DFT-only (PBE) for comparison. Both the dry zeolite and the adsorption complex with water are considered. Table 2 lists selected atomic distances in the reactants, transition states and products. The combined MP2/DFT calculations lead to smaller O–H bond lengths in the reactants and products compared to the DFT-only (PBE) results. The distance between the water molecule and the acidic site is reduced as well. DFT-only (70 Ry), DFT-only (100 Ry) and MP2/ DFT (basis set A) yield imaginary frequencies for the transition state structures at 1247i, 1202i, and 1304i cm1 , respectively, for the dry zeolite (Fig. 1b), and at 251i, 283i, and 229i cm1 , respectively, for the hydrated system (Fig. 1e). Table 3 summarises the binding energies of the water molecule to the two investigated acidic sites in HSSZ-13. The DFT-only (100 Ry) results agree with the MP2/ DFT (basis set B, BSSE-corrected) results within 1 kJ/ mol. Since the GGA tends to overestimate hydrogen bonding but cannot account for dispersive interactions this is a result of error cancellation. For basis set A the BSSE in the MP2 part is about three times as big as for basis set B. Considering the results of our systematic basis set convergence study we can make a final estimate for the water binding energy by adding 1.5 kJ/mol (see Table 1: 18.7 kJ/mol from best extrapolation minus 17.2 kJ/mol for the actual basis set) to the BSSE-corrected MP2/DFT (basis set B) values arriving at 78 and 73 kJ/ mol for the proton sites O1 and O2 , respectively. To our knowledge these are the most accurate water binding
392
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
Table 2 Selected atomic distances (pm) for the proton jump from oxygen site O1 to O2 in HSSZ-13 calculated with our hybrid MP2/DFT scheme (basis sets A, B)
Dry r(O1 –H) r(O2 –H) Hydrated r(O1 –H) r(O1 Oa ) r(O2 Oa ) r(O2 –H) a
Reactant
Transition state
Product
97.4, 96.6 (97.8, 97.3)
123.4, 123.7 (124.1, 124.4) 125.1, 125.3 (127.2, 127.0)
97.6, 96.9 (98.1, 97.6)
148.6, 246.2, 241.7, 134.9,
250.5, 251.1 (251.3, 250.8) 104.9, 103.6 (107.6, 106.5)
103.7, 102.1 (105.7, 104.6) 251.7, 253.1 (254.5, 254.0)
(150.8, (249.7, (244.5, (134.3,
147.1) 246.9) 244.0) 135.9)
Values in parentheses are DFT-only results (70, 100 Ry). Oxygen atom of the water molecule.
Table 3 Electronic binding energies of water in the zeolite HSSZ-13 for different proton positions (kJ/mol) Method
H2 O H–O1
H2 O H–O2
DFT-only (70 Ry) DFT-only (100 Ry) MP2/DFT (basis set A) MP2/DFT (basis set B)
78 75 86 (71) 82 (76)
73 71 81 (63) 77 (72)
Values in parentheses are BSSE-corrected.
energies calculated for zeolites so far. After adding zeropoint vibrational and thermal corrections, the predicted heats of adsorption (298 K) are 73 and 69 kJ/mol, respectively. Experimental estimates are rather uncertain and range from 51 to 85 kJ/mol (Table 26 in [32]). Table 4 shows the results for the electronic barriers, DEz , and the reaction energies, DE, of the proton jump. Going from 70 to 100 Ry cutoff in the DFT-only calculations does not change energy differences by more than 3 kJ/mol. Comparison with the combined MP2/
Table 4 Electronic barriers and reaction energies, DEz and DE (kJ/mol), for the proton jump from oxygen site O1 to O2 in HSSZ-13 Method
Dry DE
z
Hydrated DE
DEz
DE
DFT-only (70 Ry) DFT-only (100 Ry)
65 (54) 68 (54)
10 (11) 11 (11)
20 (16) 22 (18)
14 (14) 15 (16)
MP2/DFT (basis set A) MP2(C) DFT(S)-DFT(C) DFT(S) MP2(C)-DFT(C)
71 (55) 53 18 66 5
11 (6) )8 19 10 1
26 (16) 23 3 20 6
16 (9) 10 6 13 3
MP2/DFT (basis set B)a MP2(C) DFT(S)-DFT(C) DFT(S) MP2(C)-DFT(C)
81 (65) 64 17 68 13
13 (8) 0 13 10 3
30 (20) 27 3 22 8
17 (10) 11 6 14 3
a
151.9 247.2 242.3 136.9
Values in parentheses are ZPE-corrected. ZPE corrections taken from calculations with basis set A.
DFT results confirms the expected underestimation of reaction barriers with the GGA density functional (PBE). With basis set A the barriers increase by 6 kJ/mol compared to DFT-only (70 Ry), and the qualitatively better basis set B even raises the 100 Ry DFT-only values by 8 to 13 kJ/mol. These are significant changes (e.g., increase of up to 36% for the hydrated zeolite) which become important when predictions for rates of such processes are made (see below). We also observe a small increase of 1 to 2 kJ/mol in the calculated reaction energies on application of the combined MP2/DFT approach which is assigned to cancellation effects since reactants and products are very similar. Within our MP2/DFT scheme the electronic energy of a system can be decomposed in two different ways (Table 4). The first one (see Eq. (1)) starts from the pure MP2 value for the 3T cluster model (C) calculated at the structure obtained by the combined MP2/DFT approach, MP2(C). The long-range correction term, defined as the DFT energy difference between the total system and the cluster model, DFT(S)-DFT(C), is expected to become smaller with growing size of the cluster model. For the proton jump in dry HSSZ-13 the long-range correction increases the MP2 cluster model result for the barrier significantly, and for the reaction energy it can even change the sign. Although being constrained to the geometry of the periodic solid the 3T model itself is too small to yield reliable results. For hydrated HSSZ-13 the effect of long-range corrections is smaller (but still important) for two reasons: the additional water molecule increases the size of the cluster model, and it attenuates the deformation of the entire zeolite lattice during the proton jump. The second energy decomposition in Table 4 starts from the pure DFT results for the complete periodic system (S) at the structure obtained with the combined method. They agree within 1 kJ/mol with the results of the DFT-only calculations showing that the structures do not change much upon switching to the combined MP2/DFT approach. The correction term, MP2(C)DFT(C), represents the local MP2 correction to the full
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
periodic DFT result (see Eq. (2)). This correction clearly increases the barriers as expected. Sierka and Sauer [13] used the B3LYP hybrid functional for full periodic single-point calculations. For the electronic energy barrier DEz and reaction DE (dry zeolite) they report a value of 72 and 12 kJ/mol, respectively. This is an improvement compared with the PBE functional used here (DFT-only calculations) and in good agreement to our MP2/DFT (basis set A) results, but the barrier is still underestimated compared to our MP2/DFT (basis set B) result of 81 kJ/mol. Sierka and Sauer [16] also introduced coupled cluster (CCSD(T)) corrections to their B3LYP results for a 1T cluster model and gained additional 13 kJ/mol for the proton jump barrier. This CCSD(T) correction includes an MP2 part of 8 kJ/mol (private communication to the authors) and we conclude that CCSD(T) corrections to our MP2/DFT calculations may increase the proton jump barrier by another 5 kJ/mol (dry zeolite). Zero-point vibrational energies (ZPE, see Table 4) reduce the energy barriers by 11–14 kJ/mol for the dry zeolite and by 4 kJ/mol for the water assisted proton jump (DFT-only calculations). At the MP2/DFT (basis set A) level these values are 16 and 10 kJ/mol, respectively. This means that the ZPE-corrected energy barriers are less different than the electronic energy barriers (DFT-only vs. MP2/DFT). To illustrate the importance of reliable predictions of energy barriers for the kinetics of proton jump reactions rate constants are calculated according to classical transition state theory, see Table 5. Using a cluster model for the zeolite HZSM-5, Ryder et al. [15] showed that the presence of water significantly lowers the energy barrier and therefore increases the rate for proton jumps compared to dry systems. For HSSZ-13 we obtain differences of about eight orders of magnitude between the rate constants for the dry and the hydrated zeolite at room
Table 5 Rate constants (per site) in s1 at different temperatures for the proton jump reaction from oxygen site O1 to O2 in HSSZ-13 for the dry zeolite and in presence of one water molecule Forward Dry
Reverse Hydrated
Dry
Hydrated
DFT-only (100 Ry) 298 K 2 103 400 K 5 105 500 K 2 107
7 108 3 109 9 109
2 105 2 107 3 108
8 1011 9 1011 9 1011
MP2/DFT (basis set B)a 298 K 4 101 400 K 4 104 500 K 2 106
3 109 3 1010 1 1011
9 102 3 105 1 107
5 109 1 1010 2 1010
a
Partition functions needed for the calculation of the rate constants taken from calculations with basis set A.
393
temperature. In a DFT molecular dynamics simulation of water in the zeolite gmelinite [14] spontaneous proton transfer is observed after 1.2 ps. The characteristic time for the reverse proton jump in hydrated HSSZ-13 calculated as reciprocal rate constant from the present data is also about 1.2 ps (DFT-only (100 Ry), Table 5). However, improving the potential energy surfaces by combined MP2/DFT embedding may change rate constants and, hence, characteristic times by up to two orders of magnitude.
5. Conclusion We implemented a hybrid QM/QM approach to improve the DFT potential energy surface of large chemical systems locally, i.e., at reaction sites, by MP2. In addition to the proton jump reaction studied here, the approach can also be used to calculate van der Waals contributions for the binding of small molecules to surfaces, e.g., [33], and for biomolecules. Extension beyond MP2 is possible when plane wave DFT codes are combined with codes for CCSD(T) gradients.
Acknowledgements This work has been supported by the Deutsche Forschungsgemeinschaft (SPP 1155). The computing facilities of the Norddeutscher Verbund f€ ur Hoch- und H€ ochstleistungsrechnen (HLRN) have been used. Special thanks go to Dr. B. Kallies (ZIB) for the OpenMP parallelisation of the rimp2 module of Turbomole.
References [1] C.J. Cramer, Essentials of Computational Chemistry, Wiley, New York, 2002. [2] Q. Zhang, R. Bell, T.H. Truong, J. Phys. Chem. 99 (1995) 592. [3] S. Sadhukhan, D. Mu~ noz, C. Adamo, G.E. Scuseria, Chem. Phys. Lett. 306 (1999) 83. [4] S. Chawla, G.A. Voth, J. Chem. Phys. 108 (1998) 4697. [5] P. Knowles, M. Sch€ utz, H.-J. Werner, in: J. Grotendorst (Ed.), Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 97. [6] W. Klopper, in: J. Grotendorst (Ed.), Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 181. [7] S. Humbel, S. Sieber, K. Morokuma, J. Chem. Phys. 105 (1996) 1959. [8] P. Sherwood, in: J. Grotendorst, Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 285. [9] J. Sauer, M. Sierka, J. Comput. Chem. 21 (2000) 1470. [10] M. Krossner, J. Sauer, J. Phys. Chem. 100 (1996) 6199.
394
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
[11] A.H. de Vries, P. Sherwood, S.J. Collins, A.M. Rigby, M. Rigutto, G.J. Kramer, J. Phys. Chem. B 103 (1999) 6133. [12] J.T. Fermann, C. Blanco, S. Auerbach, J. Chem. Phys. 112 (2000) 6779. [13] M. Sierka, J. Sauer, J. Chem. Phys. 112 (2000) 6983. Benco, T. Demuth, J. Hafner, F. Hutschka, Chem. Phys. Lett. [14] L. 324 (2000) 373. [15] J.A. Ryder, A.K. Chakraborty, A.T. Bell, J. Phys. Chem. B 104 (2000) 6998. [16] M. Sierka, J. Sauer, J. Phys. Chem. B 105 (2001) 1603. [17] A.M. Vos, F. De Proft, R.A. Schoonheydt, P. Geerlings, Chem. Commun. 12 (2001) 1108. [18] U. Eichler, C.M. K€ olmel, J. Sauer, J. Comput. Chem. 18 (1997) 463. [19] CPMD: J. Hutter, W. Andreoni, P. Ballone, M. Bernasconi, P. Bl€ ochl, A. Curioni, P. Focher, E. Fois, P. Giannozzi, S. Goedecker, D. Marx, M. Parrinello, U. R€ othlisberger, M. Tuckerman, MPI f€ ur Festk€ orperforschung Stuttgart and IBM Research Division Z€ urich, 1996.
[20] R. Ahlrichs, M. B€ar, M. H€aser, H. Horn, C. K€ olmel, Chem. Phys. Lett. 162 (1989) 165. [21] F. Haase, R. Ahlrichs, J. Comp. Chem. 14 (1993) 907. [22] F. Weigend, M. H€aser, Theor. Chem. Acc. 97 (1997) 331. [23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [24] S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1996) 1703. [25] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. [26] A. Sch€afer, H. Horn, R. Ahlrichs, J. Chem. Phys. 97 (1992) 2571. [27] T.H. Dunning, J. Chem. Phys. 90 (1989) 1007. [28] R.A. Kendall, T.H. Dunning, R.J. Harrison, J. Chem. Phys. 96 (1992) 6796. [29] D.E. Woon, T.H. Dunning, J. Chem. Phys. 98 (1993) 1358. [30] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [31] T. Helgaker, W. Klopper, H. Koch, J. Noga, J. Chem. Phys. 106 (1997) 9639. [32] J. Sauer, P. Ugliengo, E. Garrone, V.R. Saunders, Chem. Rev. 94 (1994) 2095. [33] P. Ugliengo, A. Damin, Chem. Phys. Lett. 366 (2002) 683.
A hybrid MP2/planewave-DFT scheme for large chemical systems: proton jumps in zeolites Christian Tuma, Joachim Sauer
*
Humboldt-Universit€at zu Berlin, Institut f€ur Chemie, Unter den Linden 6, Berlin D-10099, Germany Received 20 January 2004; in final form 17 February 2004 Published online: 10 March 2004
Abstract We present an embedding scheme to introduce local corrections at post Hartree–Fock level to density functional theory (DFT) calculations. As a first application we study proton jump reactions in the zeolite HSSZ-13 and show that energy barriers and rate constants are significantly changed by second-order Møller–Plesset perturbation theory (MP2) corrections to plane wave based DFT calculations. Electronic energy barriers increase from 68 to 81 kJ/mol (dry zeolite), and from 22 to 30 kJ/mol (hydrated zeolite). The predicted heats of adsorption of one water molecule onto the Brønsted acidic sites O1 and O2 are 73 and 69 kJ/mol, respectively. Ó 2004 Elsevier B.V. All rights reserved.
1. Introduction In electronic structure studies of materials and molecular systems the use of density functional theory (DFT) has met with great success. The main advantage of DFT is its computational efficiency and scaling behaviour which, especially in combination with a plane wave basis set, allows the application to extended molecular systems or to solids. For chemical problems (in particular involving transition states) hybrid density functionals are often the best choice among available functionals, see e.g., [1–3]. Unfortunately, the calculation of the exact exchange part is computationally too expensive with plane waves (e.g., [4]), and one is limited to various gradient-corrected functionals. Although the plane wave basis set quality can be improved systematically by increasing the kinetic energy cutoff, DFT itself does not offer methodological hierarchies of increased accuracy as post Hartree–Fock (HF) methods [5,6] do. The motivation of our work is the desire to combine the computational efficiency of plane wave based DFT *
Corresponding author. Fax: +493020937136. E-mail address: [email protected] (J. Sauer). URL: http://www.chemie.hu-berlin.de/ag_sauer.
0009-2614/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2004.02.056
with the accuracy and systematic convergence of post HF quantum mechanical (QM) methods. We present a mechanical embedding technique for chemical systems with an important fragment, e.g., a reaction site. We use DFT in a plane wave approach as a ‘low’-level method and add local corrections at post HF level using secondorder Møller–Plesset perturbation theory (MP2). This hybrid QM/QM approach extends the pool of existing high-level/low-level combinations [7–9]. As a first application we study a proton jump reaction in the zeolite HSSZ-13 (the acidic form of the zeolite chabasite), either loaded with a single water molecule per Brønsted site or being completely dry. This type of reaction has already been studied using different approaches [10–17]. Without inclusion of exact exchange, DFT calculations underestimate reaction barriers for proton jump reactions. For example, using a 3T cluster model, Auerbach et al. [12] report 72.8, 90.2 and 99.1 kJ/ mol for the proton jump barrier in the zeolite H–Y when BLYP, B3LYP, and MP2, respectively, are used. Ryder et al. [15] investigated a proton jump in the zeolite HZSM-5 employing a 5T cluster model and found energy barriers of 85.8, 108.4 and 135.1 kJ/mol (BPW91, B3LYP, BH&HLYP, respectively). For water assisted proton jumps the cluster model study by Krossner and Sauer [10] yielded 8.1 and 16.4 kJ/mol (BP86, MP2,
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
respectively) for the energy barrier. These examples illustrate the necessity to introduce at least local improvements over GGA type density functionals in studies of proton jump barriers.
2. Method We follow the mechanical embedding scheme described in [13,18]. The total system S (e.g., unit cell or large molecule) is divided into an inner part containing the active site and an outer part. When bonds are cut on definition of the inner part, link atoms L are added to saturate the dangling bonds of the inner part. The inner part and the link atoms together form the cluster model C. The total system S and the cluster model C are calculated at ‘low’-level, in this study DFT with a plane wave basis set and atomic pseudopotentials. The cluster model C is also calculated at ‘high’-level, here MP2 with standard gaussian basis sets. The new energy hypersurface is defined by a subtraction scheme as EðS; LÞ ¼ EðCÞhigh þ ½EðSÞlow EðCÞlow ;
ð1Þ
or equivalently EðS; LÞ ¼ EðSÞlow þ ½EðCÞhigh EðCÞlow :
ð2Þ
Our embedding scheme imposes constraints on the positions of the link atoms L. This implies additional contributions to forces on atoms of broken bonds which ensures translational and rotational invariance of the derivatives of EðS; LÞ. The present embedding scheme differs from the IMOMO method of Morokuma et al. [7] and other hybrid approaches [8,9] by details of the implementation of constraints. The term in squared brackets in Eq. (1) corresponds to the ‘low’-level long-range correction to the ‘high’level cluster result. It accounts for both long-range electronic effects and structural constraints imposed by the zeolite framework on the reaction site. Small cluster models without embedding miss such corrections and therefore neglect specific characteristics of the framework and its potential impact on the problem studied. Another view is given in Eq. (2). The term in squared brackets can be considered as the local ‘high’-level correction (here MP2) to the full periodic ‘low’-level (in this study DFT) result. Our mechanical embedding scheme is implemented in the QMPOT code [13]. For the present QM/QM coupling QMPOT has been extended with an additional interface to the CPMD code [19]. Along with its interface to the Turbomole programme package [20–22] we use QMPOT for combined MP2/DFT calculations.
389
3. Computational details 3.1. DFT The GGA by Perdew, Burke, and Ernzerhof (PBE) [23] is employed in connection with a computational efficient Pade form for exchange and correlation [24]. For the silicon, aluminium, oxygen, and hydrogen atoms of the zeolite system norm-conserving pseudopotentials are used (Troullier and Martins [25]) with core radii of 1.4, 1.6, 1.15, and 0.7 a0 , respectively. The kinetic energy cutoff of the plane wave basis set is 70 Ry. This value is commonly used in DFT studies of aqueous systems with this type of pseudopotential. A series of extensive test calculations showed that the presence of oxygen requires a higher cutoff of about 100 Ry for well converged DFT results. For this reason we also perform calculations at 100 Ry cutoff. All DFT plane wave calculations use the CPMD code [19]. The unit cell parameters of HSSZ-13 (periodic boundary conditions, stoichiometry AlHSi11 O24 per unit cell) have been obtained from ion-pair shell-model potential calculations with the hydrogen atom at its energetically favoured O1 position: a ¼ 945:3 pm, b ¼ 938:7 pm, c ¼ 944:1 pm, a ¼ 94:0°, b ¼ 95:4°, and c ¼ 94:8°, see [13] for details. In our calculations these cell parameters are kept constant. The cluster model (described in Section 3.2) is treated as an isolated system, i.e., without periodic boundary conditions. To enclose the electrons of this molecular system virtually completely, a computational box bigger than the unit cell of HSSZ-13 is required. This means that the basis sets for the periodic system and the cluster model include a different number of plane waves. However, test calculations on isolated molecules show that the quality of the plane wave basis is determined by its kinetic energy cutoff and not by the number of plane waves. For the cluster model convergence in the total energy is achieved for computational boxes bigger than Under this condition total energy fluctuations are 12 A. about 1.5 and 0.2 mEh at 70 and 100 Ry cutoff, respectively, which can be neglected in the context of our embedding scheme. We use a computational box of 14 A for the water molefor the cluster model, and of 10 A cule. 3.2. MP2 The cluster model for the ‘high’-level calculations consists of three tetrahedral Si[O1=2 ]4 units (3T). It is highlighted using a ball-and-stick presentation in Fig. 1. Eight hydrogen link atoms HL have been added to saturate dangling bonds of the inner part. The HL –OT bond lengths have been fixed at 95.3 pm (T ¼ Al) and 95.4 pm (T ¼ Si), mean values from bond length optimisations on a 2T cluster model. Ahlrichs’ fully
390
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
Fig. 1. (a–f) Proton jump reaction between oxygen sites O1 and O2 in the zeolite HSSZ-13 without (a–c) and in the presence of a water molecule (d–f). Structures for initial and final states (a,c,d,f) and transition state structures (b,e) are shown, the ball-and-stick parts represent the 3T cluster model for the MP2/DFT mechanical embedding calculations (see text for details).
optimised triple-f basis set is used for oxygen and double-f basis sets for the remaining elements [26]. Polarisation functions with exponents 1.2 (O), 0.8 (H), 0.35 (Si), and 0.3 (Al) have been added. In the MP2 correlation scheme all electrons are included. This basis set can be regarded as a small one for MP2 calculations. Higher angular momentum functions are necessary to achieve reasonable MP2 results, and diffuse functions are needed to reduce the basis set superposition error (BSSE). Therefore, we made a systematic MP2 study applying Dunning’s correlation consistent (aug-)ccpVXZ basis sets [27–29] to a 2T + H2 O cluster model. Single point MP2 calculations for the water binding energy and corresponding counterpoise corrections to
account for the BSSE [30] are made using a RI (resolution of the identity) implementation [22]. Electrons in the 1s orbitals of Al, Si, and O are not explicitly correlated (‘frozen core’). The results (Table 1) show that the basis set for the oxygen atoms is decisive for the BSSE, i.e., at least triple-f quality and one set of diffuse functions are necessary. Table 1 also includes the MP2 correlation energy parts of the binding energies corrected cpc for the BSSE, DEcorr: . The results for the large basis sets show that most of the BSSE (up to 90%) arises from the cpc correlation energy. Taking the DEcorr: values for different cpc basis sets, DEcorr: ðX Þ, X ¼ 2; 3; 4 (‘D;T;Q’), we carried out 2-point extrapolations to the complete basis set 1 limit, DEcorr: [31]:
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
391
Table 1 Single point RI-MP2 calculations for water adsorption on a 2T cluster model using different basis sets n
H, Al, Si
O
DE
DEcpc
DEcorr:
cpc DEcorr:
f
ðnÞa
1 DEcorr:
1 2 3 4 5 6 7 8 9 10 11 12
cc-pVDZ
cc-pVDZ cc-pVTZ cc-pVQZ aug-cc-pVDZ aug-cc-pVTZ aug-cc-pVQZ cc-pVTZ cc-pVQZ aug-cc-pVTZ aug-cc-pVQZ cc-pVQZ aug-cc-pVQZ
88.0 72.6 67.6 64.4 66.3 67.1 74.5 68.7 65.3 65.5 69.7 65.9
53.1 57.5 59.5 54.6 58.9 61.3 58.5 60.6 59.9 61.6 61.3 61.8
21.6 20.5 20.1 17.7 20.0 21.3 22.9 21.1 20.1 20.7 22.5 21.2
5.9 12.6 15.5 11.4 15.3 17.1 13.7 16.0 15.8 17.2 16.7 17.5
0.45 0.52 0.57 0.64 0.64 0.72 0.58 0.63 0.80 0.90 0.69 0.90
(1) (2) (4) (5) (7) (9)
17.0 18.5 17.6 18.6 18.8 18.7
cc-pVTZ
cc-pVQZ
Water binding energies and corresponding BSSE-corrected values, DE and DEcpc , as well as their RI-MP2 correlation energy parts, DEcorr: and cpc 1 , are given. The fraction of the RI-MP2 correlation based BSSE on the full BSSE, f , is shown. DEcorr: is the 2-point extrapolated value for the DEcorr: complete basis set limit of the RI-MP2 correlation part of the water binding energy (all energies in kJ/mol). a Line n from which the second point is taken for extrapolation.
3
cpc cpc X 3 DEcorr: ðX Þ ðX 1Þ DEcorr: ðX 1Þ
X 3 ðX 1Þ
3
1 ¼ DEcorr:
ð3Þ
1 Table 1 shows different values for DEcorr: depending on the extrapolation points chosen. The most reliable extrapolation here (X ¼ 3 and 4, diffuse basis functions 1 on oxygen) predicts 18.7 kJ/mol for DEcorr: . When using cc-pVTZ (H, Al, Si) and aug-cc-pVQZ (O) basis sets 1 92% (17.2 kJ/mol, Table 1) of DEcorr: are recovered. Therefore we apply these basis sets to our 3T cluster model in the MP2 part of a second series of embedding calculations. To cope with the large number of basis functions (up to 1136) we use the ‘frozen core’ RIapproximation described before. The (RI)-MP2 calculations employ the Turbomole programme package [20–22].
3.3. MP2/DFT Using the QMPOT programme geometry optimisations are carried out at the combined MP2/DFT level. Two different basis set combinations (individually described above) are used: Ahlrichs’ basis sets for MP2 and 70 Ry plane wave cutoff for DFT (labelled ‘A’), as well as Dunning’s basis sets for MP2 and 100 Ry plane wave cutoff for DFT (labelled ‘B’). For comparison DFT-only calculations with CPMD are made. The convergence criterion during the geometry optimisation is a maximum cartesian component of the gradient of 1 104 Eh =a0 . Finally, numerical force constant calculations are made (cartesian displacements 0.01 a0 ) to characterise the stationary points found (DFT-only calculations with 70 and 100 Ry cutoff, and MP2/DFT with basis set A). Frequencies are calculated after projecting out the translational degrees of freedom.
4. Results and discussion For the proton jump from oxygen sites O1 to O2 both minimum structures with the proton located either on O1 or O2 and the corresponding transition state structures are located. This is done using the new combined MP2/DFT scheme and DFT-only (PBE) for comparison. Both the dry zeolite and the adsorption complex with water are considered. Table 2 lists selected atomic distances in the reactants, transition states and products. The combined MP2/DFT calculations lead to smaller O–H bond lengths in the reactants and products compared to the DFT-only (PBE) results. The distance between the water molecule and the acidic site is reduced as well. DFT-only (70 Ry), DFT-only (100 Ry) and MP2/ DFT (basis set A) yield imaginary frequencies for the transition state structures at 1247i, 1202i, and 1304i cm1 , respectively, for the dry zeolite (Fig. 1b), and at 251i, 283i, and 229i cm1 , respectively, for the hydrated system (Fig. 1e). Table 3 summarises the binding energies of the water molecule to the two investigated acidic sites in HSSZ-13. The DFT-only (100 Ry) results agree with the MP2/ DFT (basis set B, BSSE-corrected) results within 1 kJ/ mol. Since the GGA tends to overestimate hydrogen bonding but cannot account for dispersive interactions this is a result of error cancellation. For basis set A the BSSE in the MP2 part is about three times as big as for basis set B. Considering the results of our systematic basis set convergence study we can make a final estimate for the water binding energy by adding 1.5 kJ/mol (see Table 1: 18.7 kJ/mol from best extrapolation minus 17.2 kJ/mol for the actual basis set) to the BSSE-corrected MP2/DFT (basis set B) values arriving at 78 and 73 kJ/ mol for the proton sites O1 and O2 , respectively. To our knowledge these are the most accurate water binding
392
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
Table 2 Selected atomic distances (pm) for the proton jump from oxygen site O1 to O2 in HSSZ-13 calculated with our hybrid MP2/DFT scheme (basis sets A, B)
Dry r(O1 –H) r(O2 –H) Hydrated r(O1 –H) r(O1 Oa ) r(O2 Oa ) r(O2 –H) a
Reactant
Transition state
Product
97.4, 96.6 (97.8, 97.3)
123.4, 123.7 (124.1, 124.4) 125.1, 125.3 (127.2, 127.0)
97.6, 96.9 (98.1, 97.6)
148.6, 246.2, 241.7, 134.9,
250.5, 251.1 (251.3, 250.8) 104.9, 103.6 (107.6, 106.5)
103.7, 102.1 (105.7, 104.6) 251.7, 253.1 (254.5, 254.0)
(150.8, (249.7, (244.5, (134.3,
147.1) 246.9) 244.0) 135.9)
Values in parentheses are DFT-only results (70, 100 Ry). Oxygen atom of the water molecule.
Table 3 Electronic binding energies of water in the zeolite HSSZ-13 for different proton positions (kJ/mol) Method
H2 O H–O1
H2 O H–O2
DFT-only (70 Ry) DFT-only (100 Ry) MP2/DFT (basis set A) MP2/DFT (basis set B)
78 75 86 (71) 82 (76)
73 71 81 (63) 77 (72)
Values in parentheses are BSSE-corrected.
energies calculated for zeolites so far. After adding zeropoint vibrational and thermal corrections, the predicted heats of adsorption (298 K) are 73 and 69 kJ/mol, respectively. Experimental estimates are rather uncertain and range from 51 to 85 kJ/mol (Table 26 in [32]). Table 4 shows the results for the electronic barriers, DEz , and the reaction energies, DE, of the proton jump. Going from 70 to 100 Ry cutoff in the DFT-only calculations does not change energy differences by more than 3 kJ/mol. Comparison with the combined MP2/
Table 4 Electronic barriers and reaction energies, DEz and DE (kJ/mol), for the proton jump from oxygen site O1 to O2 in HSSZ-13 Method
Dry DE
z
Hydrated DE
DEz
DE
DFT-only (70 Ry) DFT-only (100 Ry)
65 (54) 68 (54)
10 (11) 11 (11)
20 (16) 22 (18)
14 (14) 15 (16)
MP2/DFT (basis set A) MP2(C) DFT(S)-DFT(C) DFT(S) MP2(C)-DFT(C)
71 (55) 53 18 66 5
11 (6) )8 19 10 1
26 (16) 23 3 20 6
16 (9) 10 6 13 3
MP2/DFT (basis set B)a MP2(C) DFT(S)-DFT(C) DFT(S) MP2(C)-DFT(C)
81 (65) 64 17 68 13
13 (8) 0 13 10 3
30 (20) 27 3 22 8
17 (10) 11 6 14 3
a
151.9 247.2 242.3 136.9
Values in parentheses are ZPE-corrected. ZPE corrections taken from calculations with basis set A.
DFT results confirms the expected underestimation of reaction barriers with the GGA density functional (PBE). With basis set A the barriers increase by 6 kJ/mol compared to DFT-only (70 Ry), and the qualitatively better basis set B even raises the 100 Ry DFT-only values by 8 to 13 kJ/mol. These are significant changes (e.g., increase of up to 36% for the hydrated zeolite) which become important when predictions for rates of such processes are made (see below). We also observe a small increase of 1 to 2 kJ/mol in the calculated reaction energies on application of the combined MP2/DFT approach which is assigned to cancellation effects since reactants and products are very similar. Within our MP2/DFT scheme the electronic energy of a system can be decomposed in two different ways (Table 4). The first one (see Eq. (1)) starts from the pure MP2 value for the 3T cluster model (C) calculated at the structure obtained by the combined MP2/DFT approach, MP2(C). The long-range correction term, defined as the DFT energy difference between the total system and the cluster model, DFT(S)-DFT(C), is expected to become smaller with growing size of the cluster model. For the proton jump in dry HSSZ-13 the long-range correction increases the MP2 cluster model result for the barrier significantly, and for the reaction energy it can even change the sign. Although being constrained to the geometry of the periodic solid the 3T model itself is too small to yield reliable results. For hydrated HSSZ-13 the effect of long-range corrections is smaller (but still important) for two reasons: the additional water molecule increases the size of the cluster model, and it attenuates the deformation of the entire zeolite lattice during the proton jump. The second energy decomposition in Table 4 starts from the pure DFT results for the complete periodic system (S) at the structure obtained with the combined method. They agree within 1 kJ/mol with the results of the DFT-only calculations showing that the structures do not change much upon switching to the combined MP2/DFT approach. The correction term, MP2(C)DFT(C), represents the local MP2 correction to the full
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
periodic DFT result (see Eq. (2)). This correction clearly increases the barriers as expected. Sierka and Sauer [13] used the B3LYP hybrid functional for full periodic single-point calculations. For the electronic energy barrier DEz and reaction DE (dry zeolite) they report a value of 72 and 12 kJ/mol, respectively. This is an improvement compared with the PBE functional used here (DFT-only calculations) and in good agreement to our MP2/DFT (basis set A) results, but the barrier is still underestimated compared to our MP2/DFT (basis set B) result of 81 kJ/mol. Sierka and Sauer [16] also introduced coupled cluster (CCSD(T)) corrections to their B3LYP results for a 1T cluster model and gained additional 13 kJ/mol for the proton jump barrier. This CCSD(T) correction includes an MP2 part of 8 kJ/mol (private communication to the authors) and we conclude that CCSD(T) corrections to our MP2/DFT calculations may increase the proton jump barrier by another 5 kJ/mol (dry zeolite). Zero-point vibrational energies (ZPE, see Table 4) reduce the energy barriers by 11–14 kJ/mol for the dry zeolite and by 4 kJ/mol for the water assisted proton jump (DFT-only calculations). At the MP2/DFT (basis set A) level these values are 16 and 10 kJ/mol, respectively. This means that the ZPE-corrected energy barriers are less different than the electronic energy barriers (DFT-only vs. MP2/DFT). To illustrate the importance of reliable predictions of energy barriers for the kinetics of proton jump reactions rate constants are calculated according to classical transition state theory, see Table 5. Using a cluster model for the zeolite HZSM-5, Ryder et al. [15] showed that the presence of water significantly lowers the energy barrier and therefore increases the rate for proton jumps compared to dry systems. For HSSZ-13 we obtain differences of about eight orders of magnitude between the rate constants for the dry and the hydrated zeolite at room
Table 5 Rate constants (per site) in s1 at different temperatures for the proton jump reaction from oxygen site O1 to O2 in HSSZ-13 for the dry zeolite and in presence of one water molecule Forward Dry
Reverse Hydrated
Dry
Hydrated
DFT-only (100 Ry) 298 K 2 103 400 K 5 105 500 K 2 107
7 108 3 109 9 109
2 105 2 107 3 108
8 1011 9 1011 9 1011
MP2/DFT (basis set B)a 298 K 4 101 400 K 4 104 500 K 2 106
3 109 3 1010 1 1011
9 102 3 105 1 107
5 109 1 1010 2 1010
a
Partition functions needed for the calculation of the rate constants taken from calculations with basis set A.
393
temperature. In a DFT molecular dynamics simulation of water in the zeolite gmelinite [14] spontaneous proton transfer is observed after 1.2 ps. The characteristic time for the reverse proton jump in hydrated HSSZ-13 calculated as reciprocal rate constant from the present data is also about 1.2 ps (DFT-only (100 Ry), Table 5). However, improving the potential energy surfaces by combined MP2/DFT embedding may change rate constants and, hence, characteristic times by up to two orders of magnitude.
5. Conclusion We implemented a hybrid QM/QM approach to improve the DFT potential energy surface of large chemical systems locally, i.e., at reaction sites, by MP2. In addition to the proton jump reaction studied here, the approach can also be used to calculate van der Waals contributions for the binding of small molecules to surfaces, e.g., [33], and for biomolecules. Extension beyond MP2 is possible when plane wave DFT codes are combined with codes for CCSD(T) gradients.
Acknowledgements This work has been supported by the Deutsche Forschungsgemeinschaft (SPP 1155). The computing facilities of the Norddeutscher Verbund f€ ur Hoch- und H€ ochstleistungsrechnen (HLRN) have been used. Special thanks go to Dr. B. Kallies (ZIB) for the OpenMP parallelisation of the rimp2 module of Turbomole.
References [1] C.J. Cramer, Essentials of Computational Chemistry, Wiley, New York, 2002. [2] Q. Zhang, R. Bell, T.H. Truong, J. Phys. Chem. 99 (1995) 592. [3] S. Sadhukhan, D. Mu~ noz, C. Adamo, G.E. Scuseria, Chem. Phys. Lett. 306 (1999) 83. [4] S. Chawla, G.A. Voth, J. Chem. Phys. 108 (1998) 4697. [5] P. Knowles, M. Sch€ utz, H.-J. Werner, in: J. Grotendorst (Ed.), Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 97. [6] W. Klopper, in: J. Grotendorst (Ed.), Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 181. [7] S. Humbel, S. Sieber, K. Morokuma, J. Chem. Phys. 105 (1996) 1959. [8] P. Sherwood, in: J. Grotendorst, Modern Methods and Algorithms of Quantum Chemistry, Proceedings, NIC Series, vol. 3, John von Neumann Institute for Computing, Forschungszentrum J€ ulich, 2000, p. 285. [9] J. Sauer, M. Sierka, J. Comput. Chem. 21 (2000) 1470. [10] M. Krossner, J. Sauer, J. Phys. Chem. 100 (1996) 6199.
394
C. Tuma, J. Sauer / Chemical Physics Letters 387 (2004) 388–394
[11] A.H. de Vries, P. Sherwood, S.J. Collins, A.M. Rigby, M. Rigutto, G.J. Kramer, J. Phys. Chem. B 103 (1999) 6133. [12] J.T. Fermann, C. Blanco, S. Auerbach, J. Chem. Phys. 112 (2000) 6779. [13] M. Sierka, J. Sauer, J. Chem. Phys. 112 (2000) 6983. Benco, T. Demuth, J. Hafner, F. Hutschka, Chem. Phys. Lett. [14] L. 324 (2000) 373. [15] J.A. Ryder, A.K. Chakraborty, A.T. Bell, J. Phys. Chem. B 104 (2000) 6998. [16] M. Sierka, J. Sauer, J. Phys. Chem. B 105 (2001) 1603. [17] A.M. Vos, F. De Proft, R.A. Schoonheydt, P. Geerlings, Chem. Commun. 12 (2001) 1108. [18] U. Eichler, C.M. K€ olmel, J. Sauer, J. Comput. Chem. 18 (1997) 463. [19] CPMD: J. Hutter, W. Andreoni, P. Ballone, M. Bernasconi, P. Bl€ ochl, A. Curioni, P. Focher, E. Fois, P. Giannozzi, S. Goedecker, D. Marx, M. Parrinello, U. R€ othlisberger, M. Tuckerman, MPI f€ ur Festk€ orperforschung Stuttgart and IBM Research Division Z€ urich, 1996.
[20] R. Ahlrichs, M. B€ar, M. H€aser, H. Horn, C. K€ olmel, Chem. Phys. Lett. 162 (1989) 165. [21] F. Haase, R. Ahlrichs, J. Comp. Chem. 14 (1993) 907. [22] F. Weigend, M. H€aser, Theor. Chem. Acc. 97 (1997) 331. [23] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [24] S. Goedecker, M. Teter, J. Hutter, Phys. Rev. B 54 (1996) 1703. [25] N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. [26] A. Sch€afer, H. Horn, R. Ahlrichs, J. Chem. Phys. 97 (1992) 2571. [27] T.H. Dunning, J. Chem. Phys. 90 (1989) 1007. [28] R.A. Kendall, T.H. Dunning, R.J. Harrison, J. Chem. Phys. 96 (1992) 6796. [29] D.E. Woon, T.H. Dunning, J. Chem. Phys. 98 (1993) 1358. [30] S.F. Boys, F. Bernardi, Mol. Phys. 19 (1970) 553. [31] T. Helgaker, W. Klopper, H. Koch, J. Noga, J. Chem. Phys. 106 (1997) 9639. [32] J. Sauer, P. Ugliengo, E. Garrone, V.R. Saunders, Chem. Rev. 94 (1994) 2095. [33] P. Ugliengo, A. Damin, Chem. Phys. Lett. 366 (2002) 683.