A Multiconfigurational Perturbation Theory and Density Functional

Feb 15, 2012 - Of the tested functionals, B3LYP performs best in reproducing the binding energy, while the PBE0 functional gives the best structures. ...
3 downloads 11 Views 912KB Size
Article pubs.acs.org/JCTC

A Multiconfigurational Perturbation Theory and Density Functional Theory Study on the Heterolytic Dissociation Enthalpy of First-Row Metallocenes Quan Manh Phung, Steven Vancoillie,* and Kristine Pierloot Department of Chemistry, University of Leuven, Celestijnenlaan 200F, B-3001 Heverlee, Belgium ABSTRACT: The heterolytic dissociation enthalpy of a series of first-row metallocenes M(C5H5)2, M = V, Mn, Fe, and Ni, was studied by (restricted) multiconfigurational perturbation theory and density functional theory. The results were compared directly to the experimental values, taking into account all necessary contributions to the relative energy. Of the tested functionals, B3LYP performs best in reproducing the binding energy, while the PBE0 functional gives the best structures. High quality multiconfigurational perturbation calculations were also carried out, demonstrating the superior performance of a larger, restricted active space. The spin crossover behavior of manganocene is correctly predicted by multiconfigurational perturbation theory as opposed to the three functionals B3LYP, PBE0, and M06, which (severely) overstabilize the high-spin with respect to the low-spin state.

1. INTRODUCTION After the discovery of ferrocene (Fe(Cp)2), now 60 years ago,1,2 a whole class of other metallocenes appeared,3−6 covering almost the whole first-row transition metal series and some second- and third-row metals. Starting from 1990, metallocenes became the subject of numerous theoretical studies. The first studies only aimed at investigating the equilibrium geometry of metallocenes. In 1991, Park and Almlöf7 predicted a value of 1.65−1.67 Å for the iron− cyclopentadienyl distance in ferrocene, using the modified coupled-pair functional method, as compared to the experimental Fe−Cp distance of 1.661 Å.8 In the same year, the first density functional theory (DFT) study, using both local and nonlocal functionals, was carried out to investigate the geometry of ferrocene,9 and in 1997, Matsuzawa et al.10 extended this study to the geometry of a variety of metallocenes. Mayor-Lopez and Weber11 systematically studied four metallocene structures M(C5H5)2 (M = V, Mn, Fe, and Ni) with the BPW91 functional. Xu et al.12 performed a detailed study of the electronic and molecular structures of the first-row transition metal series metallocenes using the B3LYP, BLYP, and BP86 functionals. In most DFT calculations, the deviation of the bond lengths and metal−carbon distance with respect to the experimental structure is smaller than 0.02 Å.13 In 1995, Pierloot et al.14 found an iron−cyclopentadienyl distance of 1.643 Å, using the complete active space selfconsistent field (CASSCF) method followed by second-order perturbation theory (CASPT2). So far, the most accurate equilibrium geometry of ferrocene was obtained by Coriani et al.13 Their computed structure at the CCSD(T) level of theory using a TZV2P+f basis set gave an iron−cyclopentadienyl distance of 1.655 Å, in excellent agreement with the experimental value of 1.661 Å.15 Despite the fact that metallocenes had been discovered a long time ago, their binding energy was not measured until the experimental study of Ryan et al.,16 who indirectly measured both the heterolytic and homolytic dissociation enthalpy of © 2012 American Chemical Society

ferrocene, vanadocene, manganocene, and nickelocene. Since then, a number of theoretical studies on binding energies of metallocenes were carried out, for both the heterolytic and homolytic dissociations. In this work, we will focus on the heterolytic dissociation enthalpy. Even until now, there have been large discrepancies between experimental and theoretical results for this quantity. The dissociation enthalpy of some metallocenes is difficult to compute accurately, as it involves a change in spin between the (usually) low-spin metallocene and the high-spin metal ion. This requires a correct treatment of the relative energies of different spin-states, which is particularly problematic for many DFT functionals17−40 but can be achieved by means of ab initio wave-function-based methods. Of these, multiconfigurational perturbation theory is the only reasonably cheap option for studying relatively large transition metal complexes. Using the computationally cheap DFT method, a systematic investigation of the binding energy of a series of metallocenes can easily be performed. Results from the study of MayorLopez and Weber11 showed that both LDA and BPW91 functionals have trouble in predicting the binding energy of metallocenes, overestimating its value by up to +98 and +28 kcal/mol, respectively. Swart41 examined the heterolytic dissociation enthalpy of metallocenes of the first-row transition metals (Sc−Zn) extended with alkaline-earth metals (Mg, Ca) and several second-row transition metals (Ru, Pd, Ag, Cd). The latter study was performed using the OPBE functional, which was found earlier to perform well for problems involving the relative energies of states with different spin.29 The trends in metal−ligand bonding in these complexes were analyzed using an energy decomposition analysis. Qualitative results of the binding energy were obtained, with errors of about +26.0, +4.8, +23.4, and +48.5 kcal/mol for the metallocenes of V, Mn, Fe, and Ni, respectively. The overestimation of the binding energy Received: December 6, 2011 Published: February 15, 2012 883

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

give reliable results for transition metal compounds;59 PBE,54,55,60,61 which is widely used for solid-state calculations; PBE0;54,55,60−62 B3LYP;54−57,63,64 and TPSSh.54,55,60,65,66 The latter three are hybrid functionals with different amounts of exact HF exchange of 25%, 20%, and 10%, respectively; B97D,67 a GGA functional which is parametrized for describing dispersion interaction; and M06, a highly parametrized hybrid meta-GGA containing 27% HF exchange, which is proposed to be appropriate for main group elements, organometallics, kinetics, and noncovalent bonds.68 After testing on ferrocene, only basis set II was used for the remaining metallocenes, and the functionals were restricted to those that performed well for the dissociation enthalpy. These are B3LYP, PBE0, and M06. To obtain the heterolytic dissociation enthalpy of the metallocenes at room temperature, the electronic binding energy ΔE (EM2+ + 2 × ECp− − EM(Cp)2) was first determined. Additionally, the zero-point energy (EZPE) and thermal (Ethermal) corrections were added, which were obtained from frequency calculations on M(Cp)2. Furthermore, a counterpoise correction69,70 was added (ECPC) to account for the basisset superposition error (BSSE). When dealing with first-row transition metals, one often neglects relativistic effects since they are assumed to provide a minor contribution. However, as can be seen in the work of Klopper and Lüthi,46 the contribution of (scalar) relativistic effects to the binding energy of ferrocene is as large as +6 kcal/mol. Therefore, in the DFT part of our study, a relativistic correction (Erel) was also taken into account by using a second-order Douglas−Kroll−Hess (DKH) Hamiltonian.71−73 It is well-known that DFT has difficulties with properly describing intermolecular interactions, in particular van der Waals forces,74 with many functionals predicting an unbound state for, e.g., rare-gas dimers75 or π−π stacking.76 Therefore, an empirical dispersion correction Edisp (DFT-D2)67 was taken into account to correctly describe the dispersion interaction between two cyclopentadienyl rings (Cp−). This DFT-D2 correction was thoroughly tested and proven to be reliable and “robust” in describing both inter- and intramolecular interactions in hundreds of different systems. For example, using the BLYP functional with dispersion correction (BLYP-D), the binding energy of the benzene dimer is predicted to be within 0.2 kcal/mol of the values obtained with the computationally more expensive CCSD(T)/ CBS method.77 The corresponding errors of M06-2X, B3LYPD, and PBE-D are 0.7, 0.5, and 0.15 kcal/mol, respectively. Hence, using the DFT-D2 correction, we expect to obtain the same accuracy when treating the cyclopentadienyl dimer system. Since the DFT-D2 parameters are only available for the B3LYP, PBE, BP86, and B97D functionals, dispersion corrections for the PBE0 and TPSSh functionals were taken from B3LYP-D. Because the M06 functional was claimed to give good results for the dispersion interaction,68,78 no dispersion correction was taken into account for this functional. All CASSCF/CASPT2 and RASSCF/RASPT2 calculations were performed with the MOLCAS 7.6 package.79 Extended ANO-RCC basis sets (with a total of 1367 basis functions) were used with the following contractions: [10s9p8d6f4g2h] for the metal atom,80 [8s7p4d3f1g] for carbon,81 and [6s4p3d1f] for hydrogen.82 This basis set was previously shown to perform well for the computation of the ferrocene binding energy, giving rise to a small BSSE.47 The Cholesky decomposition technique was used to approximate the twoelectron integrals, using a threshold of 10−6 au.83 Scalar relativistic effects were included using the standard second-

