Ab initio benchmark study for the oxidative addition of CH 4 to Pd: Importance of basis-set flexibility and polarization

To obtain a state-of-the-art benchmark potential energy surface ~PES! for the archetypal oxidative addition of the methane C-H bond to the palladium atom, we have explored this PES using a hierarchical series of ab initio methods„Hartree-Fock, second-order Møller-Plesset perturbation theory, fourth-order Møller-Plesset perturbation theory with single, double and quadruple excitations, coupled cluster theory with single and double excitations ~CCSD!, and with triple excitations treated perturbatively @CCSD~T!#... and hybrid density functional theory using the B3LYP functional, in combination with a hierarchical series of ten Gaussian-type basis sets, up to g polarization. Relativistic effects are taken into account either through a relativistic effective core potential for palladium or through a full four-component all-electron approach. Counterpoise corrected relative energies of stationary points are converged to within 0.1–0.2 kcal/mol as a function of the basis-set size. Our best estimate of kinetic and thermodynamic parameters is 28.1 ~28.3! kcal/mol for the formation of the reactant complex, 5.8 ~3.1! kcal/mol for the activation energy relative to the separate reactants, and 0.8 ~21.2! kcal/mol for the reaction energy ~zero-point vibrational energy-corrected values in parentheses !. This agrees well with available experimental data. Our work highlights the importance of sufficient higher angular momentum polarization functions, f and g, for correctly describing metal– d-electron correlation and, thus, for obtaining reliable relative energies. We show that standard basis sets, such as LANL2DZ 11 f for palladium, are not sufficiently polarized for this purpose and lead to erroneous CCSD ~T! results. B3LYP is associated with smaller basis set superposition errors and shows faster convergence with basis-set size but yields relative energies ~in particular, a reaction barrier ! that are ca. 3.5 kcal/mol higher than the corresponding CCSD ~T! values. © 2004 American Institute of Physics. @DOI: 10.1063/1.1792151 #


I. INTRODUCTION
The activation of the C-H bond in alkanes is a challenging and important goal of catalysis. It is often the first step in the catalytic conversion of the abundant and nonreactive alkanes into more useful products. 1 Chemical intuition tells us that it is difficult to let the C-H bond be activated by transition metal atoms, or to let it even be coordinated to them. Alkanes are poor electron donors and acceptors. The alkane C-H bond is strong and nonpolar. Because the highest occupied molecular orbital is low lying, it is unsuitable for electron donation, whereas the high-lying * lowest unoccupied molecular orbital is unsuitable for accepting electron density. 2 In the group of the transition metal elements palladium is one of the most important catalysts, mostly in conjunction with ligands. 3 The insertion of the palladium atom into the C-H bond in alkanes has therefore received considerable attention, both experimentally 4 -8 and theoretically. 5,7,[9][10][11][12][13][14][15][16] In this work the insertion of the Pd-d 10 atom into the C-H bond of methane as an important example of this type of reactions has been surveyed, see Chart 1. Experimental investigations on the kinetics of the reaction of palladium with alkanes have been carried out by Weisshaar and co-workers 6,7 using laser-induced fluorescence ͑LIF͒ techniques and, most recently, by Campbell. 8 These studies show that Pd forms collisionally stabilized a͒ FAX: ϩ31-20-44 47629. Electronic mail: FM.Bickelhaupt@few.vu.nl complexes with alkanes, in particular also methane, 8 and that the reaction rate is extremely small to negligible. The exponential decay of the Pd signal vs alkane pressure suggests a complexation energy of at least 8 kcal/mol for Pd-alkane complexes. 6 This provides us with an experimental boundary condition for the stability of the reactant complex of Pd ϩmethane.
The purpose of the present study is twofold. In the first place, we wish to obtain a reliable benchmark for the potential energy surface ͑PES͒ for the archetypal organometallic reaction of methane oxidative addition to Pd͑0͒. This is done by exploring this PES with a hierarchical series of ab initio and hybrid density functional methods "Hartree-Fock ͑HF͒, the B3LYP functional, second-order Møller-Plesset perturbation theory ͑MP2͒, fourth-order Møller-Plesset perturbation theory with single, double, and quadruple excitations ͑MP4SDQ͒, coupled cluster theory with single and double excitations ͑CCSD͒, and with triple excitations treated perturbatively ͓CCSD͑T͔͒… in combination with a hierarchical series of ten Gaussian-type basis sets of increasing flexibility and polarization ͑up to g functions͒. The basis set superposition error ͑BSSE͒ is accounted for by counterpoise correction ͑CPC͒. 17 Relativistic effects were shown before to be important for the present model reaction. 15 Here, they are treated either with a relativistic effective core potential for palladium or with a full four-component all-electron approach ͑see Sec. II͒. The existing computational benchmark for oxidative addition of methane to palladium was obtained by Blomberg and co-workers 7 with the parametrized configuration interaction method PCI-80, 18 in which the effect of correlation is estimated by an extrapolation procedure. The PCI-80 study arrives at a Pdϩmethane complexation energy of 5.1 kcal/mol, an activation energy of 3.6 kcal/mol with respect to separate reactants, and a reaction energy of 2.3 kcal/mol ͑see Table I͒.
These values and, in particular, the activation energies appear to be highly sensitive to the level of theory used. A spectrum of values has been computed, for example, for the activation energy of this reaction that ranges from ϩ30.5 to Ϫ3.8 kcal/mol ͑see Table I͒. In view of this situation, it is appropriate to explore to which extent the PCI-80 values are converged with respect to both the order of correlation incorporated into the theoretical method and the degree of flexibility and polarization of the basis set. The present study also serves to clarify this issue. Note also that, in addition to the extrapolation procedure associated with PCI-80, the computation of these benchmark values involves a further approximation, namely the final scaled MCPF energies of the PCI-80 study were not computed at the MCPF but, instead, the HF optimum geometry. 7,14 In the present study, we use a consistent set of equilibrium and transition-state structures that have been fully optimized at the BLYP/TZ2P ͑see Sec. II A͒ level of the generalized gradient approximation ͑GGA͒ of density functional theory ͑DFT͒.
A second purpose of our work is to find out how well standard basis sets designed for use with high-level correlated ab initio methods such as CCSD͑T͒ are suited for correctly describing correlation phenomena associated with organometallic reactions involving bond breaking processes. The activation of the methane C-H bond by palladium serves as a test case. At this point, however, we anticipate that our findings have implications beyond the scope of this model reaction.

