Article pubs.acs.org/JPCA
Unrestricted Prescriptions for Open-Shell Singlet Diradicals: Using Economical Ab Initio and Density Functional Theory to Calculate Singlet−Triplet Gaps and Bond Dissociation Curves Daniel H. Ess* and Thomas C. Cook Department of Chemistry and Biochemistry, Brigham Young University, Provo, Utah, 84602, United States S Supporting Information *
ABSTRACT: Here we present and test several computational prescriptions for calculating singlet−triplet (ST) gap energies and bond dissociation curves for open-shell singlet diradicals using economical unrestricted single reference type calculations. For ST gap energies from Slipchenko and Krylov’s atom and molecule test set (C, O, Si, NH, NF, OH+, O2, CH2, and NH2+) spin unrestricted Hartree−Fock and MP2 energies result in errors greater than 15 kcal/mol. However, spin−projected (SP) Hartree−Fock theory in combination with spincomponent-scaled (SCS) or scaled-opposite-spin (SOS) second-order perturbation theory gives ST gap energies with a mean unsigned error (MUE) of less than 2 kcal/mol. Density functionals generally give poor results for unrestricted energies and only the ωB97X-D, the M06, and the M06-2X functionals provide reasonable accuracy after spin-projection with MUE values of 4.7, 4.3, and 3.0 kcal/mol, respectively, with the 6-311++G(2d,2p) basis set. We also present a new one parameter hybrid density functional, diradical-1 (DR-1), based on Adamo and Barone’s modified PW exchange functional with the PW91 correlation functional. This DR-1 method gives a mean error (ME) of 0.0 kcal/mol and a MUE value of 1.3 kcal/mol for ST gap energies. As another test of unrestricted methods the bond dissociation curves for methane (CH4) and hydrofluoric acid (H−F) were calculated with the M06-2X, DR-1, and ωB97X-D density functionals. All three of these functionals give reasonable results for the methane C−H bond but result in errors greater than 50 kcal/mol for the H−F bond dissociation. Spin-projection is found to significantly degrade bond dissociation curves past ∼2.2 Å. Although unrestricted Hartree−Fock theory provides a very poor description of H−F bond dissociation, unrestricted SCS-MP2 and SOS-MP2 methods give accurate results.
1. INTRODUCTION Open-shell singlet diradical molecules1 are often treated with density functional methods using a broken symmetry spin unrestricted formalism2 where opposite spin electrons have separate Kohn−Sham (KS) orbitals that are allowed to polarize. Although this type of calculation results in a spin contaminated “DFT wavefunction” that is an admixture of multiple spin states and is not an eigenfunction of the S2 operator, the energy is typically superior to that of a closed-shell solution.3 In addition, density functional calculations are typically computationally very efficient compared to multireference ab initio methods. However, it remains unclear how accurate spin contaminated values are for commonly used functionals. We recently found that using unrestricted calculations with the PBEPBE and LSDA functionals resulted in errors greater than 22 kcal/mol for the singlet−triplet (ST) gap energies (ΔEST) of Slipchenko and Krylov’s4 atom and molecule set of diradicals.5 Although application of the Yamaguchi spin-projection technique6,7 or the Ziegler, Rauk, and Baerends (ZRB) energy correction8 reduced this error by about half, the errors remain too large to suggest either of these functionals for general use. One solution to improve the accuracy of calculating ST gaps with the PBEPBE and LSDA functionals is to use the fractional spin method.5,9 This method uses an ensemble of closed-shell densities10 to give the correct open-shell singlet density. For the © 2012 American Chemical Society
PBEPBE and LSDA functionals ST gaps were modeled with errors of less than 3 kcal/mol.5 However, a severe limitation to the fractional spin approach is that popular hybrid density functionals, such as B3LYP, result in significant errors because of the inclusion of Hartree−Fock exchange.9,11 Because the PBEPBE and LSDA functionals perform poorly for general use, the question arises whether any currently available functional can accurately treat open-shell singlet diradicals. This is especially important since many diradical species are too large to be treated with correlated multireference ab initio methods.12 It is important to note that an alternative methodology based on the spin-flip procedure and TDDFT can be used economically and accurately to compute ST gaps.13 There are only a few studies exploring the accuracy of recently developed density functionals for modeling ST gap energies of open-shell singlet species. Most notable are the recent studies by Saito et al. that compared B3LYP, CAMB3LYP, BHandHLYP, and LC-BLYP functionals for calculating ST gap energies of trimethylenemethane, iminoallyl, and oxyallyl diradicals, as well as antiaromatic molecules.14 Received: January 18, 2012 Revised: April 5, 2012 Published: May 11, 2012 4922
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
Table 1. ST Gap Energies, ME, and MUE for Spin Unrestricted Hartree-Fock and MP2 Methods with the 6-311++G(2d,2p) Basis Set (kcal/mol) C O Si NH NF OH+ O2 CH2 NH2+ ME MUE a
exp.a
HF
29.1 45.4 17.3 35.9 34.3 50.5 22.6 32.9 44.6
14.6 23.2 7.6 20.0 20.4 26.6 17.9 21.4 25.8
29.5 46.5 15.8 40.2 40.8 53.2 35.9 42.9 51.4
16.8 24.4 11.5 20.0 18.5 26.8 4.1 18.4 23.5
31.5 47.7 19.6 40.2 38.9 53.4 23.7 39.9 49.1
14.4 21.5 10.2 16.9 16.0 23.7 3.3 15.0 20.0
29.1 44.8 18.4 37.0 36.4 50.4 22.9 36.5 45.7
13.2 20.0 9.6 15.3 14.8 22.2 2.9 13.3 18.3
27.9 43.3 17.7 35.5 35.2 48.8 22.5 34.8 43.9
−15.0 15.0
5.2 5.2
−16.5 16.5
3.5 3.5
−19.1 19.1
0.9 1.1
−20.3 20.3
−0.3 1.0
SP-HF
MP2
SP-HF-MP2
SCS-MP2
SP-HF-SCS-MP2
SOS-MP2
SP-HF-SOS-MP2
See ref 4 for reference to experimental values and CASSCF SOCI value for NH2+.
Here we present our investigation to identify reasonably economical spin unrestricted single reference ab initio and density functional methods capable of accurately modeling ST gaps and bond dissociation curves that involve open-shell singlet diradicals. Section 3.1 reports our analysis of Hartree− Fock and MP2 theory with spin-projection and spin scaled PT2 corrections. Section 3.2 tests the accuracy of density functionals with and without spin-projection. In Section 3.3 we present a new one parameter hybrid density functional (DR-1) that gives accurate ST gap energies. Section 3.4 investigates the accuracy of M06-2X and DR-1 methods for open-shell and closed-shell states of divalent species. Lastly, in Section 3.5 we evaluate the performance of density functional and MP2 methods for modeling the bond dissociation curves of methane and hydrofluoric acid.
for only OS correlation energy along with complete neglect of SS correlation energy. This method has been termed scaledopposite-spin MP2 (SOS-MP2).17 E SCS = E HF + α(ESS) + β(EOS)
2.1. ST Gap Energies from Economical Single Reference Ab Initio Methods. Single reference Hartree− Fock theory is typically the starting point for perturbation and multireference treatments. Although Hartree−Fock theory is not expected to perform well because of the lack of electron correlation, there are some types of reactions, for example, isodesmic reactions, that can be treated because of error cancellation. Table 1 gives the spin unrestricted Hartree−Fock ST gap energies (ΔEST) for the Slipchenko and Krylov atom and molecule test set that includes C, O, Si, NH, NF, OH+, O2, CH2, and NH2+.4 The positive energy values in Table 1 correspond to the triplet ground state and the open-shell singlet excited state.18 The unrestricted Hartree−Fock ΔEST values are significantly underestimated on average by ∼15 kcal/mol. Unexpectedly, however, spin-projected Hartree−Fock theory (SP-HF) gives reasonable ΔEST values for these atoms and molecules with a mean unsigned error (MUE) of 5.2 kcal/mol. Unrestricted MP2 and spin-component-scaled MP2 (SCSMP2) energies show degradation in accuracy compared to unrestricted Hartree−Fock theory. The MUE values are 16.5 and 19.1 kcal/mol, respectively (Table 1). Accurate ST gap energies can be estimated by combining a spin−projected Hartree−Fock reference energy with PT2 correlation energies. Using unscaled PT2 energies (SP-HF-MP2) results in ΔEST energies with a MUE value of 3.5 kcal/mol while the use of Grimme’s scaling factors16 for same spin and opposite spin correlation energies (SP-HF-SCS-MP2) gives a MUE value of 1.1 kcal/mol. This demonstrates that after the singlet Hartree− Fock reference energy is spin-projected PT2 energy corrections improve accuracy by inclusion of dynamic correlation energy. The use of Head-Gordon’s 1.3 scale factor17 for only opposite spin correlation energies and neglect of same spin correlation energy combined with a spin−projected Hartree−Fock reference energy (SP-HF-SOS-MP2) gives errors (ME = −0.3 kcal/mol; MUE = 1.0 kcal/mol) nearly identical to higher Grimme’s suggested scaling. 2.2. ST Gap Energies from Density Functional Methods. Density functional theory provides the most economical platform for analysis of open-shell singlet diradicals and reaction mechanisms that involve diradical intermediates
2. METHODOLOGY DETAILS Ab initio and density functional singlet and triplet molecular geometries were optimized using GAUSSIAN 09 with either the 6-31G(d,p) or 6-311++G(2d,2p) basis set.15 Geometric values are given in the Supporting Information. All reported ST gap energies (ΔEST = ET − ES) for molecules correspond to optimized singlet and triplet geometries. All singlet atoms and molecules have ⟨S2⟩ values of 1.0 while triplet species have ⟨S2⟩ values of 2.0. The spin−projection technique used follows that suggested by Yamaguchi and co-workers.6 In this procedure the spin-projected singlet energy (ESPS) is obtained by adding to the contaminated singlet energy (ES) the difference in energy between the singlet and triplet energies scaled by χ, which is a weighting ratio of singlet to triplet spin contamination values (eqs 1−2). For open-shell singlets with ⟨S2⟩ values of 1.0 spinprojection gives equivalent results to the ZRB energy correction.8 E SPS = ES + χ[ESE T]
(1)
χ = ( S2 S / S2 T )/[1 − ( S2 S / S2 T )]
(2)
(3)
For MP2 calculations Grimme has introduced the idea of scaling the spin components of PT2 dynamic correlation energy.16 This scaling provides a convenient way to systematically improve MP2 energies. The original spin-componentscaled MP2 method (SCS-MP2) scales the same spin (SS) correlation energies by α = 1/3 and opposite spin (OS) correlation energies by β = 6/5 (eq 3). Head-Gordon has suggested a more economical method using a 1.3 scale factor 4923
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
Table 2. ST Gap Energies, ME, and MUE for Spin Unrestricted Density Functional Methods with the 6-311++G(2d,2p) Basis Set (kcal/mol) MPW1PW91
CAM-B3LYP
BH&HLYP
BMK
B2PLYP
ωB97X-D
B3LYP
X3LYP
C O Si NH NF OH+ O2 CH2 NH2+
9.0 16.5 5.1 13.0 11.9 18.4 10.2 11.6 15.8
9.1 16.5 5.1 13.0 12.0 18.5 10.4 11.5 15.8
9.2 16.7 6.5 13.6 11.2 18.6 8.9 13.0 16.7
9.9 17.0 5.1 13.5 12.5 19.0 10.6 12.0 16.4
9.3 17.1 4.8 13.4 13.2 19.2 12.0 12.2 16.7
10.8 17.2 6.7 15.0 12.9 19.0 10.7 14.8 18.0
11.0 18.6 6.7 14.8 13.7 20.8 10.8 13.2 17.9
12.7 19.3 6.3 15.9 14.7 21.5 11.7 14.9 19.2
M06-2X 14.6 20.5 6.8 19.0 18.7 24.1 14.2 18.4 23.6
M06 12.4 24.6 7.1 19.8 18.6 28.6 15.8 15.8 24.6
ME MUE
−22.3 22.3
−22.3 22.3
−22.0 22.0
−21.8 21.8
−21.6 21.6
−20.8 20.8
−20.6 20.6
−19.6 19.6
−17.0 17.0
−16.2 16.2
Table 3. Spin-Projected ST Gap Energies, ME, and MUE for Spin Unrestricted Density Functional Methods with the 6-311+ +G(2d,2p) Basis Set (kcal/mol) B3LYP
X3LYP
MPW1PW91
CAM-B3LYP
BH&HLYP
BMK
B2PLYP
ωB97X-D
M06
M06-2X
C O Si NH NF OH+ O2 CH2 NH2+
18.1 32.9 10.3 26.0 23.8 36.9 20.5 23.2 31.6
18.2 33.0 10.3 26.0 24.0 37.0 20.7 23.1 31.6
18.6 33.5 13.1 27.2 22.5 37.1 17.8 26.1 33.3
19.9 33.9 10.3 26.9 25.1 38.0 21.2 24.1 32.7
18.8 34.2 9.8 26.9 26.4 38.5 24.1 24.5 33.3
21.7 34.4 13.6 30.0 25.7 38.0 21.4 29.7 35.9
22.2 37.3 13.7 29.7 27.5 41.5 21.6 26.5 35.7
25.4 38.6 12.8 31.8 29.4 42.9 23.4 29.7 38.4
24.8 49.1 14.2 39.5 37.2 57.1 31.5 31.6 49.1
29.1 41.0 13.8 37.9 37.3 48.1 28.3 36.8 47.0
ME MUE
−9.9 9.9
−9.9 9.9
−9.3 9.3
−8.9 8.9
−8.5 8.8
−6.9 6.9
−6.3 6.3
−4.5 4.7
3.3 4.3
0.7 3.0
mol. The M06 and M06-2X functionals slightly overestimate the ST gaps with average errors of 3.3 and 0.7 kcal/mol. Their MUE values are 4.3 and 3.0 kcal/mol, respectively. The results presented in Tables 2 and 3 utilized the large 6311++G(2d,2p) basis set. However, this basis set is often intractable for large molecules. Therefore, we also tested the B3LYP, BHandHLYP, and M06-2X methods with the medium sized 6-31G(d,p) basis set (Table 4). These functionals show a slight uptick in their performance. Overall, spin-projected M062X values with a medium sized basis set provide a reasonable
(or transition states). Over the past several years there has been an intense surge in the development of new density functionals. Most notable is the M06 series by Zhao and Truhlar,19 Yanai et al.’s long-range corrected B3LYP method (CAM-B3LYP),20 the modified B97 functionals by Chai and Head-Gordon,21 and Boese and Martin’s BMK functional.22 These functionals along with the popular B3LYP, 2 3 the X3LYP variant, 2 4 MPW1PW91,25 and BHandHLYP (BH&HLYP) functionals were evaluated for modeling ST gap energies. Table 2 gives the spin unrestricted energy results of nine pure and hybrid functionals along with the B2PLYP26 double hybrid method.27 It is clear from Table 2 that spin contaminated singlet energies are highly underestimated. For example, the popular B3LYP functional has a mean error (ME) of −22.3 kcal/mol. The long-range corrected version CAM-B3LYP performs just as poorly with a ME of −21.8 kcal/mol. There is slight improvement with the use of M06-2X and M06 functionals, which show MEs of −17.0 and −16.2 kcal/mol, respectively. Clearly, unrestricted density functional energies cannot be suggested for use to model ST gap energies with open-shell singlets. In contrast, spin-projected density functional ΔEST values show a range of success. The B3LYP, X3LYP, MPW1PW91, CAM-B3LYP, and BHandHLYP functionals still underestimate the ΔEST gap energies by over 8.5 kcal/mol. BMK and the double hybrid B2PLYP method show average errors of ∼6−7 kcal/mol. Three functionals show reasonable accuracy. The ωB97X-D functional underestimates the ST gaps by an average of 4.5 kcal/mol and the MUE is only slightly worse at 4.7 kcal/
Table 4. Spin-Projected ST Gap Energies, ME, and MUE for Spin Unrestricted Density Functional Methods with the 631G(d,p) Basis Set (kcal/mol)
4924
B3LYP
BH&HLYP
M06-2X
C O Si NH NF OH+ O2 CH2 NH2+
19.7 33.9 11.3 27.2 24.2 37.6 20.9 26.0 32.2
20.5 35.1 11.2 28.2 27.0 39.3 24.5 27.6 33.7
31.5 41.6 14.6 39.6 37.5 48.6 28.1 40.5 47.2
ME MUE
−8.8 8.8
−7.3 7.7
1.8 3.7
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
Table 5. Spin-Projected ST Gap Energies for the DR-1 Method, ME, and MUE (kcal/mol) C O Si NH NF OH+ O2 CH2 NH2+ ME MUE a
exp.a
DR-1b
SP-DR-1b
SP-DR-1c
SP-DR-1d
SP-DR-1e
SP-DR-1f
29.1 45.4 17.3 35.9 34.3 50.5 22.6 32.9 44.6
15.1 22.1 8.9 18.4 15.9 24.6 12.6 16.8 21.8 −17.4 17.4
30.3 44.3 17.9 36.8 31.9 49.3 25.3 33.6 43.5 0.0 1.3
31.8 44.7 18.8 38.0 32.0 49.7 25.6 36.4 43.8 0.9 1.9
31.0 44.4 18.3 38.0 32.0 49.4 25.3 35.0 43.8 0.5 1.7
29.8 43.7 17.0 36.4 31.3 48.5 25.1 33.3 42.6 −0.5 1.5
29.5 43.7 16.7 36.0 31.3 48.4 25.0 32.8 42.4 −0.8 1.4
See ref 4. b6-311++G(2d,2p). c6-31G(d,p). dcc-pVDZ. ecc-pVTZ. faug-cc-pVTZ.
singlet (1A1),30 and triplet (3B1) energies for CH2 as a function of the HCH angle. The M06-2X and DR-1 methods give highly similar curves for the triplet and open-shell singlet states. However, there is a difference for the closed-shell singlet curves.
and likely reliable protocol for modeling ST gap energies that involve open-shell singlet diradicals. 3.3. DR-1, One Parameter Hybrid Density Functional Approximation for Singlet Diradicals. We have designed a new functional, diradical-1 (DR-1), that gives highly accurate energy differences between triplet and open-shell singlet states. This method combines robust and well-known exchange and correlation functionals in a pseudo one parameter hybrid scheme. Generally, one parameter hybrid Fock−Kohn−Sham exchange-correlation energies can be represented by eq 4. HF Exc = PE + (1 − P1)(Ex LSD + ΔEx GGA ) + Ec LSD 1 x
+ ΔEcGGA
(4)
Adamo and Barone used this one parameter approach to give the MPW1PW91 (also labeled as mPW1PW91) method that combines a modified PW exchange functional with the PW91 correlation functional.25 Truhlar and co-workers have developed similar one parameter methods such as the MPW1K method that provides accurate thermochemical kinetics.28,29 The spin-projected MPW1PW91 functional gives a MUE value of 9.3 kcal/mol for ST gap energies (Table 3, Section 3.2). A significantly more accurate functional is achieved with a simple modification of this one parameter scheme where the P1 parameter only controls the relative amounts of LSD and HF exchange (eq 5). The diradical-1 (DR-1) method uses Adamo and Barone’s modified PW exchange functional along with the PW91 correlation functional and P1 is set to 0.2 (eq 6). Table 5 gives the spin unrestricted and spin-projected ST gap energies computed with the DR-1 method. Gratifyingly, the spinprojected DR-1 functional gives a ME of 0.0 kcal/mol and MUE of 1.3 kcal/mol with the 6-311++G(2d,2p) basis set. A small error for spin-projected ST gaps is also found with medium sized 6-31G(d,p) and cc-pVDZ basis sets and larger cc-pVTZ and aug-cc-pVTZ basis sets (Table 5). HF Exc = PE + (1 − P1)(Ex LSD) + ΔEx GGA + Ec LSD 1 x
+ ΔEcGGA
(5)
Exc = 0.2Ex HF + 0.8Ex LSD + ΔEx mPW + Ec LSD + ΔEc PW91
Figure 1. Triplet (3B1), open-shell singlet (1B1), spin-projected openshell singlet (1B1−SP), restricted closed-shell singlet (1A1) and spinprojected unrestricted singlet (1A1-SP) electronic energy for CH2 as a function of HCH angle from 90° to 180°. C−H distances were fixed at their optimized bond length. The top plot is for the M06-2X method with the 6-311++G(2d,2p) basis set. The bottom plot is for the DR-1 method with the 6-311++G(2d,2p) basis set. (Energies plotted in Hartree).
(6)
3.4. Spin State Energies of Divalent XH2 species (X = C, Si, P). Analysis of CH2, SiH2, and PH2+ species by the M062X and DR-1 methods allows testing of closed-shell singlet states relative to triplet and open-shell singlet states. Figure 1 shows the calculated open-shell singlet (1B1), closed-shell 4925
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
Figure 2. Left (M06-2X/6-31G(d)) and right (DR-1/6-31G(d)) methane bond dissociation curves plotted from 0.8 Å to 2.8 Å. Vertical axis is plotted in kcal/mol. The full CI energy values were taken from ref 38.
mol above the closed-shell singlet energy, and the spinprojected singlet energy is 43.3 kcal/mol above the ground state. For PH2+ the DR-1 method again underestimates the energy difference between the closed-shell singlet and triplet (9.8 kcal/ mol) and the open-shell singlet (35.9 kcal/mol). However, the DR-1 method gives a 26.1 kcal/mol energy difference between the lowest energy triplet state and the open-shell singlet. This is only 0.9 kcal/mol below the experimental value. Analysis of CH2, SiH2, and PH2+ species confirms that M062X is highly accurate at modeling all spin states while the DR-1 method shows some deficiency for modeling closed-shell singlets. As with any functional this demonstrates that the DR-1 method may not be accurate for broad general use. However, the development of this method does demonstrate that an accurate functional can be formulated within the platform of spin-projected density functional theory to model open-shell singlet diradicals that have strong nondynamical correlation effects.34−36 3.5. CH4 and H−F Bond Dissociation Energy Curves. Singlet diradicals are also generated as closed-shell two electron bonds are dissociated. Methane and H−F provide two ideal comparisons because the methane CH bond dissociates into a hydrogen radical and localized methyl radical while the hydrofluoric acid bond dissociates into a hydrogen radical and fluorine radical with degenerate occupied orbitals.37 Dutta and Sherrill have computed the full configuration interaction (FCI) limit for the methane C−H bond dissociation curve with the 6-31G(d) basis set and the H−F bond dissociation curve with the 6-31G(d,p) basis set.38 They have also shown that for methane spin unrestricted Hartree−Fock and MP2 methods result in errors of approximately 83 and 14 kcal/mol at the complete dissociation limit. In the Dutta and Sherrill study the only density functional method tested was B3LYP.38 Figure 2 shows the methane bond dissociation curves from 0.8 Å to 2.8 Å for the M06-2X and DR-1 methods compared with the FCI results from Dutta and Sherrill.38 For the M06-2X method, the ⟨S2⟩ value is nonzero after 2.2 Å (0.5) and increases to 0.9 at 2.8 Å. Similarly, the DR-1 method shows an onset of spin contamination with an ⟨S2⟩ value of 0.3 at 2.2 Å that increases to 0.9 at 2.8 Å. At 2.8 Å the FCI dissociation energy is 106.2 kcal/mol.38 The M06-2X spin contaminated
At a 90° angle the M06-2X method gives a closed-shell singlet energy that is lower than the triplet energy. The closed-shell singlet curve crosses with the triplet energy curve at ∼105°, the open-shell singlet curve at ∼130° angle, and the spin-projected curve at ∼160° angle. In contrast, the DR-1 method gives a triplet energy that is lower than the closed-shell singlet energy at a 90° angle. However, the spin-projected unrestricted DR-1 curve shows the correct ordering. The difference in where the singlet-triplet spin state crossings occur can be attributed to how the M06-2X and DR-1 methods model the closed-shell singlet state. The M06-2X method gives a 13.0 kcal/mol ST gap between the triplet and the closed-shell singlet state. This is an overestimation by 4.0 kcal/mol compared with experiment (9.0 kcal/mol).31 The DR-1 method overestimates this ST gap by 11.0 kcal/mol relative to experiment. The SiH2 and PH2+ divalent species have closed-shell singlet ground states. The M06-2X method gives an energy difference between the SiH2 closed-shell singlet and triplet state of 18.8 kcal/mol. This is 2.2 kcal/mol lower than the experimental value of 21.0 kcal/mol.32 The spin-projected M06-2X ST gap between the ground state and the open-shell singlet state is 42.4 kcal/mol, which is only 2.1 kcal/mol lower than experiment.4 The spin-projected M06-2X prescription also accurately models the ST gap energy between the triplet state and the open-shell singlet state to within 0.1 kcal/mol. The DR-1 method gives the energy difference between the closed-shell singlet and the triplet states of SiH2 to be 13.5 kcal/mol. This is an underestimation of the experimental value by 7.5 kcal/mol because the closed-shell singlet energy is too high. The spin-projected DR-1 energy difference between the ground state and the open-shell singlet state is 36.4 kcal/mol. Again, this ST gap energy is also underestimated by 8.0 kcal/ mol. However, the spin-projected DR-1 energy difference between the triplet and the open-shell singlet states is 22.9 kcal/mol, which is only 0.7 kcal/mol smaller than experiment. The experimental energy difference between the ground state singlet and triplet states of PH2+ is 17.3 kcal/mol while the energy gap with the open-shell singlet is 44.3 kcal/mol.33 This gives an experimental triplet to open-shell singlet gap energy of 27.0 kcal/mol. The M06-2X method models all of these energy gaps accurately. The triplet energy is computed to be 15.4 kcal/ 4926
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
Figure 3. Left (M06-2X/6-31G(d,p)) and right (DR-1/6-31G(d,p)) H−F bond dissociation curves plotted from 0.7 Å to 2.8 Å. Vertical axis is plotted in kcal/mol. The full CI energy values were taken from ref 38.
neglect of SS correlation or to the specific OS correlation scaling value, we have calculated the energy of dissociation at 2.8 Å varying the OS scaling factor from 0.0 to 2.0 for when there is complete neglect of SS correlation energy and when the full SS correlation energy is used. Figure 4 shows that the OS
value is 106.9 kcal/mol. Spin-projection lowers this value to 101.6 kcal/mol.39 The DR-1 spin contaminated dissociation energy is 115.5 kcal/mol. After spin-projection the energy is lowered to 111.5 kcal/mol. Although not shown in Figure 1, the ωB97X-D functional gave spin contaminated and spinprojected dissociation energies of 111.5 and 107.4 kcal/mol. Clearly, all three of these density functionals give spin contaminated and spin-projected dissociation energy values within reasonable agreement of the FCI values. In contrast to the methane C−H bond dissociation, the H−F bond dissociation is significantly more problematic to model with density functional methods. Figure 3 shows the M06-2X and DR-1 bond dissociation curves for H−F. Although spin contamination does not set in until a partial bond length of 2.2 Å, the density functional curves significantly deviate from the FCI curve starting at ∼1.4 Å. At 2.8 Å the FCI dissociation energy is 126.8 kcal/mol. The spin unrestricted M06-2X energy at 2.8 Å is 192.9 kcal/mol, an error of 66.1 kcal/mol. Unfortunately, the spin-projected energy increases this energy to 208.6 kcal/mol. This highlights that when a bond dissociates into fragments where at least one of the radical fragments is expected to have degenerate occupied orbitals, spin-projection will likely lead to worse results than spin contaminated energies. The spin unrestricted DR-1 dissociation energy at 2.8 Å is 178.9 kcal/mol. This value has an error of 52.1 kcal/mol compared to FCI. Although this is an improvement over the M06-2X energy, it is still not reasonable. Again, spin-projection of the singlet DR-1 energies results in a degradation of the predicted dissociation curve. Since all density functionals, including our newly devised DR-1, fail to model the H−F dissociation curve, we now turn our attention to Hartree−Fock and MP2 methods with spin scaling but without spin-projection. Spin unrestricted Hartree− Fock and MP2 dissociation energies at 2.8 Å are 228.8 and 151.5 kcal/mol, respectively. This is an error of 102 and 25 kcal/mol. Using Grimme’s SCS-MP2 method gives an improved value of 139.8 kcal/mol with an error of 13.0 kcal/ mol compared to the FCI energy. Importantly, Head-Gordon’s opposite-spin (OS) scaling value of 1.3 and neglect of same spin (SS) correlation gives a dissociation energy value of 129.7 kcal/mol that is only 2.9 kcal/mol higher than the FCI value. To explore if this accurate dissociation energy value is due to
Figure 4. H−F bond dissociation energy at 2.8 Å as a function of OS scaling factor from 0.0 to 2.0. The nonspin-projected SOS-MP2 energies corresponds to complete neglect of SS correlation energy. The SOS-SS energies correspond to scaling on OS correlation energy and inclusion of the full SS correlation energy (kcal/mol).
scaling factor is most important. When there is complete neglect of SS correlation energy, the FCI energy value matches the MP2 value when the OS correlation energy is scaled by Head-Gordon’s suggested value. However, when the full SS correlation energy is included, the MP2 energy matched the FCI value when the OS correlation is scaled by a factor closer to 1.4.
4. CONCLUSION Here we have tested and developed computationally economical protocols for using spin unrestricted type calculations to model ST gaps and bond dissociation curves that involve open-shell singlet diradicals. Major conclusions are as follows: • For ab initio methods, SP-HF theory gives errors of ∼5 kcal/mol for ST gap energies. Highly accurate results with errors of less than 2 kcal/mol can be obtained using spin projected Hartree−Fock theory with scaled PT2 4927
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
• • • •
•
■
Article
T. Theor.Chim. Acta 1978, 48, 185. (c) Yamaguchi, K.; Takahara, Y.; Fueno, T.; Houk, K. N. Theor. Chim. Acta 1988, 73, 337. (d) Yamaguchi, K.; Kawakami, T.; Nagao, H.; Yamaguchi, K. Chem. Phys. Lett. 1994, 231, 25. For more recent studies see: (e) Kitagawa, Y.; Saito, T.; Ito, M.; Shoji, M.; Koizumi, K.; Yamanaka, S.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2007, 442, 445. (f) Kitagawa, Y.; Saito, T.; Ito, M.; Nakanishi, Y.; Shoji, M.; Koizumi, K.; Yamanaka, S.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Int. J. Quantum Chem. 2007, 107, 3094. (g) Saito, T.; Kitagawa, Y.; Shoji, M.; Nakanishi, Y.; Ito, M.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2008, 456, 76. (h) Saito, T.; Nishihara, S.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2009, 483, 168. (i) Kitagawa, Y.; Saito, T.; Nakanishi, Y.; Kataoka, Y.; Shoji, M.; Koizumi, K.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Int. J. Quantum Chem. 2009, 109, 3641. (j) Saito, T.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. 2010, 368, 1. (7) Wittbrodt, J. M.; Schlegel, H. B. J. Chem. Phys. 1996, 105, 6574. (8) Ziegler, T.; Rauk, A.; Baerends, E. J. Theoret. Chim. Acta. 1977, 43, 261. (9) (a) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. J. Chem. Phys. 2008, 129, 121104. (b) Cohen, A. J.; Mori-Sánchez, P.; Yang, W. Science 2008, 321, 792. (c) Zeng, X.; Hu, H.; Hu, X.; Cohen, A. J.; Yang, W. J. Chem. Phys. 2008, 128, 124510. (10) (a) Slater, J. C.; Mann, J. B.; Wilson, T. M.; Wood, J. H. Phys. Rev. 1969, 184, 672. (b) Levy, M. Phys. Rev. A 1982, 26, 1200. (c) Dunlap, B. I.; Mei, W. N. J. Chem. Phys. 1983, 78, 4997. (d) Gross, E. K. U.; Oliveria, L. N.; Kohn, W. Phys. Rev. A 1988, 37, 2809. (e) Averill, F. W.; Painter, G. S. Phys. Rev. B 1992, 46, 2498. (f) Zhao, Q.; Morrison, R. C.; Parr, R. G. Phys. Rev. A 1994, 50, 2138. (g) Becke, A. D.; Savin, A.; Stool, H. Theor. Chem. Acc. 1995, 91, 147. (h) Cramer, C. J.; Dulles, F. J.; Giesen, D. J.; Almlöf, J. Chem. Phys. Lett. 1995, 245, 165. (i) Valiev, M. M.; Fernando, G. W. Phys. Rev. B 1995, 52, 10697. (j) Wang, S. G.; Schwarz, W. H. E. J. Chem. Phys. 1996, 105, 4641. (k) Schipper, P. R. T.; Gritsenko, O. V.; Baerends, E. J. Theor. Chem. Acc. 1998, 99, 329. (l) Schipper, P. R. T.; Gritsenko, O. V.; Baerends, E. J. J. Chem. Phys. 1999, 111, 4056. (m) Filatov, M.; Shaik, S. J. Phys. Chem. A 2000, 104, 6628. (n) Visser, S. P. d; Filatov, M.; Shaik, S. Phys. Chem. Chem. Phys. 2000, 2, 5046. (o) Pérez-Jiménez, Á . j.; PérezJordá, J. M. Phys. Rev. A 2007, 75, 012503. (11) Another drawback of the fractional spin approach is that most commercial software do not allow self-consistent calculation with fractional occupation. (12) (a) Wang, J.; Yu, Z. Y.; Philpott, M. R.; Vukovic, S.; Lester, W. A.; Cui, T.; Kawazoe, Y. Phys. Chem. Chem. Phys. 2010, 12, 9839. (b) Sun, Z.; Huang, K.-W.; Wu, J. J. Am. Chem. Soc. 2011, 133, 11896. (c) Gao, X.; Hodgson, J. L.; Jiang, D-e.; Zhang, S. B.; Nagase, S.; Miller, G. P.; Chen, Z. Org. Lett. 2011, 13, 3316. (13) (a) Shao, Y.; Head-Gordon, M.; Krylov, A. I. J. Chem. Phys. 2003, 118, 4807. (b) Wang, F.; Ziegler, T. J. Chem. Phys. 2004, 121, 12191. (c) Golubeva, A. A.; Nemukhin, A. V.; Klippenstein, S. J.; Harding, L. B.; Krylov, A. I. J. Phys. Chem. A 2007, 111, 13264. (14) (a) Saito, T.; Nishihara, S.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Yamada, S.; Isobe, H.; Okumura, M.; Yamaguchi, K. Theor. Chem. Acc. 2011, 130, 739. (b) Saito, T.; Nishihara, S.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Yamada, S.; Isobe, H.; Okumura, M.; Yamaguchi, K. Theor. Chem. Acc. 2011, 130, 749. (15) Frisch, M. J.; et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (16) Grimme, S. J. Chem. Phys. 2003, 118, 9095. (17) Jung, Y.; Lochan, R. C.; Dutoi, A. D.; Head-Gordon, M. J. Chem. Phys. 2004, 121, 9793. (18) For atomic and diatomic singlet species the use of unrestricted α and β spin orbitals results in maximum spin polarization that mimics the spatial occupation of a multireference type calculation. Calculation of the unrestricted energy for methylene and nitrenium was carried out by requiring α electron occupation in the 3a1 orbital and β electron occupation of the 1b1 orbital. Importantly, this results in a different
correlation energies (SP-HF-SCS-MP2 or SP-HF-SOSMP2). For ST gap energies, density functionals generally give poor results without spin-projection. From currently available functionals, only ωB97X-D and M06-2X functionals give reasonable accuracy for ST gap energies with spin-projection. The newly introduced DR-1 functional gives highly accurate spin-projected ST gap energies. However, closed-shell singlet energies are underestimated. For the methane C−H bond most density functionals give accurate dissociation curves. In contrast, for H−F dissociation M06-2X, ωB97X-D, and DR-1 give unreasonably high errors greater than 50 kcal/mol. H−F dissociation can be modeled accurately using the SOS-MP2 method and is best if used without spinprojection.
ASSOCIATED CONTENT
S Supporting Information *
Singlet and triplet bond lengths and angles for molecules and full reference 15. This material is available free of charge via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank the Fulton Supercomputing Lab for computational support. Acknowledgment is made to the Donors of the American Chemical Society Petroleum Research Fund for support of this research (to D.H.E; 51081-DNI3). T.C.C. thanks the BYU Department of Chemistry and Biochemistry for support. D.H.E. also thanks Erin Johnson (UC Merced), Weitao Yang (Duke), and Ernest Davidson (U. of Washington) for useful conversations.
■
REFERENCES
(1) (a) Salem, L.; Rowland, C. Angew. Chem., Int. Ed. Engl. 1972, 11, 92. (b) In accord with Salem we use the word singlet diradical to describe a species with opposite spin electrons occupying degenerate or nearly degenerate orbitals. (2) Gräfenstein, J.; Cremer, D. Phys. Chem. Chem. Phys. 2000, 2, 2091. (3) (a) Pople, J. A.; Nesbet, R. K. J. Chem. Phys. 1954, 22, 571. (b) Wang, J.; Becke, A. D.; Smith, V. H., Jr. J. Chem. Phys. 1995, 102, 3477. (c) Pople, J. A.; Gill, P. M. W.; Handy, N. C. Int. J. Quantum Chem. 1995, 56, 303. (d) Grafenstein, J.; Hjerpe, A. M.; Kraka, E.; Cremer, D. J. Phys. Chem. A 2000, 104, 1748. (e) Cremer, D. Mol. Phys. 2001, 99, 1899−1940. (f) Cohen, A. J.; Tozer, D. J.; Handy, N. C. J. Chem. Phys. 2007, 126, 214104. (g) Menon, A. S.; Radom, L. J. Phys. Chem. A 2008, 112, 13225. (h) Menon, A. S.; Wood, G. P. F; Moran, D.; Radom, L. J. Phys. Chem. A 2007, 111, 13638. (i) Menon, A. S.; Radom, L. J. Phys. Chem. A 2008, 112, 13225. (j) Graham, D. C.; Menon, A. S.; Goerigk, L.; Grimme, S.; Radom, L. J. Phys. Chem. A 2009, 113, 9861. (4) Slipchenko, L. V.; Krylov, A. I. J. Chem. Phys. 2002, 117, 4694. (5) Ess, D. H.; Johnson, E. R.; Hu, X.; Yang, W. J. Phys. Chem. A 2011, 115, 76. (6) (a) Yamaguchi, K.; Yoshioka, Y.; Fueno, T. Chem. Phys. Lett. 1977, 46, 360. (b) Yamaguchi, K.; Yoshioka, Y.; Takatsuka, K.; Fueno, 4928
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
The Journal of Physical Chemistry A
Article
energy solution compared to when both α and β electrons occupy polarized 3a1 type orbitals. (19) (a) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215. (b) Zhao, Y.; Truhlar, D. G. Acc. Chem. Res. 2008, 41, 157. (20) T. Yanai, T.; Tew, D.; Handy, N. Chem. Phys. Lett. 2004, 393, 51. (21) (a) Chai, J.-D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2008, 10, 6615. (b) Chai, J.-D.; Head-Gordon, M. J. Chem. Phys. 2008, 128, 084106. (22) Boese, A. D.; Martin, J. M. L. J. Chem. Phys. 2004, 121, 3405. (23) (a) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (b) Becke, A. D. Phys. Rev. A 1988, 38, 3098. (c) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (d) Less, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (e) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. J. Phys. Chem. 1994, 98, 11623. (24) Xu, X.; Goddard, W. A., III Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 2673. (25) Adamo, C.; Barone, V. J. Chem. Phys. 1998, 108, 664. (26) Grimme, S. J. Chem. Phys. 2006, 124, 034108. (27) See Supporting Information for results for LC-ωPBE, ωB97D, M06L, M06HF, and B97D functionals. (28) Lynch, B. J.; Fast, P. L.; Harris, M.; Truhlar, D. G. J. Phys Chem. A 2000, 104, 4811. (29) The MPW1K method gives ME and MUE values of −6.1 and 6.2 kcal/mol, respectively, for singlet−triplet gaps. (30) For a discussion on using unrestricted calculations for the 1A1 state of methylene see: Kitagawa, Y.; Saito, T.; Nakanishi, Y.; Kataoka, Y.; Matsui, T.; Kawakami, T.; Okumura, M.; Yamaguchi, K. J. Phys. Chem. A 2009, 113, 15041. (31) Jensen, P.; Bunker, P. R. J. Chem. Phys. 1988, 89, 1327. (32) (a) Berkowitz, J.; Greene, J. P.; Cho, H.; Ruscic, B. J. J. Chem. Phys. 1987, 86, 1235. (b) Escribano, R.; Campargue, A. J. Chem. Phys. 1998, 108, 6249. (33) Berkowitz, J.; Cho, H. J. Chem. Phys. 1989, 90, 1. (34) (a) Becke, A. D. J. Chem. Phys. 2002, 117, 6935. (b) Becke, A. D. J. Chem. Phys. 2003, 119, 2972. (c) Becke, A. D. J. Chem. Phys. 2005, 122, 64101. (d) Proynov, E.; Shao, Y.; Kong, J. Chem. Phys. Lett. 2010, 493, 381. (35) (a) Becke, A. D.; Johnson, E. R. J. Chem. Phys. 2007, 127, 124108. (b) Johnson, E. R.; Becke, A. D. J. Chem. Phys. 2008, 128, 124105. (36) (a) Toulouse, J.; Gori-Giorgi, P.; Savin, A. Theor. Chem. Acc. 2005, 114, 305. (b) Weimer, M.; Sala, F. D.; Gorling, A. J. Chem. Phys. 2008, 128, 144109. (37) (a) Krylov, A. I. Chem. Phys. Lett. 2001, 350, 522. (b) Krylov, A. I. Chem. Phys. Lett. 2001, 338, 375. (c) Krylov, A. I.; Sherrill, C. D. J. Chem. Phys. 2002, 116, 3194. (d) Abrams, M. L.; Sherrill, C. D. J. Phys. Chem. A 2003, 107, 5611. (e) Sears, J. S.; Sherrill, C. D.; Krylov, A. I. J. Chem. Phys. 2003, 118, 9084. (f) Casanova, D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2009, 11, 9779. (38) Dutta, A.; Sherrill, C. D. J. Chem. Phys. 2003, 118, 1610. (39) Triplet energies used to perform spin-projection were evaluated using the singlet geometries.
4929
dx.doi.org/10.1021/jp300633j | J. Phys. Chem. A 2012, 116, 4922−4929