with GGA functionals can be reduced by using hybrid functionals, which contain a small admixture of Hartree−Fock (HF) exchange. Frunzke et al.42 used the B3LYP functional and obtained a binding energy that was 10 kcal/mol higher than the experimental value of the heterolytic dissociation enthalpy of ferrocene. On the other hand, employing the same functional and 6-311+G* basis set, Padma Malar43 found a value of 20 kcal/mol below the experimental one. Furche and Perdew44 benchmarked various functionals for the homolytic dissociation enthalpy of metallocenes (V, Mn, Fe, and Ni). These authors showed that the binding energies decrease in the order LSDA ≫ PBE > TPSS > TPSSh ≈ BP86 > experimental > B3LYP. However, none of these functionals yielded binding energies that are correct to within the experimental uncertainty. To the best of our knowledge, previous ab initio studies with wave function methods have only considered ferrocene. The first study was performed by Pierloot et al.14 using the CASSCF/CASPT2 approach.45 The computed heterolytic dissociation enthalpy closely agreed with experimental values. However, the large basis set superposition error (+30 kcal/mol) and the lack of zero-point energy corrections gave a relatively large uncertainty to the result. The main problem of these calculations was that the size of the basis sets had to be limited because of computational limitations. In a later study,46 it was shown that, when extrapolated to the basis set limit both CASPT2 and CCSD(T) predict a heterolytic dissociation enthalpy of 653 kcal/mol, much higher than the experimental value of 635 ± 6 kcal/mol. Recently,47 we have shown that with large basis sets and an extended active space, restricted multiconfigurational perturbation theory (RASPT2)48 yields lower values for the dissociation enthalpy, much closer to the experimental value than the CASPT2 result with the same basis set. In this study, we make use of the complete and restricted active space methods followed by the second-order perturbation theory (CASSCF/CASPT245 and RASSCF/RASPT248) to compute the heterolytic binding energies of a number of metallocenes. Considering the failure of many functionals to predict the correct heterolytic dissociation enthalpy, we also wanted to investigate the performance of a number of important functionals, including all necessary corrections that might be important, in particular (counterpoise) corrections for basis set superpostion errors, scalar relativistic effects, and corrections for the lack of dispersion interaction in most DFT functionals. We have limited our study to those first-row metallocenes for which an experimental value of the heterolytic dissociation enthalpy is available.16 These are ferrocene (Fe(Cp)2), vanadocene (V(Cp)2), manganocene (Mn(Cp)2), and nickelocene (Ni(Cp)2).

2. COMPUTATIONAL DETAILS All DFT calculations were performed with the programs Gaussian 0949 (M06 functional) and TURBOMOLE v6.0350 (all other functionals). To construct a general strategy and to test the accuracy of DFT, we first studied the binding energy of ferrocene with different basis sets and functionals. All DFT calculations were performed using the unrestricted Kohn− Sham approach, but no spin contamination was found in all cases. Three types of basis sets were used: 6-31G* for C and H atoms and def2-TZVP for the metal atom51 (basis set I); def2TZVP for C and H atoms52 and def2-QZVPP for the metal atom53 (basis set II); and def2-QZVPP for all atoms53 (basis set III). The functionals used were BP86,54−58 which is believed to 884

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

order Douglas−Kroll−Hess Hamiltonian.71−73 Since it is computationally unfeasible to optimize the metallocene structure with CASPT2 or RASPT2, the ground state geometry of the metallocenes and cyclopentadienyl was taken from the DFT structure that gave the lowest CASPT2/RASPT2 energy, i.e., the PBE0 structure. In the perturbation step, an imaginary level shift84 of 0.1 au was used to prevent intruder states, while the default IPEA shift for the zeroth-order Hamiltonian85 (0.25 au) was used. All valence electrons, including the metal (3s,3p) semicore electrons, were included in the CASPT2/RASPT2 calculations. Figure 1 shows the distribution of the electrons over the (predominant) metal 3d orbitals in the ground state of the