A. Geometries
Geometries of the stationary points were optimized with the ADF program 19,20 using the GGA to DFT ͑Ref. 21͒ at BLYP, 22,23 in combination with a large uncontracted set of Slater-type orbitals ͑STOs͒ containing diffuse functions: TZ2P. The TZ2P basis is of triple-quality and has been augmented with two sets of polarization functions: 2p and 3d on hydrogen, 3d and 4 f on carbon, and 5p and 4 f on palladium. The core shells of carbon (1s) and palladium (1s2s2p3s3p3d) were treated by the frozen-core approximation. 19 An auxiliary set of s, p, d, f, and g STOs was used to fit the molecular density and to represent the Coulomb and exchange potentials accurately in each selfconsistent-field ͑SCF͒ cycle. 19 Relativistic effects were accounted for using the zeroth-order regular approximation ͑ZORA͒. 24 Through a vibrational analysis, stationary points were confirmed to be equilibrium structures ͑no imaginary frequencies͒ or a transition state ͑one imaginary frequency͒.

B. Ab initio methods
Energies of the stationary points were computed with the program packages GAUSSIAN 25 and DIRAC 26 using the following hierarchy of quantum chemical methods: HF, MP2, 27 MP4SDQ, 28 CCSD 29,30 and CCSD͑T͒. 31 Finally, DFT calculations have been done with the B3LYP functional. 23,32 In calculations with the GAUSSIAN program, relativistic effects are described using a relativistic effective core potential ͑ECP͒ for palladium ͑vide infra͒. On the other hand, in calculations with the DIRAC program, relativistic effects are accounted for using a full all-electron four-component Dirac-Coulomb approach with a spin-free Hamiltonian ͑SFDC͒. 33 The two-electron integrals over the small components have been neglected and corrected with a simple Coulombic correction, which has been shown reliable. 34

