12716
J. Phys. Chem. A 2010, 114, 12716–12724
Homogeneous Ni Catalysts for H2 Oxidation and Production: An Assessment of Theoretical Methods, from Density Functional Theory to Post Hartree-Fock Correlated Wave-Function Theory Shentan Chen, Simone Raugei, Roger Rousseau, Michel Dupuis,* and R. Morris Bullock Center for Molecular Electrocatalysis, Chemical and Materials Sciences DiVision Pacific Northwest National Laboratory, Richland WA ReceiVed: July 21, 2010; ReVised Manuscript ReceiVed: September 23, 2010
A systematic assessment of theoretical methods applicable to the accurate characterization of catalytic cycles of homogeneous catalysts for H2 oxidation and evolution is reported. The key elementary steps involve heterolytic cleavage of the H-H bond and formation/cleavage of Ni-H and N-H bonds. In the context of density functional theory (DFT), we investigated the use of functionals in the generalized gradient approximation (GGA) as well as hybrid functionals. We compared the results with wave-function theories based on perturbation theory (MP2 and MP4) and on coupled-cluster expansions [CCD, CCSD, and CCSD(T)]. Our findings indicate that DFT results based on Perdew correlation functionals are in semiquantitative agreement with the CCSD(T) results, with deviations of only a few kilocalories/mole. On the other hand, the B3LYP functional is not even in qualitative agreement with CCSD(T). Surprisingly, the MP2 results are found to be extremely poor, in particular for the diproton Ni(0) and dihydride Ni(IV) species on the reaction potential energy surface. The Hartree-Fock reference wave function in MP2 theory gives a poor reference state description for these states that are electron rich on Ni, giving rise to erroneous MP2 energies. We present a detailed potential-energy diagram for the oxidation of H2 by these catalysts after accounting for the effects of solvation, as modeled by a polarizable continuum, and of free energy estimated at the harmonic level of theory. Introduction There is great current interest in catalysts that utilize inexpensive, earth-abundant metals1 for the interconversion of electrical and chemical energy, such as the electrochemical oxidation and production of H2 (Scheme 1) and reduction of O2. A biomimetic approach to the design of catalysts relies on the understanding of the specific chemical functionalities that are incorporated in enzymes and how the functionalities contribute to their mechanism and activity. Such an understanding offers a path to designing synthetic catalysts that perform the same reactions. In this regard, the recent elucidation of the structures of several hydrogenase enzymes with bimetallic iron or iron-nickel complexes at the active sites2-10 has provided a strong impetus for the development of inexpensive synthetic catalysts for hydrogen oxidation and production. Highly efficient electrocatalysts for both hydrogen production and hydrogen oxidation have been recently designed and synthesized.11,12 The catalyst precursors are air-stable, mononuclear nickel(II) complexes containing cyclic diphosphine ligands with nitrogen bases incorporated into the ligand (see structure I). A key feature for obtaining the high catalytic rates observed with these complexes has been shown to be the presence and positioning of the nitrogen bases near the metal center, because these bases facilitate the heterolytic cleavage or formation of the H-H bond. These pendant amines function as proton relays, accelerating intramolecular and intermolecular proton-transfer reactions, but they also
have additional roles, including facilitating the heterolytic cleavage of hydrogen and proton-coupled electron-transfer reactions.
An important step toward the design of more efficient catalysts lies in the characterization, through experimental and theoretical investigations, of the chemical properties (metal hydricities, pKa values, redox potentials) of such complexes and how variations in these properties affect the catalytic efficiency. Systematic computational characterizations of such properties have been recently reported.13-16 A complete characterization of all the intermediates in their catalytic cycle is also a key to understanding the structure/activity relationship for these comSCHEME 1
* Corresponding author.
10.1021/jp106800n 2010 American Chemical Society Published on Web 11/11/2010
Homogeneous Ni Catalysts for H2 Oxidation and Production plexes, in particular providing insights about the energy landscape associated with the elementary steps of the catalytic cycle.17,18 All of this earlier work was based on density functional theory (DFT) calculations that account for electron correlation and thus is broadly believed to be accurate, at a reasonable computing cost for even large organometallic complexes. The class of catalysts considered here is denoted Ni(PR2NR’2)22+ where P2N2 denotes 1,5-diaza-3,7-diphosphacyclooctane ligands with various substituent groups R and R’ covalently bound to each of the phosphorus and nitrogen atoms, respectively. When R and R’ are phenyl groups, Xray diffraction studies have shown that the cation of [Ni(PPh2NPh2)2(CH3CN)](BF4)2, (where PPh2NPh2 is 1,3,5,7tetraphenyl-1,5-diaza-3,7-diphosphacyclooctane) is a trigonal bipyramid with bonds to four phosphorus atoms of the two bidentate diphosphine ligands and to the nitrogen atom of an acetonitrile molecule, which is the most commonly employed solvent for studying the chemistry of these catalysts. Two of the six-membered rings formed by the diphosphine ligands and Ni have boat conformations with an average Ni-N distance to the two pendant bases of ∼3.4 Å. When R ) cyclohexyl and R’ ) benzyl, the cation in [Ni(PCy2NBz2)2](BF4)2, (where Cy ) cyclohexyl and Bz ) benzyl) is a distorted square planar complex. One conformation associated with the four sixmembered rings formed upon coordination of the diphosphine ligands to the metal has them in a boat form. In this case, the average Ni-N distance to the pendant base is ∼3.3 Å. [Ni(PPh2NPh2)2(CH3CN)]2+ is an electrocatalyst for hydrogen production in acidic acetonitrile solutions.17,18 The closely related complex [Ni(PCy2NBz2)2]2+ is an electrocatalyst for hydrogen oxidation in basic acetonitrile solutions.17 When considering that modifications to the groups R and R’ have profound impact on the overall catalytic activity, it is necessary that theoretical studies account for the effects of substitutions, a requirement which in turn implies that large numbers of atoms be explicitly included in the models. To date, this can only be practically achieved via DFT studies. Hence, it is of paramount importance to gauge how reliable DFT may be in accounting for experimentally observable quantities. This is particularly true for the part of the reaction cycle that is known to encompass the ratelimiting steps for both H2 oxidation and production, specifically the transformation of H2 into two protons, which reside on the pendant amines.11,12 In the present paper, we provide a complete characterization of the initial stages of H2 oxidation by using several levels of theories, for a model catalyst. We selected R and R’ to be methyl groups. Such a model has strong similarities with several actual catalysts and allows us to carry out benchmark calculations by using electron-correlated wave-function theories, including coupled-cluster theory with single and double excitations and perturbative treatment of triple excitations CCSD(T). We believe that this model system has many of the important features like those of the real catalysts; yet, it is small enough to be amenable to high-level electron-correlation treatements. We also report results for lower levels of theory (Moller-Plesset perturbation theory up to second and fourth order, including contributions from triple excitations). These comparative studies are aimed at establishing a baseline level of reliability, in particular from the point of view of DFT. We also assess the effects of solvation and entropy on the stability of intermediates along the reactions path. A comparison with a more realistic model of an actual catalyst (R ) cyclohexyl, R’ ) Me)19 is also provided.
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12717 Before discussing the details of our investigation, we summarize here the most important findings. The hydrogen molecule forms initially a dihydrogen20 complex with the nickel Ni(II) site, that splits H2 readily, forming a proton/hydride Ni(II) intermediate (with a protonated amine), which further evolves into a diproton Ni(0) species with two protonated amines. First, the heterolytic character of the cleavage of the H-H bond and the EPR silent character of the catalysts18 suggest that the process ought to be well described by a spin-restricted singletspin closed-shell determinant (Slater or Kohn-Sham) and that DFT should describe well the electron-correlation effects along the reaction path. Our findings show that commonly used functionals such as PBE21,22 and B3P8623,24 results are in qualitatively accord and give results in agreement with CCSD(T). However, there is a notable exception: the B3LYP23,25,26 hybrid functional (often the de facto choice for these types of studies in organometallic chemistry) predicts the reaction to be slightly endothermic, whereas the more accurate methods show the reaction to be strongly exothermic. Second, the ionic character of the bond-breaking process suggests that the selfinteraction error often associated with unpaired electrons and often seen in the description of bond breaking might be small, and therefore, hybrid functionals (which include a fraction of the Hartree-Fock (HF) exchange) might not offer substantial accuracy advantages. This is borne out by our present results. Third, MP2 theory and DFT are often believed to offer a comparable level of accuracy, even though exceptions have been reported,27 notably when involving 3d transition metals. The present case offers a catastrophic breakdown of this rule because MP2 is found to be strikingly incorrect. Not surprisingly, this observation applies as a start to the HF level of theory that yields a diproton product much higher (+60 kcal/mol) in energy than the dihydrogen reactant. In contrast, MP2 gives a highly exothermic process (-110 kcal/mol). This is in contrast with CCSD(T), which yields a small (-7.2 kcal/mol) exothermicity for the overall reaction. The catastrophic behavior of MP2 is seen for the dihydride Ni(IV) species and the diproton Ni(0) species; both of these states are electron rich at the Ni site. The latter state is formally a d10 configuration of the Ni atom, which is poorly described by the reference closed-shell HF single determinant. The technical details related to these calculations are given in Section II of this paper. Section III reports and discusses the results obtained for the model system and that we believe to be applicable to this class of catalysts and other catalysts that show similarity in the structure of the ligands. We conclude in Section IV. Methods Overview. Correlated wave-function and DFT calculations were carried out to characterize the reactant, products, intermediates, and transition states for the reaction of H2 with the cation of [Ni(PMe2NMe2)2]2+ (where PMe2NMe2 is 1,3,5,7-tetramethyl-1,5-diaza-3,7-diphosphacyclooctane). Experimentally the groups on the phosphorus nitrogen atoms are different (cyclohexyl, benzyl, phenyl). However, using methyl is computationally more tractable, given our interest in carrying out CCSD(T) calculations. The all-methyl complex possesses the essential features of the actual catalysts, a four-coordinate Ni2+ center with chelating diphosphines, each with two pendant amines. These structural features are known to facilitate heterolytic H-H bond breaking or forming,11,12 and hence, this system should be a satisfactory model of the real catalysts for the purposes of
12718
J. Phys. Chem. A, Vol. 114, No. 48, 2010
Chen et al. depicted in Figure 2, and their key structural parameters given in the Supporting Information (Table S1). The coordinates for all the intermediates and transitions states are available in the Supporting Information, along with tables summarizing the geometrical parameters obtained at different levels of theory. Structures of the intermediates and transition states are barely affected by the choice of computational protocol. In particular, the various methods discussed here are able to predict the short distance between the metal center and the pendant amine N atom (differences within 0.1 Å). Therefore, our discussion will be limited to only the analysis of the relative energy of the species involved in early stages of H2 oxidation.
Figure 1. Energy diagram for the H2 oxidation by Ni(PMe2NMe2)22+ calculated at the B3P86/bs1 level of theory.
the current study, which is to gauge the accuracy of the treatment of the electronic correlations by DFT methodologies. We determined the structures of a series of stable intermediates suggested by the experimental understanding of the H2 oxidation process: (i) a reactant structure made of the Ni(II) dication complex and an isolated H2 molecule (structure a); (ii) a Ni(II) dihydrogen dication adduct with the H2 molecule forming a long-range complex with the Ni dication (structure b); (iii) a Ni(IV) dihydride dication species with the two hydride H atoms bound to the Ni atom (structure c); (iv) a Ni(II) hydrideproton dication species with one hydride H bound to Ni and a proton H attached to an amine (structure d); (v) a Ni(0) diproton dication species with two protons bound to two amines (structure e). In addition, transition-states structures for the transformation of b to c, b to d, c to d, and d to e were also determined. They are denoted TS(bc), TS(bd), TS(cd), and TS(de), respectively. The relative energy diagram obtained at the B3P86/bs1 (see below) level of theory is shown in Figure 1 with structures
We note that, in this particular case of H2 oxidation, the ratedetermining step is the formation of the hydride-proton intermediate d. The dihydride intermediate is significantly higher in energy. The nature of the potential energy surface changes with the nature of the ligands,28 although the heterolytic character of bond cleavage or formation is retained in all cases, and so is the closed-shell (paired electrons) nature of the electronic wave function. Basis Set Selection. We assessed the effect of basis-set quality and size. Our aim was (i) to establish to what extent the basis set affects the reaction energy profile and (ii) to what extent structural parameters change with the basis set and finally (iii) to assess whether the nonessential groups (the -CH2- groups in the six-membered rings and the substituent groups R and R’ on P- and N-methyl, benzyl, phenyl, cyclohexyl-) require an extended basis set like the phosphine and amine bases, which are expected to play a crucial role in the H2 oxidation or production reactions along with the metal center. Thus, several basis sets were used. 1. The correlation-consistent series of basis sets cc-pVnZ (n ) D, T, Q) for all atoms, recently extended to transition-
Figure 2. Structures of minima and transition states involved in the H2 oxidation by Ni(PMe2NMe2)22+. All H atoms are removed for simplicity except those from the reactant H2. Color code is as follows: Ni (light blue), P (orange), N (dark blue), C (gray), H (purple).
Homogeneous Ni Catalysts for H2 Oxidation and Production
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12719
TABLE 1: Influence of the Basis on the Relative Energy (kcal/mol) of the Various Species Involved in the H2 Chemistrya a b TS(bc) c TS(cd) TS(bd) d TS(de) e ME
bs1
bs2
bs3
bs4
6-31G**
6-311G**
cc-pVDZ
cc-pVTZ
cc-pVQZ
CBS
0.0 -3.6 6.1 4.9 5.4 0.9 -7.8 -2.7 -10.7 -1.9
0.0 -2.5 8.6 7.9 8.6 3.2 -4.9 1.7 -5.2 1.3
0.0 -2.6 8.3 6.9 6.4 1.0 -7.9 -4.2 -12.5 -1.5
0.0 -4.5 5.5 4.2 4.9 0.6 -8.1 -3.0 -11.0 -2.3
0.0 -5.9 0.4 0.4 1.8 -1.1 -9.1 -1.2 -7.2 -3.7
0.0 -3.1 7.4 6.8 7.7 2.8 -4.6 3.0 -2.7 1.3
0.0 -4.2 4.8 4.0 4.7 0.1 -7.8 -2.0 -8.5 -2.0
0.0 -3.1 7.6 6.8 7.5 2.5 -6.5 -0.7 -8.7 -0.24
0.0 -2.7 7.8 7.3 8.0 2.5 -6.4 -0.5 -8.8 -0.01
0.0 -2.7 7.8 7.3 8.1 2.4 -6.4 -0.4 -8.8 0.0
a CBS is the complete basis-set limit energy extrapolated from the cc-pVnZ series. ME is the mean error with respect to the estimated complete basis-set limit. All of the calculations were performed on structures optimized at the DFT(B3P86)/bs1 level of theory.
metal atoms;29-31 using these basis sets the complete basisset limit was estimated with the three-point extrapolation formula.32 2. The 6-31G** basis set for all atoms, the double-ζ valence polarized basis set, very well established for second- and third-row atoms, with more recent extensions to transitionmetal atoms.33 3. The 6-311G** basis set for all the atoms.34,35 4. bs1, comprising the Stuttgart effective core SDD basis set for Ni with an effective potential for the core electrons36 (three d functions and one f function), combined with the 6-31G* basis set for all the ligand atoms and the 6-31G** for H2 (that includes p polarization functions on H atoms). 5. bs2, comprising of Ahlrichs’ triple-ζ VTZ basis set for Ni37 (three d functions) and the 6-311G** basis set for all other atoms, including polarization d and p functions for C, N, and P and for H, respectively. 6. bs3, which resembles bs2, but the nonessential atoms (C and H) were given a double-ζ quality basis set 3-21G, with the essential atoms (P, N, and H from H2) being given a 6-31G** basis set. 7. bs4, made up of Wachter’s basis augmented with d and f basis functions for Ni35,38 (four d functions and one f function) and the 6-311G** basis set for all other atoms (bs4 is identical to bs2 except for the Ni basis set). DFT and Post-HF Method Selections. The starting level of theory used in this work was DFT(B3P86)23,24 in conjunction with basis set bs1. The choice of the functional and basis set was based on the methodology published by Qi et al.13 who demonstrated an excellent level of accuracy for a variety of properties for a large number of Ni complexes as well as complexes involving other transition metals. The geometries of the stable intermediates and of the transition states (mentioned above) were initially optimized at this level of theory. We tested other density functionals including B3LYP,23,25,26 BLYP,25,26,39 ans PBE21,22 and the more recent BMK,40 M05,41 and M0642 functionals. The favorable computational cost for generalized gradient approximation (GGA) functionals, such as PBE, makes them attractive for example for ab initio molecular dynamics calculations to characterize the dynamics of the protons. As we will see below, the various DFT methods do not provide a consistent picture of the relative energetics of the various species considered in this work. Thus, it became necessary to carry out a set of reference calculations by using the highly correlated molecular orbital theory method CCSD(T), to establish a metric by which to gauge the reliability of the DFT schemes. We also considered the Moeller-Plesset secondand fourth-order perturbation theory methods MP2 and MP4. The calculations in this work were performed with the NWCHEM software43,44 and Gaussian 09.45 Basis sets which
are not available in the NWCHEM basis set library were obtained from the Environmental Molecular Science Laboratory (EMSL) basis-set library.46,47 Results and Discussion Assessment of Basis Sets. To establish a protocol for assessing the accuracy of the various electronic-structure methods, it is first necessary to establish the dependence of our results upon the basis set. Specifically, we aim to determine a minimum level of basis set which will ensure a consistent level of accuracy with respect to results obtained for a larger and more complete basis set. We determined the optimized structures for the various states along the reaction paths by using the DFT(B3P86)/bs1 level of theory. This is our base level of theory because it has been demonstrated to yield good thermodynamic quantities such as pKa, hydricities, and redox potentials for a series of similar compounds.13 By using these structures, we calculated singlepoint energies of the species by using the B3P86 functionals and various basis sets. The results are presented in Table 1. The best basis set used here is the correlation-consistent polarized quadruple-ζ valence basis set cc-pVQZ, which has been widely benchmarked in electron correlated thermodynamics of small molecules.48,49 As indicated above and most importantly, the cc-pVnZ basis sets include a highly accurate correlation consistent basis set for Ni. Indeed, both the cc-pVTZ and cc-pVQZ basis sets relative energies are very close to the complete basis-set limit extrapolated from the cc-pVnZ series. Among the other adopted basis sets, bs1 and bs4 give relative energies consistently close to the ones from the cc-pVQZ basis set when considering structures b, d, and e because those display similar stabilities within a little less than 2 kcal/mol. Other basis sets are somewhat more erratic in reproducing the relative stabilities of the various intermediates and transition states. In particular, the 6-311G** and bs2 basis sets underestimate significantly the stability of the diproton structure e. bs2 shares the same basis as bs4 for all the nonmetal atoms; therefore, we can assign the deficiency of the 6-311G** set to the atomic basis for Ni. Likewise, the 6-31G** basis overestimates the stability of d and underestimates the energy barriers for c to d and b to d. Given that the difference between bs1 and 6-31G** resides essentially in the description of the Ni atom, it can be again concluded that the 6-31G** basis set provides an unsatisfactory description of the metal center. Overall, these results confirm the quality of bs1 that yield good thermodynamics properties for a set of similar catalysts.12 Comparing bs1 and bs3 suggests that, not so surprisingly, the quality of the basis set on the nonessential atoms is not critical for the study of relative energies of stationary points on the potential energy surface.
12720
J. Phys. Chem. A, Vol. 114, No. 48, 2010
Chen et al.
TABLE 2: Relative Energies (kcal/mol) Calculated at DFT(B3P86) Level for Different Basis Setsa
a
structure
bs1
bs2
bs3
a b TS(bc) c TS(cd) TS(bd) d TS(de) e
0.0 -3.6 6.1 4.9 5.4 0.9 -7.8 -2.7 -10.8
0.0 -2.5 8.2 7.8 8.4 3.1 -5.1 1.6 -5.2
0.0 -2.5 8.3 6.5 6.7 1.2 -7.9 -3.1 -12.2
TABLE 4: Relatives Energies (kcal/mol) for HF, Post-HF, and DFT(B3P86) Methods Obtained by Using Basis set bs3a HF a b TS(bc) c TS(cd) TS(bd) d TS(de) e
MP2
MP4
CCD
CCSD CCSD(T) B3P86
0.0 0.0 0.0 0.0 21.6 -7.3 -4.2 4.1 69.6 -37.9 13.3 13.1 81.6 -60.5 15.7 9.5 78.5 -47.7 11.3 14.3 33.3 0.1 2.2 12.6 22.7 -17.2 -8.5 -0.1 70.2 -59.5 -7.6 7.0 62.0 -110.8 -26.6 -10.6
0.0 1.2 19.7 22.9 21.8 9.1 -2.0 13.1 7.5
0.0 -1.8 12.7 12.2 12.5 5.0 -5.3 2.8 -7.2
0.0 -2.5 8.3 6.6 6.7 1.2 -7.9 -3.1 -12.2
a All values are single-point energy calculations on DFT(B3P86)/ bs3 optimized structures.
The structures were optimized for each basis set.
We reoptimized the structures of all the species by using the DFT(B3P86) level of theory with basis sets bs1, bs2, and bs3. It can be seen in Table 2 that the relative energies obtained with the various basis sets are consistent with the single-point energies, indicating that the effects due to reoptimization of the structures for each basis set are not large, and the reoptimization efforts are not necessarily warranted. Assessment of the Effect of the Functional of the Density in DFT: Comparison with Moeller-Plesset and CoupledCluster Wave-Function Theories. It is of interest to assess the dependency of the relative energies on the functionals of the density in DFT. It is well accepted that hybrid functionals such as B3P86 and B3LYP are more accurate than GGA functionals such as BLYP and PBE in the field of chemistry, albeit at a significantly higher computational cost.50,51 There also remains an uncertainty about the performance of B3LYP when treating transition metals.27,52 The comparison between functionals is reported in Table 3. It can be seen that for the present system, B3P86 and PBE give qualitatively the same relative energies for all of the reaction intermediates and transition states (differences within 2 kcal/ mol), except for the species e and the transition state connecting d to e, TS(de), for which PBE yields an appreciably lower energy (by about 5 and 6 kcal/mol, respectively). M06 is also consistent with B3P86, except for species c which is significantly higher in energy. In contrast, B3LYP, BLYP, BMK, and M05 predict far higher relative energies (in some cases, as high as 10 and 20 kcal/mol for B3LYP and M05, respectively). In particular, B3LYP and BLYP do not yield structure b at all. In short, these functionals yield a qualitatively different potentialenergy landscape than B3P86 and PBE. B3LYP relative energies for the B3P86 optimized structures show the same trends. Different basis sets also yield similar observations. Hence, the differences can be attributed solely to the functional. The consistency between BLYP and B3LYP energies strongly suggests that the root cause of the discrepancy with B3P86 may lie, in part, in the LYP description of electron correlation, which,
for this class of compounds, shows appreciable differences with the correlation functionals derived by Perdew and co-workers (P86 and PBE). Given the discrepancies discussed above among DFT descriptions of this system, a more stringent evaluation of the role of electron correlations in this reaction was required through comparison with electron-correlated molecular orbital theories. It is broadly accepted that DFT(B3LYP) is of an accuracy similar to that of second-order Moeller-Plesset perturbation theory MP2 for nonmetal containing species and that the CCSD(T) coupled-cluster level of theory often approaches chemical accuracy. Thus, it is of great interest to assess the accuracy of the DFT results for this class of organometallic complexes with respect to the MPn theories, namely, MP2 and MP4, and coupled-cluster theories, specifically CCSD(T). The relative energies of the various structures at these levels of theory, calculated with the bs3 basis set, are given in Table 4. We note that, as before, basis set bs3 is qualitatively consistent with larger basis sets; yet, it is small enough to allow MP and CC calculations. However, it must be stressed that, with this basis set, the correlated wave-function theory methods used here are farther away from their intrinsic limits with respect to basisset convergence than the DFT methods are. Therefore, we are looking only at semiquantitative consistency among methods. The B3P86 relative energies are in overall qualitative agreement with the CCSD(T) relative energies. In particular, B3P86 seems to describe well the weakly bound H2 complex b. The largest difference lies in the relative stability of structure e, which B3P86 predicts to be almost twice as stable as CCSD(T). Structure e is the diproton species with a Ni(0) site in a formally d10 configuration with the two protons attached to the amines. When comparing the CCSD and CCSD(T) relative energies, we observe that electron-correlation effects captured by the triple excitations are significant (Table 4). However, these effects do not qualitatively alter the nature of the potential energy surface in that the overall reaction is predicted to be exothermic
TABLE 3: Relative Energies (kcal/mol) for Different DFT Functionalsa a b TS(bc) c TS(cd) TS(bd) d TS(de) e a
B3P86//B3P86
PBE//PBE
0.0 -3.6 6.1 4.9 5.4 0.9 -7.8 -2.7 -10.7
0.0 -4.3 2.0 -1.1 1.8 -1.5 -8.9 -8.5 -15.3
BLYP//BLYP
B3LYP//B3LYP
B3LYP//B3P86
BMK//B3P86
M05//B3P86
M06//B3P86
0.0
0.0
12.4 9.8 10.1 8.2 0.1 2.8 -4.5
16.5 15.5 16.7 10.2 1.0 9.1 0.2
0.0 2.3 14.3 13.3 14.5 8.1 -1.0 6.7 -2.0
0.0 1.4 16.5 17.1 21.1 8.7 -0.6 24.0 21.5
0.0 0.0 18.7 20.5 20.6 9.2 1.5 6.8 -3.9
0.0 -7.5 9.1 11.8 13.0 2.3 -6.1 2.4 -8.2
All calculations performed by using the bs1 basis set.
Homogeneous Ni Catalysts for H2 Oxidation and Production TABLE 5: Lowering of Interaction Energy (in kcal/mol) from Electron Correlation Calculated by Using Basis Set bs3a ∆(CCSD(T)) ∆(CCSD) ∆(CCD) ∆ (MP4) ∆ (MP2) POP a b c d e
0.00 -23.4 -69.4 -28.0 -69.2
0.00 -20.4 -58.7 -24.7 -54.5
-17.5 -20.7 -72.1 -22.8 -72.6
0.00 -25.8 -65.9 -31.2 -88.6
0.00 -28.9 -142.1 -39.9 -172.8
8.64 8.68 8.83 8.74 9.35
a All values are single-point energy calculations on DFT(B3P86)/ bs3 optimized structures. The Ni d-orbital population (POP) based on NBO analysis performed by using HF/bs3 is also reported.
TABLE 6: HOMO and LUMO Energies and HOMO-LUMO Gap (in a.u) of the HF Wave Function of the Various Statesa a a TS(bc) c TS(bd) TS(cd) d TS(de) e
ε(HOMO)
ε(LUMO)
∆ε(HOMO-LUMO)
-0.5758 -0.5756 -0.5677 -0.5400 -0.5628 -0.5286 -0.5483 -0.4797 -0.4773
-0.1318 -0.0978 -0.1140 -0.1171 -0.0957 -0.1074 -0.0907 -0.0898 -0.0670
0.4440 0.4778 0.4537 0.4229 0.4671 0.4212 0.4576 0.3899 0.4103
a All values are single-point energy calculations on DFT(B3P86)/ bs3 optimized structures.
by 7 kcal/mol. We remark that M06 seems to describe better the high-energy dihydride c species as seen in a comparison of Tables 3 and 4. It is not surprising that the HF relative energies are the least accurate, because correlation effects are known to be critical for describing transition-metal complexes. This is evident in all the stable intermediates b, d, and e as well as the accompanying transition-state structures. With this in mind, it is very noticeable (Table 4) that the MP2 level of theory is incorrect by a large margin. MP2 suggests that e is about 110 kcal/mol more stable than a, in contrast to about 7 kcal/mol predicted by CCSD(T). In Table 5, we list the lowering of the interaction energy calculated at the HF level due to electron correlation calculated at the CCSD(T), CCSD, CCD, MP4, and MP2 level of theory. All these wave-function theories use HF as the reference configuration in their treatment of electron correlation. Thus, it is particularly convenient and appropriate to assess how these different electron-correlation treatments differ and vary among the intermediates on the potential-energy pathway. It can be seen, for example, that the lowering of the interaction energies calculated with CCD, CCSD, and MP4 follows semiquantitatively the lowering observed with CCSD(T) for all the intermediates. It is notable that, for MP2, the lowering in interaction energy for the c and e intermediates is much larger than the lowering obtained with CCSD(T) (142.1 and 172.8 kcal/mol for MP2 compared to 69.4 and 69.2 kcal/mol for CCSD(T) for c and e intermediates, respectively). Thus, the two problematic intermediates for the MP2 level of theory are c, the dihydride Ni(IV) state, and e, the diproton Ni(0) state. The large HOMO-LUMO gaps of the reference HF wave functions (Table 6) are not suggestive of the presence of low-lying excited states or near-degeneracy effects as the cause of an overestimation of the correlation energy and the failure of MP2 (i.e., breakdown of the perturbation theory ansatz of a single reference state). However, it is notable in Table 6 that the HOMO and
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12721 TABLE 7: Relative Gas-Phase Energy (∆E), Gas-Phase Free Energy [∆G(g)], and Solution (Acetonitrile) Free Energy [∆G(g+s)] for the H2 Oxidation by [Ni(PMe2NMe2)2]2+ Calculated at the B3P86/bs1 Level of Theory (in kcal/mol) a b TS(bc) c TS(cd) TS(bd) d TS(de) e
∆E
∆G(g)
∆G(g+s)
0.0 -3.6 6.1 4.9 5.4 0.9 -7.8 -2.7 -10.8
0.0 8.5 17.5 16.3 16.8 12.7 6.8 10.8 6.7
0.0 8.2 17.0 16.8 17.3 12.2 3.7 9.8 5.2
LUMO states of the diproton species e are shifted upward considerably from what they are in a, b, c, and d. Species e is characterized as a d10 state formally and also on the basis of electron-population analysis of the HF wave function. Perhaps more than for any other state, we believe that HF gives a poor description of e. The in-out electron correlation effects53,54 in this d-shell-electron-rich state are especially significant, making HF a poor reference state for perturbation treatment. MP2 overshoots electron correlation from double excitations included by perturbation treatment. Higher excitations in MP4, CCD, CCSD, and CCSD(T) correct for this MP2 behavior (Table 5). Species c (dihydride) is also electron rich at Ni (in spite of being d6 formally) and suffers from the same failure as species e, even if slightly less. It is notable that CCD theory, which includes in addition to double exictations (as MP2 does), also doubles of doubles (quadruples), behaves qualitatively similarly to CCSD(T). The comparison between MPn, CCSD(T), and DFT/B3LYP allows us to draw some additional conclusions on the discrepancy between B3P86 and PBE, on one side, and other DFT descriptions, on the other side. First, we note that BLYP often performs worse than hydrid-GGA functionals55 even though exceptions exist.56 The inclusion of a fraction of the HF exchange interactions, as done in B3LYP, improves usually the accuracy in many circumstances.57 In the current situation however, adding part of the exact exchange does not help but rather degrades the quality of the results as one goes from BLYP to B3LYP. Likewise, the BMK functional contains a sizable fraction of HF exchange (42%), the errors of which are not properly compensated for by the correlation term. In contrast, the Perdew correlation functional used in the B3P86 functional, or its improved extension in the PBE functional, seems to perform markedly better, giving energies in qualitative accord with the CCSD(T) values. Free-Energy Reaction Pathways. So far, we have calculated energies and discussed relative energies. For comparison with experimental data, it is necessary to calculate free energies. In what follows, we account for the solvation and vibrational contributions to the free energies and assess how the energy landscape is modified accordingly. Vibrational contributions in the harmonic approximation to the free energy were calculated via standard equations.58 Solvation free energies for acetonitrile, the typical solvent employed in electrocatalytic studies of these complexes, were calculated by using a dielectric continuum description of solvation (C-PCM,59,60 Bondi radii,61 acetonitrile, dielectric constant ε ) 35.688) at the B3P86/bs1 level of theory. The results are reported in Table 7 and Figure 3. Two pathways lead to H2 splitting. Both of them involve the formation of a loosely bound Ni(II)/dihydrogen adduct b with a high free energy (∆G(b) ) +8.2 kcal/mol). The formation of
12722
J. Phys. Chem. A, Vol. 114, No. 48, 2010
Chen et al. TABLE 8: Comparison of the Relative Solution (Acetonitrile) Free Energy for the H2 Oxidation by [Ni(PMe2NMe2)2]2+ and [Ni(PCy2NMe2)2]2+ Calculated at the B3P86/bs1 Level of Theory (in kcal/mol) a b TS(bc) c TS(cd) TS(bd) d TS(de) e
Figure 3. Lowest (free-) energy reaction pathway for the oxidation of H2 catalyzed by Ni(PMe2NMe2)22+ at the B3P86/bs1 level of theory. The blue line is the relative gas-phase energy (∆E) pathway, the green line is the gas-phase free-energy [∆G(g)] pathway, and the red line is the solution free-energy [∆G(g+s)] pathway.
the dihydrogen adduct induces a distortion of the nickel center resulting in a twist angle between the P-Ni-P of 57.6°, see Table S1 in the Supporting Information. The complex is characterized by a Ni-H distance of 1.780 Å and an H-H distance of 0.801 Å, with the H-H axis oriented toward the N atom of the two pendant amines. The lowest-free-energy pathway (∆G‡(bd) ) 12.2 kcal/mol) involves the formation of a Ni(II)/proton-hydride species d (∆G(d) ) 3.7 kcal/mol), which results from an asymmetric heterolytic H-H splitting as a consequence of the strong interaction between the lone pair orbital of the amine with the antibonding orbital σ*(H2) and at the same time the interaction of the bonding orbital σ(H2) with an empty Ni d orbital. The d species is characterized by a N-H(δ+) · · · (δ-)H-Ni interaction (H(δ+) · · · H(δ-) distance of about 1.646 Å). This proton-hydride species features a trigonal bipyramidal metal center with the hydride in an equatorial position (Ni-H distance of 1.5 Å). In a final but non rate-limiting step, the proton-hydride d generates a tetrahedal Ni(0)/diproton species, e, upon crossing a small activation free energy (6 kcal/mol). The proton-hydride d can be obtained also from the adduct b via a higher free-energy pathway characterized by a symmetric H2 splitting to form a cis dihydride intermediate (∆G(c) ) 16.8 kcal/mol). This intermediate lies in an extremely shallow basin on the free-energy surface and can convert back to the H2 adduct, b (∆G‡(cb) ) 0.2 kcal/mol), or easily evolve into d (∆G‡(cd) ) 0.5 kcal/mol). For the current model compound, this pathway is ∼5 kcal/mol higher in free energy than the heterolytic bond cleavage pathway, and therefore, it is expected to be only a minor reaction channel at room temperature. We note however that the relative free energies of the two routes will depend on the fine details of the structure, and in some instances, the homolytic route has been observed.28 The major effect of zero-point energy and entropy corrections is to make all of the reaction intermediates unstable with respect to the reactants, resulting in an overall endothermic reaction. Decomposition of the free energy indicates that the H2 translational entropy loss upon its coordination to the metal center is the major free-energy bottleneck, contributing as much as 11 kcal/mol to the reaction barrier. Similar effects have been observed in molecular association reactions and in particular in enzymatic activity.19,62 This energetic term would make this particular (model) compound a catalyst for H2 production but not for H2 oxidation. However, it is important to recognize that,
[Ni(PMe2NMe2)2]2+
[Ni(PCy2NMe2)2]2+
0.0 8.2 17.0 16.8 17.3 12.2 3.7 9.8 5.2
0.0 9.8 19.2 16.0 13.4 5.8 8.3 2.7
once the H2 molecule enters the coordination sphere of the metal, it could react with very small energy barriers to form a Ni(II)/ proton-hydride species, which would immediately evolve into a Ni(0)/diproton intermediate e. Solvation stabilizes the intermediates and transition states as they have positively charged protons residing on the periphery of the molecule where solvation can be most influential at stabilizing this charge. However, this relative stabilization is negligible for all of the species, except for the proton-hydride d and, to a lesser extent, the diproton species e. This can be explained in terms of the large electric dipole moment, µ, that characterizes these two species and in particular e (µ ) 6.7 D for d and µ ) 5.6 D for e), in contrast to a and b (µ ) 0.0 D and µ ) 0.1 D, respectively). Comparison with a More Realistic Model Catalyst. Given the sensitivity of the catalytic activity (preferentially H2 oxidation or H2 production) to the nature of the R and R’ substituents on P and N, it is of interest to compare the free-energy profile characterized above for the all-methyl Ni(P2MeN2Me)22+ model catalyst with that of a real catalyst for H2 oxidation. To this end, intermediates and transition states were also characterized for Ni(P2CyN2Me)22+19 at the same DFT(B3P86)/bs1 level of theory. The results are reported in Table 8. By and large, the computed reaction-energy landscape is in good agreement with that obtained for the cyclohexyl/methyl species and the methyl/methyl species, especially with respect to the fact that the rate-limiting step remains the heterolytic cleaveage of the H2 bond which has a similar energy barrier for the two species. In addition, the heterolytic splitting pathway (b f d) is found to be lower in energy than the homolytic pathway (b f c f d) by about the same amount. Note however that, for Ni(P2CyN2Me)22+, we could not locate a transition state connecting the high-energy dihydride species c to the protonhydride d. However, there are a few notable changes in the energy landscape. The H2 binding energy is slightly lower for Ni(P2CyN2Me)22+. The proton-hydride species d is appreciably higher in energy for Ni(P2CyN2Me)22+, whereas the diproton species e is significantly lower in energy, with respect to the corresponding species in Ni(P2MeN2Me)22+. The higher stability of the diproton species e in Ni(P2CyN2Me)22+ is consistent with its activity in oxidizing H2. It is worth noting that the reaction mechanism discussed here is a simplification of the real reaction mechanism, because it takes into account only one of the possible P2N2 ligand conformers of the catalysts and neglects the complex conformational changes occurring upon H2 splitting. This complexity will be analyzed in a forthcoming publication. Nevertheless, the comparison between the two complexes with different ligands is an important overall validation for the computational DFT methodology discussed throughout this paper and that we have employed to study the Ni(P2N2)22+/H2 chemistry.
Homogeneous Ni Catalysts for H2 Oxidation and Production Conclusions A comparative analysis of different quantum-chemistry protocols for the study of the early stages of hydrogen oxidation by a model system representative of a family of homogeneous Ni(II) catalysts has been presented. First, it was noted that basisset errors can be significant, and it is paramount that the description of the Ni center be handled with great care. Second, different DFT methodologies were tested against post-HF schemes such as MP2 and MP4, CCSD and, in particular, CCSD(T). The results showed that correlation effects are extremely important for this Ni complex and that a simple MP2 treatment is not able to capture these effects because of its inability to correct the correlation energy for deficiencies in the HF reference ground state for all intermediates. Rather, a more complete treatment of electron-correlation effects is required. In contrast, the hybrid functionals M06 and B3P86 and the GGA functional PBE are qualitatively similar and consistent with CCSD(T), predicting an overall exothermic reaction energy for the formation of the intermediate diproton species. Notably, the B3LYP hybrid functional, commonly employed in the study of metal complexes, yields far higher energies for all of the reaction intermediates, with an overall mildly endothermic reaction. Additional comparisons with the BLYP functional suggest that the LYP correlation functional is not able to describe the correlation effects in this type of complex. Harmonic free energies suggest that the loss of H2 translational entropy coordination to the metal center is a major contributor to the major free-energy bottleneck for H2 oxidation. Nevertheless, once the H2 molecule gets into the coordination sphere of the metal center, it quickly evolves to a Ni(II)/protonhydride species, which immediately yields a Ni(0)/diproton intermediate. Acknowledgment. We thank Dan DuBois, Jenny Yang, Ping Yang, Niri Govind, Karol Kowalski, Piotr Pecuch, and VassilikiAlexandra Glezakou for useful discussions. This research was carried out in the Center for Molecular Electrocatalysis, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences. Pacific Northwest National Laboratory is operated for the U.S. Department of Energy by Battelle under Contract No. DEAC06-76RLO 1830. Computational resources were provided at W. R. Wiley Environmental Molecular Science Laboratory (EMSL), a national scientific user facility sponsored by the Department of Energy’s Office of Biological and Environmental Research located at Pacific Northwest National Laboratory and the National Energy Research Scientific Computing Center (NERSC) at Lawrence Berkeley National Laboratory. Supporting Information Available: This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Catalysis Without Precious Metals; Bullock, R. M., Ed.; WileyVCH: Weinheim, 2010. (2) Fontecilla-Camps, J. C.; Volbeda, A.; Cavazza, C.; Nicolet, Y. Chem. ReV. 2007, 107, 4273. (3) Garcin, E.; Vernede, X.; Hatchikian, E. C.; Volbeda, A.; Frey, M.; Fontecilla-Camps, J. C. Structure 1999, 7, 557. (4) Higuchi, Y.; Ogata, H.; Miki, K.; Yasuoka, N.; Yagi, T. Structure 1999, 7, 549. (5) Nicolet, Y.; de Lacey, A. L.; Verned`e, X. J. Am. Chem. Soc. 2001, 123, 1596. (6) Pereira, A. S.; Tavares, P.; Moura, I.; Moura, J. J. G.; Huynh, B. H. J. Am. Chem. Soc. 2001, 123, 2771. (7) Peters, J. W. Curr. Opin. Struct. Biol. 1999, 9, 670.
J. Phys. Chem. A, Vol. 114, No. 48, 2010 12723 (8) Peters, J. W.; Lanzilotta, W. N.; Lemon, B. J.; Seefeldt, L. C. Science 1998, 282, 1853. (9) Volbeda, A.; Fontecilla-Camps, J. Dalton Trans. 2003, 4030. (10) Volbeda, A.; Garcin, E.; Piras, C.; de Lacey, A. L.; Fernandez, V. M.; Hatchikian, E. C.; Frey, M.; Fontecilla-Camps, J. C. J. Am. Chem. Soc. 1996, 12989. (11) Rakowski DuBois, M.; DuBois, D. L. Acc. Chem. Res. 2009, 42, 1974. (12) Rakowski DuBois, M.; DuBois, D. L. Chem. Soc. ReV. 2009, 38, 62. (13) Qi, X.-J.; Fu, Y.; Liu, l.; Guo, Q.-X. Organometallics 2007, 26, 4546. (14) Kova´cs, G.; Pa´pai, I. Organometallics 2006, 25, 820. (15) Qi, X.-J.; Liu, l.; Fu, Y.; Guo, Q.-X. Organometallics 2006, 25, 5879. (16) Nimlos, M. R.; Chang, C. H.; Curtis, C. J.; Miedaner, A.; Pilath, H. M.; DuBois, D. L. Organometallics 2008, 27, 2715. (17) Wilson, A. D.; Newell, R. H.; McNevin, M. J.; Muckerman, J. T.; Rakowski DuBois, M.; DuBois, D. L. J. Am. Chem. Soc. 2006, 128, 358. (18) Wilson, A. D.; Shoemaker, R. K.; Miedaner, A.; Muckerman, J. T.; DuBois, D. L.; DuBois, M. R. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 6951. (19) (a) Yang, J. Y.; Chen, S.; Dougherty, W. G.; Kassel, W. S.; Bullock, R. M.; DuBois, D. L.; Raugei, S.; Rousseau, R.; Dupuis, M.; Rakowski DuBois, M. Chem. Comm. 2010, DOI: 10.1039/c0cc03246h. (b) Yang, J. Y. et al. Manuscript regarding the work of Ni(P2CyN2Me)22+ is in preparation. (20) Kubas, G. J. Dihydrogen and σ-Bond Complexes: Structure, Theory, and ReactiVity; Kluwer Academic/Plenum Publishers: New York, 2001. (21) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (22) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 78, 1396. (23) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (24) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (25) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B 1988, 37, 785. (26) Miehlich, B.; Savin, A.; Stoll, H.; Preuss, H. Chem. Phys. Lett. 1989, 157, 200. (27) Niu, S.; Hall, M. B. Chem. ReV. 2000, 100, 353. (28) Yang, J. Y.; Bullock, R. M.; Shaw, W. J.; Twamley, B.; Fraze, K.; DuBois, M. R.; DuBois, D. L. J. Am. Chem. Soc. 2009, 131, 5935. (29) Balabanov, N. B.; Peterson, K. A. J. Chem. Phys. 2005, 123, 064107. (30) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (31) Woon, D. E.; Dunning, T. H., Jr. J. Chem. Phys. 1993, 98, 1358. (32) Feller, D. J. Chem. Phys. 1992, 96, 6104. (33) Rassolov, V. A.; Pople, J. A.; Ratner, M. A.; Windus, T. L. J. Chem. Phys. 1998, 109, 1223. (34) Hay, P. J. J. Chem. Phys. 1977, 66, 4377. (35) Wachters, A. J. H. J. Chem. Phys. 1970, 52, 1033. (36) Andrae, D.; Ha¨uβermann, U.; Dolg, M.; Stoll, H.; Preuβ, H. Theor. Chem. Acc. 1990, 77, 123. (37) Schafer, A.; Horn, H.; Ahlrichs, R. J. Chem. Phys. 1992, 97, 2571. (38) Bauschlicher, C. W., Jr.; Langhoff, S. R.; Barnes, L. A. J. Chem. Phys. 1989, 91, 2399. (39) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (40) Boese, A. D.; Martin, J. M. L. J. Chem. Phys. 2004, 121, 3405. (41) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Phys. 2005, 123, 194101. (42) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215. (43) Straatsma, T. P.; Apra`, E.; Windus, T. L.; Bylaska, E. J.; Jong, W. d.; Hirata, S.; Valiev, M.; Hacklr, M.; Pollack, L.; Harrison, R.; Dupuis, M.; Smith, D. M. A.; Nieplocha, J.; Tipparaju, V.; Krishnan, M.; Auer, A. A.; Brown, E.; Cisneros, G.; Fann, G.; Fru¨chtl, H.; Garza, J.; Hirao, K.; Kendall, R.; Nichols, J.; Tsemekhman, K.; Wolinski, K.; Anchell, J.; Bernholdt, D.; Borowski, P.; Clark, T.; Clerc, D.; Dachsel, H.; Deegan, M.; Dyall, K.; Elwood, D.; Glendening, E.; Gutowski, M.; Hess, A.; Jaffe, J.; Johnson, B.; Ju, J.; Kobayashi, R.; Kutteh, R.; Lin, Z.; Littlefield, R.; Long, X.; Meng, B.; Nakajima, T.; Niu, S.; Rosing, M.; Sandrone, G.; Stave, M.; Taylor, H.; Thomas, G.; Lenthe, J. v.; Wong, A.; Zhang, Z.; Pacific Northwest National Laboratory: Richland, Washington. (44) Kendall, R. A.; Apra`, E.; Bernholdt, D. E.; Bylaska, E. J.; Dupuis, M.; Fann, G. I.; Harrison, R. J.; Ju, J.; Nichols, J. A.; Nieplocha, J.; Straatsma, T. P.; Windus, T. L.; Wong, A. T. Comput. Phys. Commun. 2000, 128, 260. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.;
12724
J. Phys. Chem. A, Vol. 114, No. 48, 2010
Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian, Inc.: Wallingford CT, 2009. (46) Feller, D. J. Comput. Chem. 1996, 17, 1571. (47) Schuchardt, K. L.; Didier, B. T.; Elsethagen, T.; Sun, L.; Gurumoorthi, V.; Chase, J.; Li, J.; Windus, T. L. J. Chem. Inf. Model. 2007, 47, 1045. (48) Hajgato´, B.; Szieberth, D.; Geerlings, P.; Proft, F. D.; Deleuze1, M. S. J. Chem. Phys. 2009, 131, 224321. (49) Tekarli, S. M.; Drummond, M. L.; Williams, T. G.; Cundari, T. R.; Wilson, A. K. J. Phys. Chem. A 2009, 113, 8607. (50) Bu¨hl, M.; Reimann, C.; Pantazis, D. A.; Bredow, T.; Neese, F. J. Chem. Theory Comput. 2008, 4, 1449. (51) Rappoport, D.; Crawford, N. R. M.; Furche, F.; Burke, K. In Computational Inorganic and Bioinorganic Chemistry; Solomon, E. I., King, R. B., Scott, R. A., Eds.; Wiley: Chichester, 2008.
Chen et al. (52) Zhao, Y.; Truhlar, D. G. J. Chem. Phys. 2006, 124, 224105. (53) Davidson, E. R.; Jones, L. L. J. Chem. Phys. 1962, 37, 1918. (54) Goddard, W. A., III.; Harding, L. B. Annu. ReV. Phys. Chem. 1978, 29, 363. (55) Sousa, S. F.; Fernandes, P. A.; Ramos, M. J. J. Phys. Chem. A 2007, 111, 10439. (56) Furche, F.; Perdew, J. P. J. Chem. Phys. 2006, 124, 044103. (57) Siegbahn, P. E. M.; Blomberg, M. R. A. Chem. ReV. 2000, 100, 421. (58) McQuarrie, D. A. Statistical Mechanics; Happer & Row: New York, Evanston, San Francisco, London, 1976. (59) Barone, V.; Cossi, M. J. Phys. Chem. A 1998, 102, 1995. (60) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. J. Comput. Chem. 2003, 24, 669. (61) Bondi, A. J. Phys. Chem. 1964, 68, 441. (62) Siebert, X.; Amzel, L. M. Proteins 2003, 54, 104.
JP106800N