not strongly interact. Metal−Cp bonding occurs through charge donation from the Cp e1″(π) into the metal e1″ (dxz,dyz) orbitals, counteracted by backdonation from the metal e2′(dxy,dx2 − y2) into the Cp e2′(π*) orbitals. The Cp e1″ (π) and e2′ (π*) cyclopentadienyl orbitals should therefore be added to the metal d orbitals in the active space, yielding a total of nine active orbitals. On top of these nine orbitals, extra d′ orbitals should be included in the active space for the d5 Mn(Cp)2, d6 Fe(Cp)2, and d8 Ni(Cp)2 complexes. However, the virtual e2′ shell already present there is not pure Cp π* but rather contains a mixture of Cp π* and metal (3dxy ′ , 3dx′ 2−y2 ) character and may therefore (in an economical calculation) serve to account for both types of correlation. Furthermore, no correlating 3d′ orbital is strictly necessary if the corresponding 3d orbital is empty, as is the case for the e1″(dxz,dyz) shell in ferrocene. This leaves only the dz2′ orbital to be added, resulting in an active space of 10 orbitals including 10 electrons. On the other hand, in the 6A1′ state of both manganocene and nickelocene, the dxz and dyz orbitals are not empty, and their corresponding d′ orbitals were included, in the case of manganocene both for the 6 A1′ and 2E2′ states, giving rise to 12 active orbitals occupied by nine (Mn(Cp)2) and 12 (Ni(Cp)2) electrons, respectively. Vanadocene was given a somewhat special treatment at the CASSCF level. Since there is no double-shell effect here, the “basic” active space needed is limited to only nine orbitals, thus leaving room for including extra cyclopentadienyl orbitals. Both the e1′ (π) and e2″ (π*) were included, leading to an active space of 13 orbitals containing 11 electrons. In a previous study,47 we have shown that a more accurate description of the dissociation enthalpy of ferrocene may be obtained from a RASPT2 treatment, based on an extended active space of 18 orbitals, obtained by extending the (10,10) active space of this molecule with an extra virtual e2′(d′) and e1″(d′) shell to improve the description of the 3d double-shell effect, as well as with the cyclopentadienyl e1′ (π) and e2″ (π*) orbitals. With a total of 18 active orbitals, CASSCF becomes computationally unaffordable, and one has to turn to RASSCF instead. In the latter method, the active space is further subdivided into three subspaces RAS1, RAS2, and RAS3. All possible occupancies within RAS2 are still allowed, but the number of holes in the initially doubly occupied RAS1 orbitals and the number of electrons in the initially empty RAS3 orbitals is restricted to two. The bonding−antibonding combinations of the metal (3dxy, 3dx2−y2) and Cp π* orbitals in representation e2′ and of the Cp π and metal (3dxz, 3dyz) orbitals in representation e1″ were put in RAS2, which should furthermore include any singly occupied orbital.47 This then gives an eight orbital RAS2 in ferrocene and nickelocene. In the 6 A1′ state of manganocene, the a1′(dz2) orbital becomes singly occupied. This orbital was therefore included in RAS2 for both calculated states 6A1′ and 2E2′ of this molecule, increasing the number of RAS2 orbitals to nine. Of the remaining orbitals, the occupied orbitals go in RAS1 and the empty orbitals in RAS3. The subdivision of the 18-orbital active space of ferrocene is presented in Figure 2. For vanadocene, a different RASPT2 procedure was used. Since there cannot be a significant double-shell effect in this d3 molecule and since all important Cp π orbitals were already included in the CASSCF/CASPT2 calculations, the results of the CASPT2 calculations on vanadocene should be of comparable quality to the RASPT2 results for the other three molecules. However, as we shall see further (e.g., in Figure 3),

Figure 1. Electron configuration of metallocenes.

different metallocenes. The d5 configuration of the metal ion in manganocene may give rise to a low-spin 2E′2 or high-spin 6A′1 ground state, depending on the temperature. Therefore, both states were considered. In a traditional CASPT2 calculation, keeping the global size of the active space limited to at most 14−16 orbitals is an important consideration. For a 3dn metallocene complex, an appropriate complete active space should then at least consist of (a) the metal 3d orbitals and those Cp− π orbitals that are involved in covalent interactions with these d orbitals and (b) in case the 3d shell is more than half-filled, an additional shell of d orbitals, denoted as 3d′, to account for the so-called 3d double-shell effect. The 3d doubleshell effect is related to the presence of a large number of electrons in a compact 3d shell, giving rise to important radial correlation effects. These correlation effects cannot be adequately described by means of second-order perturbation theory and should therefore be accounted for in the zerothorder wave function of a CASPT2 or RASPT2 calculation by including an additional, more diffuse, d-shell (called either 3d′ or 4d) in the active space. The effect was first observed in a series of CASPT2 calculations on the Ni atom86 and has since then been found to systematically affect the results obtained from CASPT2 for first-row transition metal complexes. For a more extensive discussion of the double-shell effect, see refs 87−89. Within D5h symmetry, the metal d orbitals belong to the irreducible representations (irreps) a1′ (dz2), e2′ (dxy,dx2 − y2), and e1″ (dxz,dyz), whereas the cyclopentadienyl carbon pz orbitals form symmetry-adapted combinations giving rise to six occupied π orbitals belonging to the a1′, a2″, e1′, and e1″ irreps and four empty π* orbitals belonging to the e2′ and e2″ irreps. Only those Cp π and π* within the same irreps as the metal d orbitals can interact to form covalent metal−ligand combinations, i.e., the a1′, e2′, and e1″ irreps. Since the metal dz2 and Cp− π orbitals in irrep a1′ are energetically well separated, they do 885

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