C. Basis sets
For carbon and hydrogen, we used Dunning's correlation consistent augmented double-͑cc-aug-pVDZ͒ and triple-͑cc-aug-pVTZ͒ basis sets 35 in both GAUSSIAN and DIRAC calculations. For palladium, two different types of basis sets are used for the two programs, which leads to two series of basis sets for our model system: A1-A4 in the GAUSSIAN calculations and B1-B6 in the DIRAC calculations ͑see Table II͒. The series A1-A4 in the GAUSSIAN calculations are based on the Gaussian-type LANL2DZ basis set of Hay and Wadt for palladium. 36 This basis set involves a relativistic ECP that accounts for mass-velocity and Darwin terms. Basis set A1 corresponds to cc-aug-pVDZ for C and H and the standard LANL2DZ basis set for Pd in which, however, the original valence p shell has been decontracted to provide an independent function for the empty 5p orbital, which was shown to be important for accuracy. 37 As a first extension, in basis set A2, one set of 4 f polarization functions was added with an exponent of 1.472, as suggested by Ehlers et al. 38 In basis set A3, the cc-aug-pVDZ basis set ͑double-͒ for C and H is replaced by cc-aug-pVTZ ͑triple-͒, and for Pd the LANL2DZ basis set of double-quality is replaced by the LANL2TZ basis set of triple-quality, with the same primitives but further decontracted, which according to Torrent and co-workers, 39 leads to an increased accuracy. Finally, the largest basis set in this series, A4, was created by substituting the single set of 4 f polarization functions of Ehlers et al. 38 by four sets of 4 f polarization functions, as reported by Langhoff and co-workers, 40 with exponents 3.611 217, 1.295 41, 0.554 71, and 0.237 53. They were contracted as 211, resulting in three contracted 4 f functions.
The series B1-B6 in the DIRAC calculations are based on an uncontracted, Gaussian-type basis set (24s16p13d) for palladium, which is of triple-quality, and has been developed by Faegri. 41 Furthermore, Dunning's 35 cc-aug-pVDZ and cc-aug-pVTZ basis sets for C and H were used in uncontracted form because it is technically difficult to use contracted basis sets in the kinetic balance procedure in DIRAC. 42 Basis set B1 corresponds to cc-aug-pVDZ for C and H and the (24s16p13d) basis set for Pd. As a first extension, in basis set B2, one set of 4 f polarization functions was added with an exponent of 1.472 as reported by Ehlers et al. 38 In basis set B3, the cc-aug-pVDZ basis set for C and H is replaced by cc-aug-pVTZ. In basis set B4, the single set of 4 f polarization functions of Ehlers et al. 38 was substituted by four sets of 4 f polarization functions as reported by Langhoff and co-workers 40 with exponents 3.611 217, 1.295 41, 0.554 71, and 0.237 53 that, at variance with the situation in basis set A4, were kept uncontracted. Thereafter, going to basis set B5 an additional set of diffuse p functions was introduced with exponent 0.141 196, as proposed by Osanai et al. 43 Finally, the largest basis set of this series, B6, was created by adding a set of g functions, with an exponent of 1.031 690 071. This value is close to but not exactly equal to the exponent of the g functions optimized by Osanai. Instead, it is equal to the value of one of the exponents of the d set of Faegri, which reduces the computational costs.

A. Geometries of stationary points
First, we examine the stationary points along the reaction coordinate of the oxidative insertion of Pd into the C-H bond of methane. The geometries of these species were fully optimized at the ZORA-BLYP/TZ2P level of relativistic DFT ͑where ZORA stands for zeroth-order regular approximation͒ and agree well with earlier relativistic DFT studies ͑see Fig.  1͒. 15,16 The reaction proceeds from the reactants via formation of a stable reactant complex of C 2v symmetry, in which methane coordinates in an 2 fashion to Pd, followed by the transition state of C s symmetry, and finally a stable product, Completely uncontracted.
also of C s symmetry. All species have been verified through a vibrational analysis to represent equilibrium structures ͑no imaginary frequencies͒ or a transition state ͑one imaginary frequency: 778i cm Ϫ1 ). Thus, we have a set of geometries that, for all stationary points along the reaction coordinate, have been optimized consistently at the same level of relativistic DFT ͑at ZORA-BLYP/TZ2P͒ without any structural or symmetry constraint. In the following, these geometries are used in the series of high-level ab initio calculations that constitute our benchmark study of the PES for the oxidative addition reaction of PdϩCH 4 .