effect of including the (3s,3p) orbitals into RAS1 was found to be negligibly small ( BP86 > TPSSh > B97D > PBE0 > M06 ≈ exptl. > B3LYP. Comparing the calculated results for the heterolytic dissociation enthalpy with the experimental value, we find that the BP86, PBE, TPSSh, and B97D functionals give poor results, with deviations larger than 50 kcal/mol for BP86 and PBE, while the B3LYP, PBE0, and M06 results closely agree with experimental results, with deviations within the uncertainty of the experimental result. The strong overestimation of the binding energy with BP86 and PBE may be, at least partially,

Table 1. Bond Lengths (in Å) and the Heterolytic Dissociation Enthalpy (in kcal/mol) of Ferrocene with the B3LYP Functional Fe−Cp (Å) C−C (Å) ΔE ΔEdisp ΔECPC ΔEZPE ΔErel ΔEthermal ΔH298° errora

basis set I

basis set II

basis set III

1.684 1.428 654.9 12.0 −13.7 −4.9 5.6 1.6 655.5 20.5

1.691 1.423 624.7 11.8 −5.6 −8.6 5.5 1.8 629.4 −5.6

1.691 1.422 619.8 11.9 −2.2 −7.9 5.7 2.2 629.5 −5.5

As compared to the experimental value of ΔH°298 = 635 ± 6 kcal/ mol.16

a

metallocenes. The more extensive basis sets II and III reduce the BSSE to about 5 and 2 kcal/mol, respectively. After including ΔECPC, the binding energy obtained with both basis sets is the same. This indicates that the rather small BSSE with both basis sets may be adequately treated with the counterpoise method, and that basis set II is suitable to obtain quantitative results for the binding energy. Basis II will therefore be used to investigate the binding energy of the other metallocenes. We note from Table 1 that both the dispersion (ΔEdisp) and relativistic corrections (ΔErel) to the B3LYP binding energy are important: together they increase the dissociation enthalpy by around 17 kcal/mol! The calculated structures of ferrocene with basis set II and different functionals are given in Table 2. The best Table 2. Experimental and Calculated Structures of Ferrocenea bond lengths (Å) Fe−Cp BP86 BP86-D PBE PBE-D B3LYP B3LYP-D PBE0 TPSSh B97D M06 CCSD(T)c exptl.c

1.649 1.630 1.642 1.628 1.691 1.669 1.652 1.643 1.638 1.652 1.655 1.661

C−C

angle (deg) C−H

1.434

1.087

1.434

1.087

1.423

1.078

1.420 1.426 1.432 1.419 1.433 1.435(2)

1.080 1.080 1.084 1.080 1.077 1.080(7)

Cp−Hb 0.94 1.93 1.08 1.79 1.00 1.92 1.48 1.27 1.89 1.72 1.03 1.7(2)

a

Bond lengths in Ångstroms, angles in degrees. bAngle between hydrogens and cyclopentadienyl ring, positive values means hydrogen is bent towards metal. cFrom ref 13.

experimental values in the literature were obtained from vibrationally corrected gas electron diffraction (GED) and neutron diffraction (ND) experiments,15,91,92 while the best computational results for the geometry were computed at the CCSD(T) level.13 For all functionals tested, the computed bond lengths are in good agreement with both the CCSD(T) and experimental results, within 0.05 Å. We first focus on the performance of different functionals in predicting the Fe−Cp distance. The best result, only 0.009 Å 887

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

Table 3. Calculated Results for the Heterolytic Dissociation Enthalpy (in kcal/mol) of Ferrocene ΔE ΔEdisp ΔECPC ΔEZPE ΔErel ΔEthermal ΔH298° errore

BP86

PBE

B3LYP

PBE0

TPSSh

B97D

M06

CASPT2a

RASPT2a

678.5 12.7 −5.4 −8.3 5.8 2.1 685.4 50.4

686.4 9.1 −6.2 −8.3 6.0 2.0 689.0 54.0

624.5 11.8 −5.6 −8.6 5.5 1.8 629.4 −5.6

633.1 11.8b −5.0 −8.7 6.3 1.9 639.4 4.4

658.5 11.8b −4.9 −8.6 6.2 1.9 664.9 29.9

661.8

641.1 c −5.4 −8.3 6.4 1.9 635.7 0.7

664.5

649.7

−3.1 −8.7

−3.8 −8.7

1.9 654.6 19.6

1.9 639.1 4.1

−5.4 −8.3d 6.0 ≃2.0d 655.4 21.1

a

CASPT2 and RASPT2 calculations are computed at the PBE0 structure. bDispersion correction estimated from the B3LYP result. cM06 was claimed to give good results of dispersion interaction.68,78 dTaken from PBE (analytical frequencies are not implemented in Turbomole for the B97D functional). eAs compared to the experimental value16 ΔH298° = 635 ± 6 kcal/mol.

and a ΔS° value of 20 ± 1 cal/mol K and that the spin crossover temperature of manganocene is 212 K. Since the value of ΔH° is small, a mixture of both high spin and low spin structures is found at room temperature. Before further discussing the spin state energetics in manganocene, we first consider the results obtained for the structures of the different metallocenes. Table 4 shows the M−

traced back to the fact that GGA functionals systematically tend to favor low-spin states (such as the 1A1′ ground state in ferrocene) over high-spin states (the 5D ground state in Fe2+).17−40 Hybrid functionals perform better in this respect because of their admixture of Hartree−Fock exchange, giving rise to a relative stabilization of the high-spin state on the metal ion and therefore to a decrease of the dissociation energy. An exception is the B97D functional. This functional performs much better than the other two GGA functionals considered here, and even better than the hybrid TPSSh functional. Turning to the results from multiconfigurational perturbation theory, we find that CASPT2 overestimates the heterolytic dissociation energy by about 20 kcal/mol, while the results are drastically improved in RASPT2, after extending the active space so as to describe the full metal 3d double-shell effect and the frontier cyclopentadienyl π−π* correlation. The CASPT2/ RASPT2 results presented here closely agree with the results obtained for ferrocene in a recent benchmark study of the RASPT2 method for (first-row) transition metal complexes,47 where these results are discussed in more detail and RASPT2 results based on alternative RAS space are presented. We note, however, that the dissociation enthalpy reported here is slightly higher, by 5−6 kcal/mol, than the RASPT2 results obtained there with the same active space and similar basis sets (basis III−IV in ref 47). The difference is due to the different structure used for the ferrocene molecule in both studies, the PBE0 structure used here giving a RASPT2 energy that is lower by 5.5 kcal/mol than the experimental structure used by Vancoillie et al.47 Nevertheless, with both structures, an excellent result is obtained from RASPT2 for the heterolytic dissocation enthalpy of ferrocene, accurate to within the experimental uncertainty. 3.2. Vanadocene, Manganocene, Nickelocene. The heterolytic dissociation enthalpy of the other metallocenes (M(C5H5)2 with M = V, Mn, Ni) was investigated at the DFT level, using the same procedure as for ferrocene but employing only the B3LYP, PBE0, and M06 functionals and basis set II. The ground state of vanadocene is a 4A2′ state, while nickelocene has a 3A2′ ground state. Manganocene is a spincrossover complex. At high temperatures, the ground state is a high-spin 6A1′ state (μ = 5.5 μB at 373 K), whereas at lower temperatures or in the solid state, manganocene forms chainlike structures and has a low-spin 2E2′ ground state (μ = 1.99 μB at 193 K).96 Within the latter state, D5h manganocene spontaneously distorts to C2v symmetry to remove the 2E2′ degeneracy. A deuterium NMR spectroscopy measurement of Köhler and Schlesinger97 showed that the equilibrium 2E2′ ⇌ 6 A1′ is characterized by a ΔH° value of 2.89 ± 0.19 kcal/mol

Table 4. Bond Lengths (in Å) of the Metallocenes with Different Functionals V(Cp)2 4

A′2

B3LYP B3LYP-D PBE0 M06 exptl.a

1.972 1.957 1.934 1.906 1.928

B3LYP B3LYP-D PBE0 M06 exptl.a

1.418 1.419 1.415 1.414 1.434

Mn(Cp)2 6

A′1

M−Cp (Å) 2.093 2.081 2.061 2.024 2.050 C−C (Å) 1.419 1.419 1.415 1.414 1.429

Fe(Cp)2 1

A′1

Ni(Cp)2 3

A′2

1.691 1.669 1.652 1.652 1.661

1.877 1.856 1.837 1.821 1.829

1.423 1.424 1.420 1.419 1.435

1.419 1.420 1.416 1.415 1.439

a

Obtained with gas electron diffraction (GED): refs 99 (vanadocene), 98 (manganocene), 15 and 100 (ferrocene), and 101 (nickelocene).

Cp and C−C optimized distances obtained with the different functionals B3LYP, B3LYP-D, PBE0, and M06, and compared to the experimental data available. For manganocene, the structure considered is for the high-spin 6A1′ state, for which an experimental gas phase structure is available.98 It can be seen that the reliability of the results depends on the complex. Even more than for ferrocene, the M−Cp distance is overestimated by B3LYP for the other metallocenes, while with M06 the same distance is underestimated more strongly than in ferrocene. Including dispersion with B3LYP-D shortens the M−Cp distance, although it is still too long with respect to experimental results, and longer than the distance obtained with PBE0 or M06. For PBE0, the error varies: while the M− Cp distance is too short for ferrocene, it is too large for highspin manganocene, vanadocene, and nickelocene. Of the tested functionals, PBE0 is the best in describing the M−Cp distance, with the largest error for manganocene (−0.011 Å). All functionals also have difficulties in describing the C−C distances in all metallocenes. All C−C distances are too short, as they were in ferrocene. 888

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

kcal/mol, whereas at the RASPT2 level, this error is slightly further reduced to 1.35 kcal/mol. The dispersion correction is important and accounts for about 3 kcal/mol of relative doublet stabilization. The significant effect of the dispersion correction is related to the shorter Mn−Cp distance in the doublet structure compared to the sextet structure, because the importance of the attractive force between the cyclopentadienyl rings increases as the rings get closer to each other. We next turn to the discussion of the dissociation enthalpy. The results obtained from DFT are presented in Tables 6−8,

We now take a closer look at the calculated spin state energetics in manganocene. In Table 5, we have collected the Table 5. Relative energy difference between the 2E2′ and 6A1′ states of manganocene ΔEHS−LS

B3LYP-D

PBE0-Da

M06

CASPT2

RASPT2

exptl.b

−3.70

−7.05

−17.20

5.77

4.93

3.58

a

PBE0 with dispersion correction as the difference between B3LYP and B3LYP-D. bCalculated (gas-phase) −ΔGtherm(T c) at the experimental (solution) crossover temperature Tc = 212K

doublet-sextet gaps computed with all methods considered. In order to confront these energy differences with experimental results, it is useful to determine an “experimental” estimation of the electronic doublet−sextet gap, starting from the experimental spin crossover temperature of 212 K, measured by Köhler and Schlesinger.97 At this temperature, the Gibbs free energy ΔG = 0. ΔG(T) is equal to ΔEHS−LS + ΔGtherm(T), with ΔEHS−LS corresponding to the electronic energy difference E( 6A1′) − E(2E2′), and ΔGtherm(T) being the thermal contributions to the free energy. The latter may be obtained from a DFT calculation of vibrational frequencies. The PBE0 functional was used for this purpose (but this quantity is in fact quite insensitive to the choice of functional). At 212 K, a value of ΔGtherm(212) of −3.58 kcal/mol was obtained. This value may therefore serve as the “experimental” electronic energy difference ΔEHS−LS = 3.58 kcal/mol. However, since the spincrossover temperature was obtained in solution, the value in the gas phase, and thus the relative stability of the spin states, might differ. As can be seen from Table 5, and conforming with previous findings,21,29 all three hybrid functionals have trouble predicting the correct ground state of manganocene. B3LYP-D predicts the high-spin 6A1′ state at 3.7 kcal/mol below the 2E2′ low-spin ground state, while the corresponding values for PBE0-D and M06 are even higher: 7.0 and 17.2 kcal/mol, respectively. The tendency of hybrid functionals to overstabilize the high-spin with respect to the low-spin state in manganocene could be attributed to a too high percentage of HF exchange. A possible way to resolve this would then be to reduce the contribution of HF exchange. This has been done, e.g., with B3LYP to improve the description of the spin state energetics in a series of ferrous complexes,19 leading to the B3LYP* functional with only 15% instead of 20% HF exchange. However, no general rules may be formulated in this respect. For example, in a recent study on manganese porphyrins, it was shown that the contribution of HF exchange should in fact be increased to as much as 55% to provide a correct description of the energy gap between the high spin (sextet) and intermediate spin (quartet) state in these complexes.102 The overstabilization in manganocene of the high-spin 6A1′ with respect to the low-spin 2E′ state by the hybrid functionals was also observed by Swart, to such an extent that no spin crossover behavior was predicted. When using the OPBE functional instead, Swart obtained the correct relative spin-state ordering, but with an energy difference of ΔEHS−LS = 8.2 kcal/mol, thus instead overestimating the stability of the doublet state by almost 5 kcal/mol. More accurate results may be obtained from multiconfigurational wave function methods. As Table 5 indicates, with a ΔEHS−LS splitting of 5.77 kcal/mol, CASPT2 correctly predicts a spin crossover compound, overstabilizing the doublet state by 2.2

Table 6. Heterolytic Dissociation Enthalpy (in kcal/mol) of Selected Metallocenes with the B3LYP-D Functional V(Cp)2 A2′

4

ΔE ΔECPC ΔEZPE ΔErel ΔEthermal ΔH298° exptl.16 error

615.6 −5.5 −7.4 2.4 1.9 607.0 606 ± 6 1.0

Mn(Cp)2 A1′

6

581.1 −5.1 −5.7 2.5 0.4 573.2 572 ± 6 1.2

Fe(Cp)2 1

A1′

636.3 −5.6 −8.6 5.5 1.8 629.4 635 ± 6 −5.6

Ni(Cp)2 A2′

3

667.4 −5.1 −6.6 4.1 0.9 660.7 652 ± 6 8.7

Table 7. Heterolytic Dissociation Enthalpy (in kcal/mol) of Selected Metallocenes with the PBE0-Da Functional V(Cp)2 A2′

4

ΔE ΔEdisp ΔECPC ΔEZPE ΔErel ΔEthermal ΔH298° exptl16 error

623.4 7.5 −5.0 −7.3 2.0 2.0 622.6 606 ± 6 16.6

Mn(Cp)2 A1′

6

582.6 6.9 −4.6 −5.6 2.7 1.1 583.1 572 ± 6 11.1

Fe(Cp)2 1

A1′

633.1 11.8 −5.0 −8.7 6.3 1.9 639.4 635 ± 6 4.4

Ni(Cp)2 A2′

3

667.1 8.7 −4.5 −6.6 3.8 1.1 669.6 652 ± 6 17.6

a

PBE0 with dispersion correction as the difference between B3LYP and B3LYP-D.

Table 8. Heterolytic Dissociation Enthalpy (in kcal/mol) of Selected Metallocenes with the M06 Functional V(Cp)2 A2′

4

ΔE ΔECPC ΔEZPE ΔErel ΔEthermal ΔH298° exptl.16 error

630.2 −5.0 −7.0 2.9 1.5 622.6 606 ± 6 16.6

Mn(Cp)2 A1′

6

585.6 −4.7 −5.2 3.4 0.5 579.6 572 ± 6 7.6

Fe(Cp)2 1

A1′

641.1 −5.4 −8.3 6.4 1.9 635.7 635 ± 6 0.7

Ni(Cp)2 A2′

3

677.2 −4.6 −6.6 5.1 1.0 672.1 652 ± 6 20.1

whereas the CASPT2/RASPT2 results may be found in Tables 9 and 10, respectively. An overview of the errors obtained with the different methods is provided by Figure 3. For manganocene, all data refer to the high-spin 6A1′ state. Comparing first the three DFT functionals, we find that (in contrast to the geometry results) B3LYP-D gives much better 889

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

the dissociation enthalpy as compared to CASPT2. For manganocene, ferrocene, and nickelocene, the difference between the binding energy computed with RASPT2 and CASPT2 increases with an increasing number of 3d electrons. This is primarily related to the increasing importance of the 3d double-shell effect in the same order. Due to the limited CAS space, the double-shell effect is only partially included in the CASSCF reference wave function of the CASPT2 treatment. After extending the active space to include a full 3d′ shell, this important correlation effect is more accurately (variationally, up to doubles) treated in the RASSCF wave function, thus providing an improved zeroth-order description leading to more accurate results from the perturbational treatment. The reduction in binding energy is very small (1.1 kcal/mol) in manganocene (d5) but increases to 15.5 kcal/mol in ferrocene (d6) and further to 16.5 kcal/mol in nickelocene (d8). In the latter case, RASPT2 actually seems to overshoot, giving a dissociation enthalpy which is too low (although just slightly outside the error bars of the experimental measurement). These results clearly illustrate the importance of properly dealing with the double-shell correlation effect in systems with a high 3d occupation number. On the other hand, for vanadocene, with just three 3d electrons, the double-shell effect becomes unimportant, but here we rather have to deal with important core−valence correlation effects between the (3s,3p) semicore and 3d valence electrons.87,90 In the CASPT2 calculations on vanadocene, these correlation effects are treated perturbationally in the PT2 step, whereas in RASPT2 they are treated variationally in the reference wave function, by extending the RAS space with the (3s,3p) orbitals. This reduces the binding energy by 5.2 kcal/mol. However, as can be seen from Figure 3, the error on the dissociation enthalpy obtained with RASPT2 for vanadocene is still exceptionally large, almost 10 kcal/mol. The origin of this large discrepancy between RASPT2 and experimental results is at the moment unclear. Increasing the level of excitations of RAS1 to RAS2 to quadruples for the metallocene only affects the dissociation energy by 0.1 kcal/mol, confirming that the RASPT2 calculations are indeed converged. Further investigations on structure, basis sets effects, active spaces, etc. will be needed to resolve this point. However, other than for vanadocene, the results obtained from RASPT2 are in excellent agreement with the experimental data, lying well within the experimental error bars for manganocene and ferrocene, and just outside for nickelocene.

Table 9. Heterolytic Dissociation Enthalpy (in kcal/mol) of Some Metallocenes with CASPT2 Computed at PBE0 Geometry V(Cp)2 4

ΔE ΔECPC ΔEZPEa ΔEthermal ΔH298° exptl.16 error a

A2′

628.5 −2.4 −7.3 2.0 620.8 606 ± 6 14.8

Mn(Cp)2 A1′

6

580.1 −2.0 −5.6 1.1 573.6 572 ± 6 1.6

Fe(Cp)2 A1′

1

664.5 −3.1 −8.7 1.9 654.6 635 ± 6 19.6

Ni(Cp)2 A2′

3

669.9 −2.5 −6.6 1.1 661.9 652 ± 6 9.9

ZPE corrections taken from PBE0.

Table 10. Heterolytic Dissociation Enthalpy (in kcal/mol) of Some Metallocenes with RASPT2 Computed at PBE0 Geometry V(Cp)2 4

ΔE ΔECPC ΔEZPEa ΔEthermal ΔH298° exptl.16 error a

A2′

623.7 −2.8 −7.3 2.0 615.6 606 ± 6 9.6

Mn(Cp)2 A1′

6

579.4 −2.4 −5.6 1.1 572.5 572 ± 6 0.5

Fe(Cp)2 A1′

1

649.7 −3.8 −8.7 1.9 639.1 635 ± 6 4.1

Ni(Cp)2 A2′

3

653.9 −3.0 −6.6 1.1 645.4 652 ± 6 −6.6

ZPE corrections were taken from PBE0.

results for the binding energy than PBE0-D or M06. With the exception of ferrocene, the latter functionals strongly overestimate the dissociation enthalpies. Considering that these two functionals were also found to overestimate the relative stability of the high-spin with respect to the low-spin state (for manganocene, see Table 5), it is quite likely that their good performance for ferrocene is in fact due to a cancellation of errors. We note that ferrocene is the only complex considered in this work for which the dissociation is accompanied by a spin flip on the metal. Overstabilization of the high-spin 5D state of the Fe2+ ion with respect to the low-spin 1A1g state of the ferrocene molecule may counteract a too high “intrinsic” binding energy (i.e., along a spin-conserving path). It is also noteworthy that the results obtained from PBE0-D and M06 for the heterolytic dissociation energy are very similar, with a difference less than 4 kcal/mol for all four molecules. As compared to a previous relativistic corrected DFT study41 using the OPBE functional with a largest error of 48.5 kcal/mol in Ni(Cp)2, the present results for the heterolytic dissociation enthalpy obtained with B3LYP-D are clearly superior. In other B3LYP studies of ferrocene42,43 in which dispersion and relativistic corrections were not included, errors of 20 kcal/mol were observed. With the exception of nickelocene, the remaining errors on the B3LYP-D results are smaller than the experimental uncertainty. For the latter complex, the binding energy is overestimated by 8.7 kcal/mol with B3LYP-D. This error is, however, still considerably smaller than the errors obtained with PBE0-D, +17.4 kcal/mol, and M06, +20.1 kcal/ mol. Considering the multiconfigurational results, we find that CASPT2 systematically overestimates the dissociation enthalpy (although the error for manganocene is quite small, only 1.6 kcal/mol). RASPT2, with its larger active space, always reduces

4. CONCLUSIONS In this work, we have studied the heterolytic dissociation enthalpy of metallocenes of a number of first-row transition metals (V, Mn, Fe, and Ni). Both DFT and high-level ab initio calculations (CASPT2 and RASPT2) were performed. Different basis sets were tested for DFT, and an elaborate analysis was performed of all important contributions to the heterolytic dissociation enthalpy. We showed that it is essential to take into account dispersion and relativistic corrections, and that sizable basis sets and a BSSE correction are necessary for obtaining quantitative results. Among the functionals considered, B3LYP (including an empirical dispersion correction) is obviously superior in providing a consistently accurate description of the heterolytic dissociation enthalpy of the four complexes. PBE0 and M06 also perform well for ferrocene but grossly overestimate the binding energy for the other three metallocenes. This 890

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

(11) Mayor-Lopez, M. J.; Weber, J. Chem. Phys. Lett. 1997, 281, 226− 232. (12) Xu, Z.-F.; Xie, Y.; Feng, W.-L.; Schaefer, H. F. J. Phys. Chem. A 2003, 107, 2716−2729. (13) Coriani, S.; Haaland, A.; Helgaker, T.; Jørgensen, P. ChemPhysChem 2006, 7, 245−249. (14) Pierloot, K.; Persson, B. J.; Roos, B. O. J. Phys. Chem. 1995, 99, 3465−3472. (15) Haaland, A.; Nilsson, J.-E. Chem. Commun. (London) 1968, 88− 89. (16) Ryan, M. F.; Eyler, J. R.; Richardson, D. E. J. Am. Chem. Soc. 1992, 114, 8611−8619. (17) Kozlowski, P. M.; Spiro, T. G.; Bérces, A.; Zgierski, M. Z. J. Phys. Chem. B 1998, 102, 2603−2608. (18) Paulsen, H.; Duelund, L.; Winkler, H.; Toftlund, H.; Trautwein, A. X. Inorg. Chem. 2001, 40, 2201−2203. (19) Reiher, M.; Salomon, O.; Artur Hess, B. Theor. Chem. Acc. 2001, 107, 48−55. (20) Reiher, M. Inorg. Chem. 2002, 41, 6928−6935. (21) Salomon, O.; Reiher, M.; Hess, B. A. J. Chem. Phys. 2002, 117, 4729−4737. (22) Baranović, G. Chem. Phys. Lett. 2003, 369, 668−672. (23) Ghosh, A.; Vangberg, T.; Gonzalez, E.; Taylor, P. R. J. Porphyrins Phtalocyanins 2001, 5, 345−356. (24) Ghosh, A.; Persson, B. J.; Taylor, P. R. J. Biol. Inorg. Chem. 2003, 8, 507−511. (25) Ghosh, A.; Tangen, E.; Ryeng, H.; Taylor, P. R. Eur. J. Inorg. Chem. 2004, 2004, 4555−4560. (26) Paulsen, H.; Trautwein, A. X. Top. Curr. Chem. 2004, 235, 197− 219. (27) Harvey, J. N. Struct. Bonding (Berlin) 2004, 112, 151−183. (28) Deeth, R. J.; Fey, N. J. Comput. Chem. 2004, 25, 1840−1848. (29) Swart, M.; Groenhof, A. R.; Ehlers, A. W.; Lammertsma, K. J. Phys. Chem. A 2004, 108, 5479−5483. (30) Fouqueau, A.; Casida, S. M. M. E.; Daku, L. M. L.; Hauser, A.; Neese, F. J. Chem. Phys. 2004, 120, 9473−9486. (31) Fouqueau, A.; Casida, M. E.; Daku, L. M. L.; Hauser, A.; Neese, F. J. Chem. Phys. 2005, 122, 044110. (32) Daku, L. M. L.; Vargas, A.; Hauser, A.; Fouqueau, A.; Cassida, M. E. ChemPhysChem 2005, 6, 1393−1410. (33) Ganzenmüller, G.; Berkaïne, N.; Fouqueau, A.; Casida, M. E. J. Chem. Phys. 2005, 122, 234321. (34) Pierloot, K.; Vancoillie, S. J. Chem. Phys. 2006, 125, 124303. (35) Strickland, N.; Harvey, J. N. J. Phys. Chem. B 2007, 111, 841− 852. (36) Pierloot, K.; Vancoillie, S. J. Chem. Phys. 2008, 128, 034104. (37) Radoń, M.; Pierloot, K. J. Phys. Chem. A 2008, 112, 11824− 11832. (38) Khvostichenko, D.; Choi, A.; Boulatov, R. J. Phys. Chem. A 2008, 112, 3700−3711. (39) Oláh, J.; Harvey, J. N. J. Phys. Chem. A 2009, 113, 7338−7345. (40) Vancoillie, S.; Radoń, M.; Zhao, H.; Pierloot, K. J. Chem. Theory Comput. 2010, 6, 576−582. (41) Swart, M. Inorg. Chim. Acta 2007, 360, 179−189. (42) Frunzke, J.; Lein, M.; Frenking, G. Organometallics 2002, 21, 3351−3359. (43) Padma Malar, E. Eur. J. Inorg. Chem. 2004, 2004, 2723−2732. (44) Furche, F.; Perdew, J. P. J. Chem. Phys. 2006, 124, 044103. (45) Andersson, K.; Malmqvist, P.-Å.; Roos, B. O. J. Chem. Phys. 1992, 96, 1218−1226. (46) Klopper, W.; Lüthi, H. P. Chem. Phys. Lett. 1996, 262, 546−552. (47) Vancoillie, S.; Zhao, H.; Tran, V. T.; Hendrickx, M. F. A.; Pierloot, K. J. Chem. Theory Comput. 2011, 7, 3961−3977. (48) Malmqvist, P.-Å.; Pierloot, K.; Shahi, A. R. M.; Cramer, C. J.; Gagliardi, L. J. Chem. Phys. 2008, 128, 204109. (49) 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.;

ambiguous behavior has been explained by a cancellation of errors in ferrocene. Dissociation of the two Cp−’s is (only) in this case attended by a spin flip between a low-spin ground state in the complex and a high-spin ground state in Fe2+, and the too high “intrinsic” binding energy predicted by both functionals is in this case counteracted by an overstabilization of the high-spin ionic state. The best metal−cyclopentadienyl distances are obtained with PBE0, while they are slightly underestimated with M06. On the other hand, B3LYP consistently predicts too long metal− cyclopentadienyl distances, even after considering dispersion effects on the structure. Single point CASPT2 and RASPT2 calculations were performed, making use of the PBE0 structures, and employing extensive basis sets so as to minimize the BSSE. With CASPT2, based on an “economical” active space, the binding energy of metallocenes is systematically overestimated. In the three metallocenes M(Cp)2 containing five or more 3d electrons, i.e., with M = Mn, Fe, or Ni, the error may be largely attributed to the inadequate treatment of the double-shell effect with an active space that does not include a full second 3d′ shell. Going to RASPT2 with a more extended active space, this problem is resolved, and the resulting RASPT2 binding energies for these three molecules closely agree with experimental results, with a maximum error of −6.6 kcal/mol for nickelocene. On the other hand, vanadocene behaves more problematic, with a remaining error of almost +10 kcal/mol on the dissociation enthalpy, even after being provided with an improved description of core− valence correlation effects by including the (3s,3p) orbitals in the active space of the RASPT2 treatment. The exceptionally large error for this molecule seems counterintuitive and remains unexplained. The present findings form the basis for an investigation of the binding energy of other metallocenes of second-row transition metals in a future study.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This investigation has been supported by grants from the Flemish Science Foundation (FWO) and from the Concerted Research Action of the Flemish Government (GOA).



REFERENCES

(1) Kealy, T. J.; Pauson, P. L. Nature 1951, 168, 1039−1040. (2) Miller, S. A.; Tebboth, J. A.; Tremaine, J. F. J. Chem. Soc. 1952, 632−635. (3) Birmingham, J. M.; Fischer, A. K.; Wilkinson, G. Naturwissenschaften 1955, 42, 96. (4) Wilkinson, G.; Cotton, F. A.; Birmingham, J. M. J. Inorg. Nucl. Chem. 1956, 2, 95−113. (5) Pfab, W.; Fischer, E. O. Z. Anorg. Allg. Chem. 1953, 274, 316− 322. (6) Wilkinson, G.; Cotton, F. A. Chem. Ind. 1954, 11, 307−308. (7) Park, C.; Almlöf, J. J. Chem. Phys. 1991, 95, 1829−1833. (8) Haaland, A. Top. Curr. Chem. 1975, 53, 1. (9) Fan, L.; Ziegler, T. J. Chem. Phys. 1991, 95, 7401−7408. (10) Matsuzawa, N.; Seto, J.; Dixon, D. A. J. Phys. Chem. A 1997, 101, 9391−9398. 891

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892

Journal of Chemical Theory and Computation

Article

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.; 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 09, Revision A.2; Gaussian Inc.: Wallingford, CT, 2009. (50) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kolmel, C. Chem. Phys. Lett. 1989, 162, 165−169. (51) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (52) Weigend, F.; Häser, M.; Patzelt, H.; Ahlrichs, R. Chem. Phys. Lett. 1998, 294, 143−152. (53) Weigend, F.; Furche, F.; Ahlrichs, R. J. Chem. Phys. 2003, 119, 12753−12762. (54) Dirac, P. A. M. Proc. R. Soc. London, Ser. A 1929, 123, 714−733. (55) Slater, J. C. Phys. Rev. 1951, 81, 385−390. (56) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200− 1211. (57) Becke, A. D. Phys. Rev. A 1988, 38, 3098−3100. (58) Perdew, J. P. Phys. Rev. B 1986, 33, 8822−8824. (59) Cramer, C. J.; Truhlar, D. G. Phys. Chem. Chem. Phys. 2009, 11, 10757−10816. (60) Perdew, J. P.; Wang, Y. Phys. Rev. B 1992, 45, 13244−13249. (61) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (62) Perdew, J. P.; Ernzerhof, M.; Burke, K. J. Chem. Phys. 1996, 105, 9982−9985. (63) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785−789. (64) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (65) Tao, J.; Perdew, J. P.; Staroverov, V. N.; Scuseria, G. E. Phys. Rev. Lett. 2003, 91, 146401. (66) Staroverov, V. N.; Scuseria, G. E.; Tao, J.; Perdew, J. P. J. Chem. Phys. 2003, 119, 12129−12137. (67) Grimme, S. J. Comput. Chem. 2006, 27, 1787−1799. (68) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (69) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553−566. (70) Simon, S.; Duran, M.; Dannenberg, J. J. J. Chem. Phys. 1996, 105, 11024−11031. (71) Hess, B. A. Phys. Rev. A 1986, 33, 3742−3748. (72) Reiher, M.; Wolf, A. J. Chem. Phys. 2004, 121, 2037−2047. (73) Reiher, M.; Wolf, A. J. Chem. Phys. 2004, 121, 10945−10956. (74) Koch, W.; Holthausen, M. C. A Chemist’s Guide to Density Functional Theory, 2nd ed.; Wiley-VCH: New York, 2001; pp 236− 238. (75) van Mourik, T.; Gdanitz, R. J. J. Chem. Phys. 2002, 116, 9620− 9623. (76) Swart, M.; van der Wijst, T.; Fonseca Guerra, C.; Bickelhaupt, F. J. Mol. Model. 2007, 13, 1245−1257. (77) Pitoňaḱ , M.; Neogrády, P.; R̆ ezác,̌ J.; Jurečka, P.; Urban, M.; Hobza, P. J. Chem. Theory Comput. 2008, 4, 1829−1834. (78) Sherrill, C. D.; Takatani, T.; Hohenstein, E. G. J. Phys. Chem. A 2009, 113, 10146−10159. (79) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P.Å.; Neogrády, P.; Pedersen, T. B.; Pitoňaḱ , M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. J. Comput. Chem. 2010, 31, 224−247. (80) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2005, 109, 6575−6579. (81) Roos, B. O.; Lindh, R.; Malmqvist, P.-Å.; Veryazov, V.; Widmark, P.-O. J. Phys. Chem. A 2004, 108, 2851−2858. (82) Widmark, P.-O.; Malmqvist, P.-Å.; Roos, B. O. Theor. Chem. Acc. 1990, 77, 291−306.

(83) Aquilante, F.; Malmqvist, P.-Å.; Pedersen, T. B.; Ghosh, A.; Roos, B. O. J. Chem. Theory Comput. 2008, 4, 694−702. (84) Forsberg, N.; Malmqvist, P.-Å. Chem. Phys. Lett. 1997, 274, 196−204. (85) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. Chem. Phys. Lett. 2004, 396, 142−149. (86) Andersson, K.; Roos, B. O. Chem. Phys. Lett. 1992, 191, 507− 514. (87) Roos, B. O.; Andersson, K.; Fülscher, M. P.; Malmqvist, P.-Å.; Serrano-Andrés, L.; Pierloot, K.; Merchán, M. In Advances in Chemical Physics: New Methods in Computational Quantum Mechanics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons: New York, 1996; Vol. XCIII, pp 219−332. (88) Pierloot, K. In Computational Organometallic Chemistry; Cundari, T. R., Ed.; Marcel Dekker, Inc.: New York, 2001; pp 123− 158. (89) Pierloot, K. Int. J. Quantum Chem. 2011, 111, 3291−3301. (90) Pierloot, K.; Tsokos, E.; Roos, B. O. Chem. Phys. Lett. 1993, 214, 583−590. (91) Takusagawa, F.; Koetzle, T. F. Acta Crystallogr., Sect. B 1979, 35, 1074−1081. (92) Brock, C. P.; Fu, Y. Acta Crystallogr., Sect. B 1997, 53, 928−938. (93) Jiménez-Hoyos, C. A.; Janesko, B. G.; Scuseria, G. E. J. Phys. Chem. A 2009, 113, 11742−11749. (94) Jensen, K. P.; Roos, B. O.; Ryde, U. J. Chem. Phys. 2007, 126, 014103. (95) Reddy, A. S.; Vijay, D.; Sastry, G. M.; Sastry, G. N. J. Phys. Chem. B 2006, 110, 2479−2481. (96) Switzer, M. E.; Wang, R.; Rettig, M. F.; Maki, A. H. J. Am. Chem. Soc. 1974, 96, 7669−7674. (97) Köhler, F. H.; Schlesinger, B. Inorg. Chem. 1992, 31, 2853− 2859. (98) Haaland, A. Inorg. Nucl. Chem. Lett. 1979, 15, 267−269. (99) Gard, E.; Haaland, A.; Novak, D.; Seip, R. J. Organomet. Chem. 1975, 88, 181−189. (100) Haaland, A.; Lusztyk, J.; Novak, D. P.; Brunvoll, J.; Starowieyski, K. B. J. Chem. Soc., Chem. Commun. 1974, 54−55. (101) Hedberg, L.; Hedberg, K. J. Chem. Phys. 1970, 53, 1228−1234. (102) Kepenekian, M.; Calborean, A.; Vetere, V.; Le Guennic, B.; Robert, V.; Maldivi, P. J. Chem. Theory Comput. 2011, 7, 3532−3539.

892

dx.doi.org/10.1021/ct200875m | J. Chem. Theory Comput. 2012, 8, 883−892