B. Energies of stationary points
As pointed out in the Introduction, the relative energies of stationary points along the reaction profile of Pd insertion into the methane C-H bond, especially the activation energy, appear to be highly sensitive to the level of theory used, as witnessed by the large spread in values computed earlier ͑see Table I͒. Here, we report the first systematic investigation of the extent to which the various thermodynamic and kinetic parameters depend on the quality of the method and the basis set as well as the extent to which these values are converged at the highest level of theory used. The results of our ab initio computations are collected in Tables III and IV ͑relative energies and BSSE͒ and graphically displayed in Figs. 2-5 ͑reaction profiles and BSSE͒. Table S.I in the supplementary material 44 shows the total energies in a.u. of all species occurring at the stationary points as well as the total energies of the corresponding Pd and methane fragments, with and without the presence of the other fragment as ghost. In this way, we can calculate the BSSE and carry out a CPC. In Table S.I, one can note a difference of approximately 5000 a.u. in the total energies of the palladium atom computed with GAUSSIAN and DIRAC. The origin of this difference is the use of an ECP in the GAUSSIAN calculations whereas all-electron ͑SFDC͒ calculations have been carried out with DIRAC.
We proceed with examining the reaction profiles for oxidative insertion of Pd into the methane C-H bond, that is, the energies of the stationary points relative to the reactants Pd ϩCH 4 , which are collected in Table III and, for CCSD͑T͒ and B3LYP, displayed in Figs. 2 and 3. At all levels of theory except HF, the reaction profiles are characterized by the formation of a stable reactant complex ͑RC͒, which leads via the transition state for insertion ͑TS͒ to the oxidative-addition product ͑P͒. Three striking observations can be made: ͑i͒ the spread in values of computed relative energies, depending on the level of theory and basis set, is enormous, up to nearly 70 kcal/mol; ͑ii͒ the size of the BSSE is also remarkably large, up to ca. 30 kcal/mol; ͑iii͒ most strikingly, convergence with basis-set size of the computed energies is still not reached with standard basis sets used routinely in CCSD͑T͒ computations on organometallic and coordination compounds. The lack of any correlation leads to a complete failure at the HF level, which yields an unbound RC, a strongly exaggerated activation barrier of ca. 45 kcal/mol, and a reaction energy that differs by only a few kcal/mol from the activation energy. In other words, the process is highly endothermic and has essentially no reverse barrier at the HF level for all basis sets used. The failure of HF for describing the PES of our model reaction is not unexpected because electron correlation, which is not contained in this approach, is important. 45,46 The activation energy drops significantly when electron correlation is introduced. Along HF, CCSD and CCSD͑T͒ in combination with basis set A1, for example, the activation barrier decreases from 45.0 to 7.7 to 3.6 kcal/mol. But also the correlated CCSD͑T͒ values obtained with standard basis sets, such as LANL2DZ or LANL2TZ with one or three f functions ͑A2-A4 in Table II͒ are questionable, as they are obviously not converged as a function of the basis-set size. For example, the activation energy of 3.6 kcal/mol at CCSD͑T͒/A1, which involves LANL2DZ for Pd, agrees exactly, that is, it coincides with that of the present benchmark of Siegbahn and coworkers obtained with PCI-80 ͑see Sec. I͒. This agreement is fortuitous. The CCSD͑T͒ value for the activation energy drops from 3.6 to Ϫ3.6 and further to Ϫ18.8 kcal/mol along basis sets A1, A2, and A3, as one f polarization function is added and, then, the flexibility of the basis set is increased from double-to triple-. Thereafter, going from basis set A3 to A4, the activation energy increases again from Ϫ18.8 to Ϫ9.4 kcal/mol as polarization is increased from one to three f functions ͑see Tables II and III͒. This is illustrated by Fig. 2, which shows the CCSD͑T͒ reaction profiles and how they vary along basis sets A1-A4. Obviously, the energy values have not reached convergence. Also, an activation energy of Ϫ18.8 or Ϫ9.4 kcal/mol is not only much lower than the present benchmark value of 3.6 kcal/mol but it is also too low for a reaction that essentially does not proceed. Similar behavior is observed for other correlated ab initio methods ͑MP2, MP4SDQ, CCSD͒ both in the relativistic ECP calculations with GAUSSIAN with basis sets A1-A4 and in the SFDC calculations with DIRAC with basis sets B1-B4 ͑see Table III͒. On the other hand, the HF calculations in which electron correlation is not accounted for are relatively insensitive toward increasing the flexibility and polarization of the basis set along A1-A4 or B1-B4 ͑see Table III͒.
Next, we note that the BSSE takes on large values in the correlated ab initio methods, whereas it is negligible if correlation is completely neglected, i.e., in HF ͑see Fig. 4͒. The BSSE values are summarized in Table IV and     sets B1-B4 are smaller than those obtained with A1-A4, but they display a similar trend as the latter ͑compare Figs. 4 and 5͒. The BSSE increases along the reaction coordinate, i.e., going from RC to TS to P ͑Figs. 4 and 5͒. The reason for this is that along this series of stationary points, the carbon and hydrogen atoms and, thus, their basis functions come closer too and begin to surround the palladium atom. This effectively improves the flexibility and polarization of the basis set and thus the description of the wave function in the region of the palladium atom. Note also that the BSSE stems nearly entirely from the improvement of the stabilization of palladium as methane ghost functions are added ͑Table IV͒. The energy lowering of methane due to adding palladium ghost functions is in all cases small, i.e., less than 1 kcal/mol. The above points out the prominent role that electron correlation plays in our model systems. And, more importantly, it also reveals the inadequacy of basis sets A1-A4 and B1-B4 for describing it. This may be somewhat surprising in view of earlier reports that basis sets of a quality comparable to that of A3, A4 and B3, B4 yield satisfactory energies for organometallic and coordination compounds ͑see, for example, the excellent reviews by Frenking et al. 45 and by Cundari and co-workers 46 ͒. On the other hand, it is consis-tent with the large variation of values for the thermodynamic and kinetic parameters obtained in earlier theoretical studies on the present model reaction ͑see Table I͒. One reason for the increased sensitivity that we find toward the quality of the theoretical approach is that the presence of f polarization functions is only the minimum requirement for describing the electron correlation of palladium 4d electrons. In this respect, the palladium basis sets in A3, A4 and B3, B4 should be considered minimal and cannot be expected to have achieved convergence. Furthermore, the consequences of any inadequacy in the basis set shows up more severely in processes, such as ours, which involve a bare, uncoordinated transition metal atom as one of the reactants because here the effect of the additional assistance of basis functions on the substrate is more severe than in situations where the transition metal fragment is already surrounded by, e.g., ligands. This shows up in the large BSSE values. Note that, in line with this, the energies of the transition state and the product vary much less along the basis set if they are computed relative to the reactant complex ͑RC͒ instead of the separate reactants ͑R͒, see Table III and Figs. 2  atoms than Pd changes much less than along R ͑no assis-tance͒ to RC ͑significant assistance͒.
We have been able to achieve virtual convergence of the CCSD͑T͒ relative energies by further increasing the flexibility and polarization of the Pd basis set and by correcting for the BSSE through CPC; see Table III and Figs. 2 and 3. Let us first point out why the CCSD͑T͒/A4 and CCSD͑T͒/B4 values cannot be trusted without further scrutiny. This is an important issue because inspection of Table III and Figs. 2 and 3 suggests that the counterpoise-corrected energies do converge from A3 to A4 and from B3 to B4. For example, the counterpoise-corrected activation energies computed with A3 and A4 at CCSD͑T͒ are equal within less than 2 kcal/mol. Note however that the BSSE of 14 -20 kcal/mol is still larger than the effect we wish to compute. It is therefore necessary to explore the behavior of the reaction profile if the basis set is further increased. In particular, we wish to achieve a situation in which the BSSE at least becomes smaller than the effects that are to be computed, i.e., the relative energies. Thus, we have introduced an additional diffuse p function ͑going from B4 to B5͒ and a g polarization function ͑going from B5 to B6͒. We have chosen the B series of basis sets ͑based on Faegri's 24s16p13d basis for Pd͒ for further improvements because they are superior to the A series ͑based on Hay and Wadt's LANL2TZ basis for Pd͒ in the sense that they yield a significantly smaller BSSE ͑compare A1-A4 with B1-B4 in Figs. 4 and 5͒. Indeed, along B3-B6, the BSSE in, for example, the CCSD͑T͒ activation energy decreases monotonically from 19.8 to 6.2 to 4.0 to 2.7 kcal/ mol and is thus clearly smaller than the relative energies that we compute ͑see Table III and Fig. 5͒. The counterpoisecorrected relative energies at CCSD͑T͒ are converged within a few tenths of kcal/mol along B1-B6. For example, the counterpoise-corrected activation energy at CCSD͑T͒ amounts to 14.7, 10.3, 6.9, 6.1, 5.9, and 5.8 kcal/mol. Our best estimate, obtained at CCSD͑T͒/B6 with CPC, for the kinetic and thermodynamic parameters of the oxidative insertion of Pd into the methane C-H bond is Ϫ8.1 kcal/mol for the formation of the reactant complex, 5.8 kcal/mol for the activation energy relative to the separate reactants, and 0.8 kcal/mol for the reaction energy. If we take into account zero-point vibrational energy ͑ZPE͒ effects computed at BLYP/TZ2P, this yields Ϫ8.3 kcal/mol for the formation of the reactant complex, 3.1 kcal/mol for the activation energy relative to the separate reactants, and Ϫ1.2 kcal/mol for the reaction energy. Our benchmark values agree well with those obtained by Siegbahn and co-workers at PCI-80, and therefore further consolidate the theoretical reaction profile. They also agree well with the experimental result, in fact slightly better so than PCI-80, that the reactant complex is bound by at least 8 kcal/mol. The fact that the experimental reaction rate is extremely small to negligible in spite of a moderate energy barrier of 3.1 kcal/mol is consistent with an important statistical or entropic bottleneck 15 ͑associated with the decrease in the number of available quantum states as one goes from reactants to transition state͒ and the extremely short lifetime of the internally hot reactant complex that has been invoked to explain why this complex has not been observed in experiments. 6 Finally, we note that the BSSE is small not only in uncorrelated ab initio calculations ͑HF͒ but also in the DFT calculations ͑B3LYP͒. This robustness of DFT is due to the way in which the correlation hole is described in this method rather than to the absence of correlation as in HF. In general, correlated ab initio methods depend more strongly on the extent of polarization of the basis set because the polarization functions are essential to generate the configurations through which the wave function can describe the correlation hole. In DFT, on the other hand, the correlation hole is built-in into the potential and the energy functional and po-larization functions mainly play the much less delicate role of describing polarization of the electron density. Interestingly, in the HF and B3LYP calculations with GAUSSIAN, we observe a small but non-negligible BSSE of 2-4 kcal/mol, which does not decrease with increasing basis-set size along A1-A4 ͑see Table IV and Fig. 4͒. In the DIRAC calculations, however, the BSSE for both HF and B3LYP is essentially zero ͑less than 0.2 kcal/mol͒ for all basis sets B1-B6. This difference between the GAUSSIAN and DIRAC calculations is ascribed to the fact that an ECP for palladium is used in the former whereas the latter are all-electron calculations. Ideally, the ECP should account for the fact that the valence orbitals must be orthogonal to the core orbitals. It is likely however that, effectively, the ECP used with the LANL2 basis sets of Pd is not able to completely project out the palladium-core components of the methane orbitals. In the DIRAC all-electron calculations this problem is of course not present. We conclude that the best B3LYP reaction profile with an activation energy of 9.3 kcal/mol is obtained at the relativistic SFDC-B3LYP/B6 level with counterpoise correction.

IV. CONCLUSIONS
We have computed a benchmark for the archetypal oxidative addition of the methane C-H bond to palladium that FIG. 4. Basis set superposition errors ͑BSSE͒ for stationary points along the reaction coordinate of oxidative insertion of PdϩCH 4 ͑in terms of these reactants͒ obtained using GAUSSIAN with various methods and basis sets. derives from a hierarchical series of relativistic ab initio methods and highly polarized basis sets. Our best estimate of kinetic and thermodynamic parameters is Ϫ8.1 kcal/mol for the formation of the reactant complex, 5.8 kcal/mol for the activation energy relative to the separate reactants, and 0.8 kcal/mol for the reaction energy. This is obtained at the counterpoise-corrected, four-component spin-free Dirac-Coulomb CCSD͑T)/(24s16p13dϩ4 f ϩpϩg) level, which is virtually converged with respect to the basis-set size.
Our benchmark values agree well with those obtained by Siegbahn and co-workers with the parameterized PCI-80 method and slightly better than the latter with experimental data. The agreement and, importantly, the fact that our CCSD͑T͒ benchmark PES derives from a converged hierarchical series of basis sets consolidates the theoretical PES for oxidative addition of CH 4 ϩPd.
Our findings stress the importance of sufficient higherangular momentum polarization functions, f and g, as well as counterpoise correction for obtaining reliable activation energies. We show that standard basis sets, such as LANL2DZϩ1 f for palladium, are not sufficiently polarized for this purpose and lead to erroneous results at CCSD͑T͒.