Predicting Singlet–Triplet Energy Splittings with Projected Hartree

Article pubs.acs.org/JPCA

Predicting Singlet−Triplet Energy Splittings with Projected Hartree− Fock Methods Pablo Rivero,*,† Carlos A. Jiménez-Hoyos,† and Gustavo E. Scuseria†,‡ †

Department of Chemistry, Rice University, Houston, Texas 77005-1892, United States Department of Physics and Astronomy, Rice University, Houston, Texas 77005-1892, United States

‡

ABSTRACT: Hartree−Fock (HF) and density functional theory (DFT) methods are known for having problems in predicting singlet−triplet energy splittings when the system displays significant diradical character. Multireference methods are traditionally advocated to deal with the spin-contamination problem inherent in broken-symmetry mean-field methods. In the present work, spin-contamination is rigorously eliminated by means of a symmetry projection approach, carried out in a variation-after-projection fashion, recently implemented in our research group. We here explore the performance of a variety of projected Hartree−Fock (PHF) approaches (SUHF, KSUHF, SGHF, and KSGHF) in predicting singlet−triplet energy gaps in a broad set of diradical systems: small diatomic molecules, carbenes and silenes, and a few larger molecules (trimethylenemethane and benzyne isomers). For most of these systems, accurate experimental data is available in the literature. Additionally, we assess the quality of the geometrical parameters obtained in SUHF-based optimizations for some of the systems considered. Our results indicate that PHF methods yield highquality multireference wave functions, providing a good description of the ground state potential surface as well as an accurate singlet−triplet splitting gap, all within a modest mean-field computational cost.

1. INTRODUCTION A plethora of quantum mechanical systems can be characterized by their singlet−triplet splitting ΔEST,1 which is the difference between their lowest-lying singlet and triplet states. This energy gap usually provides insight into the relative stability and reactivity of organic systems.2,3 They can further be used to characterize optoelectronic devices4 or as the fundamental parameter in semiconductor spin-qubits.5 In chemistry, diradicals are molecular species with two electrons occupying two degenerate (or nearly degenerate) molecular orbitals. The two unpaired electrons may be aligned ferro- or antiferromagnetically depending on the magnitude of the coupling constant. Four possible states are available for a diradical system: a singlet and a set of degenerate triplet states. The accurate theoretical description of diradical systems necessitates a two-configuration wave function for the singlet and the m = 0 triplet states. The m = ±1 states, however, may be described by means of a single configuration. Because of this, multireference methods must be used to provide an accurate theoretical estimate of the singlet−triplet gaps. In this sense, the simplest approach to accurately handle the static correlation among the set of degenerate orbitals is the CASSCF(2,2) wave function. In some cases, nevertheless, accounting for dynamical correlation is also of paramount importance in order to obtain an accurate singlet−triplet splitting. Evidently, the exact answer would be obtained by performing an exact diagonalization of the Hamiltonian among all the possible configurations spanned © 2013 American Chemical Society

by the chosen single-particle basis. Such FCI approach is, however, impractical for all but the smallest of the systems due to the combinatorial scaling of the method. We note that symmetry-breaking in a single-determinantal unrestricted framework allows one to account for a fraction of the strong correlation that otherwise necessitates a multireference treatment. In these nearly degenerate cases, allowing the up- and down-electrons to have different spatial distributions, that is, partial orbital localization, lowers the energy by minimizing the electron−electron repulsive interactions and reopens the orbital gap. The use of symmetry-broken solutions (both Hartree− Fock and Kohn−Sham) to describe open-shell singlet states has become common practice among quantum chemists. Several authors have assessed the ability of many different multireference ab initio approaches to accurately predict singlet−triplet gaps. The CASSCF approach, which correctly accounts for strong correlations in the active space but lacks the important dynamical correlations outside this space, has been widely used.6−8 Other work has been carried out using methods that build correlations outside the active space like MRCISD,9,10 CASPT2,11,12 MCQDPT2,13 DDCI,14 and variants of multireference coupled-cluster (MRCC),15,16 among others. In general, these traditional multireference Received: June 10, 2013 Revised: July 11, 2013 Published: July 18, 2013 8073

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

We here consider four different PHF variants. In SUHF and KSUHF, we use a UHF-type underlying determinant (breaking S2); in the former, a full spin projection is carried out, while in the latter an additional complex-conjugation projection is performed.33,34 Additionally, we consider SGHF and KSGHF, in which a generalized HF 40 (GHF)-type underlying determinant (breaking S2 and Sz) is used. To approximate the true ΔEST in the case of spin contaminated wave functions (UHF or DFT), a semiempirical correction equation (eq 1) is commonly used. This equation was initially proposed by Noodleman41,42 to describe magnetic coupling constants using broken-symmetry Slater determinants.

approaches have a steep computational cost and become impractical to study large systems. During the last few decades, density functional theory (DFT)17,18 has become the method of choice in quantum chemical studies because of its comparatively low computational cost. Indeed, different DFT approaches have been used by several research groups to compute singlet−triplet splittings.9,19 However, approximate functionals do not always provide reliable results. Additionally, DFT is typically used in a symmetry-broken (unrestricted) formulation, which yields an admixture of multiple spin states. We note that the evaluation of ⟨S2⟩ using the underlying noninteracting reference system is not rigorously justified.20 Attempts to merge multireference methods with DFT approximations, such as the restricted open-shell DFT, the spin-restricted ensemble-referenced Kohn−Sham21 method (REKS) ROSS-DFT,22 the completeactive-space DFT, CAS-DFT23 method, and our own recent efforts,24 typically face the challenge of avoiding a doublecounting of electron correlations. Recent work has explored the consequences of removing spin contamination from the underlying determinants in the prediction of singlet−triplet energy gaps. Tsuchimochi et al.25 carried out a constrained variation, limiting the amount of spin contamination introduced in the underlying broken symmetry wave function. This work showed that singlet−triplet energy gaps are largely improved. Moreover, a full projection from such a wave function is more easily accomplished. On other work, Saito et al.26 performed singlet−triplet calculations in diradical species using the approximate spin-projection method described by Yamaguchi;27,28 their results show that singlet− triplet splittings are significantly improved after the spin contamination is approximately removed. The approximate projection scheme has attracted increasing interest from the quantum chemistry community (see, e.g., refs 29−32.). The common theme in all these works is the importance of removing spin contamination from symmetry-broken calculations. In this work, we investigate the merits of our recently implemented projected Hatree−Fock methods33 (PHF), where a full spin projection is carried out in a variation-afterprojection (VAP) approach, that is, where the optimization of the underlying wave function (a single Slater determinant) is done in the presence of the projection operators. This is better and fundamentally different from the projection-after-variation (PAV) approach predominantly used in previous work. We have demonstrated that VAP symmetry-projected wave functions can account for a large fraction of both static and dynamic correlations. Additional merits are its black-box character, as well as its mean-field computational cost that permits the study of large systems. Evidence of its potential has been shown in recent applications to copper oxide cores,34 oneand two-dimensional Hubbard models,35,36 and polycyclic aromatic hydrocarbons.37 PHF has, as any other approximation, some disadvantages; perhaps the main one is its lack of size consistency. Although PHF has a long history in quantum chemistry,38,39 the full potential of the method was never explored because it was not efficiently formulated as an orbital optimization problem. Fundamentally, our approach avoids using Lowdin’s many-body spin projection operators in favor of one-body spin rotation operators that make the wave function rotationally invariant is spin space and exploits group theory techniques that are rooted on generalized coherent states of Lie algebras.

N ΔEST =2

ES − E T 2

⟨S ⟩S − ⟨S2⟩T

(1)

where ES and ET are the energies of the spin contaminated singlet and triplet determinants and ⟨S2⟩S and ⟨S2⟩T correspond to expectation values of S2. For comparison purposes, we here include results obtained with the above approximate formula. The goal of this article is to assess PHF variants for describing ground state geometries (only with SUHF) and vertical and adiabatic singlet−triplet splittings. The spinprojected wave function handles well the static correlation in diradical systems and adds a fraction of the dynamical correlations. Complex-conjugation projection appears to account mostly for residual dynamic correlations. We are also interested in comparing the performance of the VAP and the PAV approaches for singlet−triplet splittings. We consider four small diatomic molecules, eight carbenes, and two silenes, as well as a small set of larger molecules (trimethylenemethane and the three benzyne isomers). The availability of experimental splittings for most of these systems allow us to systematically assess the ability of PHF for describing this important electronic property. The article is organized as follows: in section 2, computational details are described; in section 3, the geometries and singlet−triplet energy gaps are reported; finally, in section 4, we present our conclusions.

2. COMPUTATIONAL DETAILS For convenience, we separate the studied systems into three groups: 4 small diatomic molecules, 10 carbenes and silenes, and 4 larger molecules. For diatomic systems, we have used the experimental geometries. In other cases, we use geometries optimized with PBE0,43,44 with SUHF, or as described below. The SUHF geometry optimizations have been carried out via numerical differentiation of the energy functional. The first group consists of four small diatomic molecules: NF, NH, O2, and OH+. In this case, we compute adiabatic and vertical singlet−triplet splittings using the experimental geometries. In addition, we have computed adiabatic and vertical splittings for PBE0 and SUHF using optimized geometries with these methods. All calculations have been carried out using the cc-pVTZ basis set. The next group consists of eight carbenes (CFBr, CHBr, CBr2, CCl2, CHCl, CF2, CH2, and HCF) and 2 silenes (SiH2 and SiF2). For these systems, we first compare the available experimental geometrical parameters with those predicted by PBE0 and SUHF using the cc-pVDZ basis. Given the good agreement among all sets of geometrical parameters, we have used SUHF geometries to evaluate singlet−triplet splittings of all carbenes. We report adiabatic splittings computed with the 8074

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

cc-pVDZ and cc-pVTZ basis for all methods, as well as the ccpVQZ basis for SUHF and KSUHF. Finally, we have considered a small set of larger molecules (trimethylenemethane and o,- m-, and p-benzynes). In these cases, we use geometries extracted from ref 25 and report adiabatic splittings predicted using SUHF, KSUHF, and SGHF with the cc-pVDZ and cc-pVTZ basis sets. We also include KSGHF splittings evaluated with the cc-pVDZ basis set. All calculations have been carried out using a locally modified version of the Gaussian09 suite of programs.45 We have used enough grid points in the spin-projected PHF methods to obtain ⟨S2⟩ correct to a precision of 10−5 or better, fully eliminating spin contamination. Our initial guess for SUHF calculations is either a true UHF minimum (for systems for which RHF is unstable) or a UHF-like determinant prepared from a restricted solution using standard orbital mixing schemes. In order to prepare an initial guess for GHF-based calculations, we use a standard Gaussian procedure (guess = mix) for breaking the Sz symmetry of the converged SUHF solution. We break complex conjugation symmetry by mixing the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) with complex coefficients. We have carried out CASSCF(2,2) calculations for all systems under study using the Gaussian09 suite of programs. Evidently, CASSCF(2,2) for a triplet state reduces to an ROHF triplet. However, for the biradical species, we have used UHF natural orbitals as an initial guess for the CAS calculations, as suggested by Pulay and Hamilton.46 For TMM, we have used CASSCF(4,4) instead of CASSCF(2,2) in order to correlate the four p-orbitals associated to the carbon centers and perpendicular to the plane of the molecule. In this work, the singlet−triplet splitting is defined as the energy of the singlet minus the energy of the triplet state. Therefore, a system with a singlet ground state will exhibit a negative energy splitting.

Table 1. Optimized Bond Lengths (in Å) for Diatomic Molecules; Experimental Values Were Taken from Ref 47 molecule O2

NF

OH+

NH

method

singlet

triplet

SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl

1.222 1.201 1.216 1.276 1.301 1.308 1.008 1.033 1.029 1.019 1.039 1.034

1.164 1.194 1.208 1.301 1.304 1.317 1.015 1.033 1.029 1.026 1.039 1.036

Table 2. ⟨S2⟩ Values for UHF and PBE0 at the Experimental Geometries Using the cc-pVTZ Basis Set; for All PHF Methods Considered in This Work, ⟨S2⟩ Is Correctly Predicted to Be s(s + 1) method

ms

O2

NF

OH+

NH

UHF

0 1 0 1

1.0255 2.0433 1.0054 2.0092

1.0130 2.0210 1.0052 2.0074

1.0079 2.0139 1.0035 2.0054

1.0106 2.0153 1.0035 2.0055

PBE0

considered. They all predict ⟨S2⟩ ≈ 1, which corresponds approximately to an equal admixture of pure singlet and triplet states. However, the triplet state does not represent a significant problem for the single-determinantal methods, providing in this case a value of ⟨S2⟩ close to 2. All the PHF methods used in this work involve the total spin projection; therefore, there is no spin contamination in any of the PHF wave functions. The adiabatic and vertical singlet−triplet splittings predicted (using the cc-pVTZ basis) with a variety of methods (PHF, UHF, and PBE0) are shown in Table 3. Here, SUHF(PAV) refers to the projection-after-variation SUHF method, where a single-shot symmetry projection is performed on the optimized UHF states. We also include for comparison CASSCF(2,2) and experimental values taken from the literature. We have evaluated the mean absolute error (MAE) with respect to experimental values for both the vertical and adiabatic energies. All methods correctly predict the triplet to be the ground state. The lowest mean absolute error (MAE) is obtained with the SUHF method, yielding results even better than those from CASSCF(2,2). We have verified that our results are well converged with respect to the basis set size for PHF methods (results with the cc-pVQZ basis were carried out, but they are not shown). PBE0 and UHF fail to provide a good description of the singlet−triplet gaps, yielding MAE values of 19.4 and 16.2 kcal/mol, respectively. This large error is mostly due to the massive spin contamination that these methods possess for the singlet states. We have also included results obtained with the SUHF(PAV) strategy. The low MAE obtained with such a strategy suggests that the restoration of symmetry may indeed provide accurate singlet−triplet splittings even when the full variational optimization has not been carried out. Nevertheless, the variation-after-projection (VAP) strategy should be preferred in all cases, as a PAV approach may lead to artificial results. Singlet−triplet splittings using the empirical formula of

3. RESULTS AND DISCUSSION In this section, we report the performance of PHF methods in the prediction of singlet−triplet splittings. Results obtained are compared with the well-known hybrid density functional PBE0 and the CASSCF(2,2) multireference method as well as some experimental values. We divide our discussion in three independent subsections according to the groups described above. 3.1. Diatomic Molecules. We start by considering a small set of diatomic molecules: NF, NH, O2, and OH+. The optimized bonding distances for these systems computed at the SUHF and PBE0/cc-pVDZ level are summarized in Table 1. Results show the optimized bond lengths for the lowest singlet and triplet states; in all cases, the triplet state corresponds to the ground state. We provide experimental values47 extracted from the literature as a reference. All systems but O2 exhibit larger interatomic distances for the triplet ground state. We observe generally good agreement between the PBE0 and SUHF optimized geometries, as well as with experimental vaues. The few outliers (triplet O2 and singlet NF) may require a more complete treatment of electron correlations in a larger basis set. We will employ SUHF geometries for systems with few geometrical degrees of freedom for the rest of this work. In Table 2 the ⟨S2⟩ value for PBE0 and UHF methods is shown. We point out that both methods exhibit a large spin contamination in the singlet states of all the diatomic systems 8075

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

CHCl, CCl2, CHBr, SiH2, and SiF2. In some cases, high-level computational predictions of singlet−triplet gaps are in disagreement with the available experimental data. In the case of CH2, FCI/TZVP predicts a gap of 11.97 kcal/mol,49 whereas the experimental value is slightly smaller, at 9 kcal/mol. This is most likely a basis set incompleteness error. However, for CBr2, there is a clear disagreement between theory and experiment, with CCSD(T)/cc-pVQZ predicting a −16.26 kcal/mol gap50 as opposed to the experimental −2 ± 3 kcal/mol value.51 Lastly, no experimental result is found for CFBr, but MRCI +DAV/cc-pVTZ and QCISD(T)/6-311+G predict −32.6450 and −31.5952 kcal/mol, respectively. In Table 4, we report the optimized and experimental geometries calculated for the systems under study. As expected, the pairwise repulsion between the unpaired electrons in triplet states lead to larger bond angles than in singlet states. We observe that SUHF optimized geometries are in good agreement with PBE0 and existing experimental data. We present in Table 5 the spin contamination obtained with UHF and PBE0 for the carbenes and silenes under consideration. UHF shows large spin contaminations for all singlet states (ms = 0) but for SiF2. PBE0, however, only shows large spin contaminations for CH2 and CHBr. In light of these, we anticipate much better singlet−triplet gaps with PBE0 than with UHF. Table 6 shows the singlet−triplet gaps computed with a variety of methods. We note that UHF fails to yield accurate gaps, as discussed above. In some cases, it is not even capable of predicting the right spin ground state (CHBr, CBr2, and HCF). Interestingly, PBE0 also fails to predict the correct ground state for CHCl. All other methods (PHF and CASSCF(2,2)) predict the correct ground state: a triplet for CH2 and a singlet for the rest of carbenes and silenes. Nevertheless, the predicted gaps vary widely among the different methods. PBE0 provides a MAE of 6.39 kcal/mol, much better than UHF (18.04 kcal/ mol) using the cc-pVTZ basis, due in part to the low spin contamination exhibited. Contrary to what was observed for the diatomics, the approach using Noodleman’s formula (eq 1) fails to improve the singlet−triplet splittings obtained with UHF and PBE0 as the MAE values are larger than without the semiempirical correction. The set of PHF methods tested, which do not have any spin contamination, provide singlet− triplet splittings in good agreement with experimental values. Judging by the mean absolute error, PHF methods offer a better description than PBE0 or CASSCF(2,2) for these diradical-like systems. We note that at least a triple-ζ quality basis set is necessary to converge singlet−triplet gaps with respect to basis size. Although all PHF methods yield an accurate description of singlet−triplet gaps, KSGHF provides the best agreement with experimental values for the carbenes and silenes under consideration. In the case of CF2 or HCF, the error with respect to experiment is below 5%. By breaking and restoring the Sz symmetry of the reference determinant, which results in a noncollinear spin representation, these methods have the flexibility of incorporating more dynamical correlation. This appears to be important to yield a more accurate description of singlet−triplet gaps. In our current implementation, which is far from ideal, SGHF and KSGHF are computationally more expensive than SUHF, a shortcoming that will be addressed in future work. 3.3. TMM and p-, o-, and m-Benzynes. Trimethylenemethane (TMM) is an example of a non-Kekulé diradical

Table 3. Singlet−Triplet Splittings (in kcal/mol) for Diatomic Molecules Evaluated Using the cc-pVTZ Basis Seta method UHF UHF† SUHF(PAV) SUHF SUHF* KSUHF SGHF KSGHF PBE0 PBE0† PBE0* CASSCF(2,2)c exptld

v/ab

O2

NF

OH+

NH

MAE

v a v a v a v a v a v a v a v a v a v a v a a

16.53 17.02 33.06 34.04 28.51 28.89 22.75 22.82 25.26 25.12 23.86 23.69 7.41 7.49 17.05 16.96 10.29 10.02 20.48 20.04 10.25 10.30 21.02 22.6

19.89 19.65 39.78 39.30 28.86 29.38 32.61 32.26 31.96 32.07 31.10 31.01 26.31 26.34 30.18 30.01 12.30 12.38 24.60 24.72 12.25 12.26 40.65 34.3

25.86

19.42 19.42 38.84 38.84 39.37 39.33 33.69 33.65 33.62 33.59 31.59 31.59 33.33 33.35 32.14 32.15 13.87 13.88 27.74 27.76 13.87 13.87 40.95 39.0

16.18 22.58 4.25 5.43 3.88 3.83 2.96 3.08 3.79 3.74 4.74 4.72 8.52 8.49 5.40 5.47 22.72 22.76 8.86 18.47 22.74 19.44 3.12

51.72 46.71 45.80 45.71 45.69 43.42 55.72 45.42 19.08 38.16 19.06 53.57 50.5

a

We have used experimental geometries for all methods except SUHF* and PBE0* for which self-consistent geometries have been used. UHF† and PBE0† refer to the splitting energies using eq 1. b Vertical (v) or adiabatic (a) splittings. cResults taken from ref 25. UMP2/TZ2P geometries have been used. dExperimental values are taken from ref 47.

eq 1 for UHF and PBE0 provide a notable improvement in the MAE values for these diatomic systems. It is interesting to note that results seem to deteriorate when more elaborate PHF methods are used when compared to SUHF. Our computed MAE values should only be taken as a reference; comparison with experimental values is difficult due to zero-point and thermal effects in the latter case. Because we are only computing electronic energies, results should be compared with FCI results, which are not available to us. 3.2. Carbenes and Silenes. We now turn our attention to the group consisting of eight carbenes (CFBr, CHBr, CBr2, CCl2, CHCl, CF2, CH2, and HCF) and 2 silenes (SiH2 and SiF2). Methylene (CH2), which has been intensely studied both theoretically and experimentally, exhibits a triplet ground state and a singlet−triplet splitting of 9.0 kcal/mol.48 With the exception of methylene, the rest of carbenes and silenes that we have studied have a singlet ground state, with a broad range of singlet−triplet gaps. The singlet state in carbenes and silenes can have either closed shell or diradical character if both electrons occupy different orbitals. The incomplete octet shell around the carbon center confers the carbenes and silenes with reactive properties that make them very useful for organic synthesis and organometallic chemistry. Singlet−triplet energy splittings for carbene compounds have been studied by several experimental and theoretical approaches. Experimental singlet−triplet gaps are available for almost all carbenes and silenes studied in this article. Well known experimental values were determined for HCF, CF2, 8076

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

Table 4. Optimized Bond Lengths (in Å) and Angles (in Degrees) for Carbenes and Silenes in Calculations with the cc-pVDZ Basisa singlet

triplet

molecule

method

R−X

R−Y

∠

R−X

R−Y

∠

CFBr

SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl SUHF PBE0 exptl

1.288 1.292

1.939 1.942

106.4 106.4

1.302 1.311

1.852 1.851

123.6 124.1

1.114 1.125 1.110 1.891 1.894 1.865 1.730 1.732 1.711 1.711 1.707

1.881 1.862 1.857 1.891 1.894 1.865 1.730 1.732 1.711 1.114 1.125

100.7 100.3 101.0 111.4 110.0 110.7 110.0 109.1 109.2 102.0 101.2

1.086 1.096

1.834 1.815

126.6 123.7

1.841 1.833

1.841 1.833

130.1 129.4

1.689 1.683

1.689 1.683

127.5 127.7

1.678 1.670

1.085 1.095

125.7 126.1

1.291 1.306 1.297 1.107 1.126 1.107 1.107 1.127 1.138 1.522 1.541 1.514 1.634 1.652 1.590

1.291 1.306 1.297 1.107 1.126 1.107 1.302 1.301 1.305 1.522 1.541 1.514 1.634 1.651 1.590

104.8 104.3 104.9 103.1 100.1 102.4 102.9 102.2 104.1 94.0 91.2 92.0 98.7 99.7 100.8

1.303 1.316

1.303 1.316

118.7 119.6

1.082 1.090 1.085 1.299 1.305

1.082 1.090 1.085 1.073 1.089

134.7 134.2 135.5 122.2 121.9

1.481 1.499

1.481 1.499

118.6 119.0

1.633 1.651

1.633 1.651

112.9 115.1

CHBr

CBr2

CCl2

CHCl

CF2

CH2

CHF

SiH2

SiF2

a Here, R−X corresponds to the bond distance between the carbon (or silicon) center and the lightest of the bonded atoms, while R−Y corresponds to the other bond length. Experimental data is shown in cases where it is available.

Table 5. ⟨S2⟩ Values for UHF and PBE0 at the Experimental Geometries Using the cc-pVTZ Basis Set; for All PHF Methods Considered in This Work, ⟨S2⟩ Is Correctly Predicted to Be s(s + 1) method

ms

CFBr

CHBr

CBr2

CCl2

CHCl

CF2

CH2

HCF

SiH2

SiF2

UHF

0 1 0 1

0.355 2.015 0.000 2.006

0.611 2.022 0.415 2.009

0.585 2.027 0.000 2.088

0.507 2.020 0.000 2.009

0.599 2.017 0.000 2.008

0.104 2.007 0.000 2.003

0.713 2.015 0.000 2.066

0.495 2.009 0.330 2.044

0.408 2.005 0.000 2.003

0.000 2.003 0.000 2.001

PBE0

molecule. In a D3h configuration, the two unpaired electrons distribute over two degenerate π orbitals to yield a 3A′2 ground state. However, the system undergoes Jahn−Teller distortion in the singlet excited state to a C2v configuration. We note that there are two competing low-energy singlet states in the TMM potential energy surface. As suggested by the experimental results of Wenthold et al.,57 we have used the A1 planar state, which has a slightly higher energy at the CASPT2(10,10)/ccpVTZ58 level than the twisted state, in our calculations. We use the geometries reported by Cramer and Smith for the singlet (A1) and triplet (A2) states of TMM obtained at the CASSCF(10,10)/cc-pVDZ level of theory.58 Benzynes are highly reactive species derived from an aromatic ring by removal of two (ortho, meta, or para) substituents. In this work, we concentrate on the simplest benzynes derived from the benzene molecule. These compounds are known to be

involved as intermediates in numerous important reactions and are of great interest for organic chemists and biochemists. The singlet−triplet splittings of benzyne molecules have been studied both experimentally and theoretically. All benzyne isomers exhibit a singlet ground state transforming as the totally symmetric irreducible representation. The o- and m-benzynes display C2v symmetry, whereas p-benzyne shows D2h symmetry. We follow ref 25 and use UB3LYP/6-31G(d) optimized geometries for the singlet and triplet states of all considered benzynes. We note, in particular, that there are two competing structures in the singlet state of m-benzyne. In ref 25, a (highly strained) bicyclic closed-shell structure was used, and the values of the singlet−triplet splittings computed with it were in large discrepancy with experimental values. In this article, we consider the alternative minimum, best described as a monocyclic biradical species converged with UB3LYP. We 8077

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

Table 6. Singlet−Triplet Gaps (in kcal/mol) for Carbenes and Silenes Using PHF Methods and SUHF Optimized Geometriesa molecule CFBr

CHBr

CBr2

CCl2

CHCl

CF2

CH2

HCF

SiH2

SiF2

MAE

basis

UHF

UHF†

SUHF

KSUHF

SGHF

KSGHF

PBE0

PBE0†

CAS(2,2)

D T Q D T Q D T Q D T Q D T Q D T Q D T Q D T Q D T Q D T Q D T Q

−10.70 −10.31

−12.89

−32.19 −34.19

−27.18 −26.38

−27.18

−24.77

13.90

−2.31 −4.25

−4.07 −8.29

−1.02 −0.99

−1.25

0.41

3.62 4.02

4.63

−10.47 −12.26

−14.10 −19.27

−10.92 −10.20

−10.92

−7.27

−16.26c

−0.43 −0.05

−0.57

−15.29 −17.56

−18.91 −22.64

−14.52 −13.86

−14.52

−12.35

−20.5d

8.72 8.58

12.47

−3.89 −5.87

−5.34 −7.10

0.89 0.94

0.89

−0.91

−6.4e

−32.31 −32.72

−33.92

−48.75 −50.14

−54.19 −55.70

−48.96 −48.92

−48.96

−47.53

−56.7e

28.6 28.32

−44.04

12.2 10.97

9.95 6.74

18.11 17.73

18.11

11.96

9.0e

0.96 1.23

1.28

−10.72 −11.68

−14.87 −15.56

−7.07 −8.52

−8.27

−8.20

−14.7e

−6.79 −7.10

−8.49

−17.87 −20.49

−20.12 −24.15

−14.80 −14.97

−14.80

−18.02

−21,−18f

−54.12 −51.88

−54.12

−69.09 −65.37

−72.03 −71.06

−69.29 −68.28

−69.29

−54.99

−75.2,−77.2g

17.80 18.04

20.38

−29.76 −30.52 −30.53 0.18 −1.44 −1.78 −8.69 −10.35 −9.040 −13.59 −15.71 −16.00 −1.51 −2.70 −58.90 −58.36 −58.62 −58.90 15.89 13.94 13.49 −12.13 −12.11 3.17 −17.23 −20.54 −20.35 −78.46 −76.77 −76.67 4.17 2.98 2.84

−35.77 −38.11

9.79 9.95

−32.09 −32.13 −32.31 0.61 −1.63 −1.91 −6.95 −8.67 −9.08 −12.60 −14.88 −15.22 −2.06 −3.26 −50.47 −48.75 −50.38 −50.47 12.53 15.45 15.06 −9.65 −10.53 3.70 −17.13 −20.43 −20.35 −52.9 −79.38 −79.40 5.17 4.24 4.03

4.37 3.54

1.59 2.43

6.42 6.39

6.28

6.82

exptl −32.6

−5.7b

a Here, D, T, and Q in the basis column refer to the cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets, respectively.53 UHF† and PBE0† refer to the splitting energies using eq 1. bReference 53. cReference 50. CCSD(T)/cc-pVQZ level. dReference 54. eReference 55. fReference 48. gReference 56.

KSUHF, however, yield splittings that are far from the experimental values for all benzynes, although the right ground state is predicted. We speculate that these methods incorporate a fraction of dynamical correlation in an unbalanced way such that the splittings obtained are worse than those from small CAS calculations. However, SGHF and KSGHF yield consistently good splittings for all the systems in this test set. We note that the use of Noodleman’s formula (eq 1) with UHF or PBE does not provide an improvement when comparing to the experimental values. We close this section by summarizing all computed singlet− triplet mean absolute errors in Figure 1. Here, we have used ccpVTZ calculations for all systems with the following exceptions: in carbenes and silenes, we use cc-pVDZ for CASSCF, and in benzynes and TMM, we use cc-pVDZ for both CASSCF and KSGHF. For the diatomic molecules, we have used the vertical splittings in preparing the figure. Several interesting features are observed: (i) UHF is, by far, the most inaccurate of the methods considered. A consistently large MAE is obtained, which is mostly due to spin contamination exhibited in the singlet states; (ii) PBE0 yields low MAE values in some test sets, but its performance is poor for the diatomic molecules, and we thus deem it unreliable; (iii) SUHF and KSUHF provide very accurate singlet−triplet splittings (considering

note that the latter is slightly higher in energy at the B3LYP level, yet it gives splittings in better agreement with experiment for methods such as CASSCF. We point the interested reader to the work of Winkler59 and Sander on the predicted DFT structures for m-benzyne. Spin contaminations computed for these systems using UHF and PBE0 are presented in Table 7. Table 7. ⟨S2⟩ Values for UHF and PBE0 Using the cc-pVDZ Basis Set; for All PHF Methods Considered in This Work, ⟨S2⟩ Is Correctly Predicted to Be s(s + 1) method

ms

TMM

p-benzyne

o-benzyne

m-benzyne

UHF

0 1 0 1

1.116 2.216 1.007 2.050

1.755 2.403 0.975 2.008

1.346 2.404 0.000 2.009

1.195 2.754 0.556 2.033

PBE0

In most cases, UHF and PBE0 yield large spin contaminations for the singlet states. It is interesting to observe that UHF, but not PBE0, shows large spin contaminations for all systems in the triplet state. Table 8 presents the predicted singlet−triplet splittings for TMM and all three benzyne isomers. CASSCF and PBE0 yield good agreement with experiment for all systems. SUHF and 8078

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

Article

Table 8. Singlet−Triplet Splittings (in kcal/mol) Computed with Several PHF Methods, UHF and PBE0, for TMM and the p-, o-, and m-Benzynesa UHF TMM p-benzyne o-benzyne m-benzyne MAE a b

D T D T D T D T D T

24.11 23.66 −10.28 −10.08 −14.34 −15.77 10.12 9.41 17.19 16.50

UHF†

SUHF

KSUHF

SGHF

KSGHF

43.88

19.82 19.1 −28.76 −28.7 −50.4 −51.44 −4.47 −5.73 14.53 14.28

18.56 18.04 −25.77 −25.57 −53.67 −54.67 −6.27 −7.14 13.83 13.69

16.69 15.93 −6.22 −6.82 −31.99 −33.35 −12.97 −13.99 4.14 3.59

16.20

−31.55 −27.10 13.05 24.96

−12.43 −43.91 −22.95 4.27

PBE0 19.81 14.40 −2.27 −2.23 −26.03 −28.01 −9.94 −10.36 6.94 5.85

PBE0†

CASSCF(2,2)

exptl

38.92

19.5

16.1b

−4.54

−1.27

−3.8c

−26.03

−26.94

−37.5c

−13.71

−12.28

−21.0c

10.29

6.30

Here, D and T refer to results with the cc-pVDZ and cc-pVTZ basis sets, respectively. UHF† and PBE0† refer to the splitting energies using eq 1. Reference 57. cReference 60.

Figure 1. Mean absolute error for all systems and methods computed.

MAE values) with low computational cost in diatomic molecules and carbenes and silenes. Larger systems, which presumably require an accurate treatment of dynamic correlations (TMM and benzynes), are problematic for these methods; (iv) SGHF and KSGHF display a low MAE for all systems and groups; they can be as reliable as CASSCF for molecular species with biradical character.

determinant, and therefore, the connection to the singleparticle picture is not completely lost. Our results show that even the simplest of the PHF approaches can yield accurate singlet−triplet splittings for diatomic molecules as well as carbenes and silenes. Larger molecules with significant spin contamination, such as TMM and the benzyne isomers, seem to require a balanced treatment of strong and weak correlations. In our PHF hierarchy, these can be accessed by breaking and restoring further symmetries (such as Sz) in a variation-after-projection approach. SGHF and KSGHF provide consistently good performance, with better agreement to experimental values than PBE0 or CASSCF(2,2); the latter is the prototypical method to describe diradical systems. Finally, we have also performed a numerically based geometry optimization for some of the systems with the SUHF method. The optimized geometries are in good agreement with the known experimental parameters as well as with PBE0. We conclude that PHF methods, due to their low computational cost and their inherent multireference character, can be a reliable alternative to predict singlet−triplet splittings in systems where Hartree−Fock or DFT methods converge to solutions with large spin contamination.

4. CONCLUSIONS In this work, we have explored the performance of projected Hartree−Fock methods for the prediction of singlet−triplet splittings of different organic molecules. The accurate determination of singlet−triplet gaps is known to be challenging as this is a property sensitive to subtle electron correlation effects. Previous works have suggested that symmetry-preserving (spin) wave functions are a must in order to predict these gaps accurately. The projected Hartree− Fock methods are an interesting alternative to carry out this task: by optimizing a broken-symmetry determinant in the presence of projection operators, one is guaranteed to both account for strong correlations and preserve the correct spin quantum number. Our recent work has shown that the variation-after-projection scheme can be carried out efficiently, with mean-field computational scaling. Lastly, the wave function in PHF methods is fully determined by a single Slater 8079

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080

The Journal of Physical Chemistry A

■

Article

(27) Yamaguchi, K.; Jensen, F.; Dorigo, A.; Houk, K. N. Chem. Phys. Lett. 1988, 149, 537−542. (28) Yamaguchi, K.; Takahara, Y.; Fueno, T.; Houk, K. N. Theor. Chim. Acta 1988, 73, 337−364. (29) Saito, T.; Thiel, W. J. Phys. Chem. A 2012, 116, 10864−10869. (30) Hratchian, H. P. J. Chem. Phys. 2013, 138, 101101. (31) Saito, T.; Nishihara, S.; Nakanishi, Y.; Kitagawa, Y.; Kawakami, T.; Yamanaka, S.; Okumura, M.; Yamaguchi, K. J. Phys. Chem. A 2010, 114, 7967−7974. (32) Saito, T.; Nishihara, S.; Kataoka, Y.; Nakanishi, Y.; Matsui, T.; Kitagawa, Y.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 2009, 483, 168−171. (33) Jiménez-Hoyos, C. A.; Henderson, T. M.; Tsuchimochi, T.; Scuseria, G. E. J. Chem. Phys. 2012, 136, 164109. (34) Samanta, K.; Jiménez-Hoyos, C. A.; Scuseria, G. E. J. Chem. Theory Comput. 2012, 8, 4944−4949. (35) Rodríguez-Guzmán, R.; Schmid, K. W.; Jiménez-Hoyos, C. A.; Scuseria, G. E. Phys. Rev. B 2012, 85, 245130. (36) Rodríguez-Guzmán, R.; Jiménez-Hoyos, C. A.; Schustki, R.; Scuseria, G. E. Phys. Rev. B 2013, 87, 235129. (37) Rivero, P.; Jiménez-Hoyos, C. A.; Scuseria, G. E. J. Phys. Chem. B 2013, DOI: 10.1021/jp401478v. (38) Löwdin, P.-O. Phys. Rev. 1955, 97, 1509−1520. (39) Mayer, I. Adv. Quantum Chem. 1980, 12, 189−262. (40) Jiménez-Hoyos, C. A.; Henderson, T. M.; Scuseria, G. E. J. Chem. Theory Comput. 2011, 7, 2667−2674. (41) Noodleman, L. J. Chem. Phys. 1981, 74, 5737. (42) Noodleman, L.; Davidson, E. R. Chem. Phys. 1987, 109, 131. (43) Ernzerhof, M.; Scuseria, G. E. J. Chem. Phys. 1999, 110, 5029− 5036. (44) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158−6159. (45) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; et al. Gaussian 09, revision H.1; Gaussian, Inc.: Wallingford, CT, 2009. (46) Pulay, P.; Hamilton, T. P. J. Chem. Phys. 1988, 88, 4926. (47) Huber, K. P.; Herzberg, G. Molecular Spectra and Molecular Structure. IV. Constants of Diatomic Molecules; Van Nostrand Reinhold: New York, 1979. (48) Berkowitz, J.; Greene, J. P.; Cho, H.; Rušcǐ ć, B. J. Chem. Phys. 1987, 86, 1235. (49) Bauschlicher, C. W.; Taylor, P. R. J. Chem. Phys. 1986, 85, 6510. (50) Standard, J. M.; Steidl, R. J.; Matthew, C. B.; Robert, W. Q. J. Phys. Chem. A 2011, 115, 1243. (51) Schwartz, R. L.; Davico, G. E.; Ramond, T. M.; Lineberger, W. C. J. Phys. Chem. A 1999, 103, 8213−8221. (52) Schwartz, M.; Marshall, P. J. Phys. Chem. A 1999, 103, 7900− 7906. (53) Chang, W.-Z.; Hsu, H.-J.; Chang, B.-C. Chem. Phys. Lett. 2005, 413, 25. (54) Irikura, K. K.; Goddard, W. A.; Beauchamp, J. L. J. Chem. Soc. 1992, 114, 48−51. (55) Jensen, P.; Bunker, P. R. J. Chem. Phys. 1988, 89, 1327. (56) Rao, D. R. J. Chem. Mol. Spectrosc. 1970, 34, 284. (57) Wenthold, P. G.; Hu, J.; Squires, R. R.; Lineberger, W. C. J. Am. Chem. Soc. 1996, 118, 475−476. (58) Cramer, C. J.; Smith, B. A. J. Phys. Chem. 1996, 100, 9664− 9670. (59) Winkler, M.; Sander, W. J. Phys. Chem. A 2001, 105, 10422− 10432. (60) Wenthold, P. G.; Squires, R. R.; Lineberger, W. C. J. Am. Chem. Soc. 1998, 120, 5279−5290.

AUTHOR INFORMATION

Corresponding Author

*(P.R.) E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Department of Energy, Office of Basic Energy Sciences, Grant No. DE-FG02-09ER16053. G.E.S. is a Welch Foundation Chair (C-0036). We are grateful to the two anonymous referees for their very helpful suggestions.

■

REFERENCES

(1) Hubbard, J. Proc. R. Soc. London, Ser. A 1963, 276, 238−257. (2) Gronert, S.; Keeffe, J. R.; O’Ferrall, R. A. M. J. Am. Chem. Soc. 2011, 133, 3381−3389. (3) Knipe, A. C. Organic Reaction Mechanisms; John Wiley & Sons: London, U.K., 2008. (4) Yu, J.; Lammi, R.; Gesquiere, A. J.; Barbara, P. F. J. Phys. Chem. 2005, 109, 10025−10034. (5) Shi, Z.; Simmons, C. B.; Prance, J. R.; Gamble, J. K.; Friesen, M.; Savage, D. E.; Lagally, M. G.; Coppersmith, S. N.; Eriksson, M. A. Appl. Phys. Lett. 2011, 99, 233108. (6) Zhao, Z. X.; Hou, C. Y.; Zhang, H. X.; Sun, C. J. Chem. Phys. 1981, 74, 4566. (7) Roos, B. O.; Linse, P.; Siegbahn, P. E. M.; Blomberg, M. R. A. Chem. Phys. 1982, 66, 197−212. (8) Knowles, P. J.; Sexton, G.; Handy, N. C. Chem. Phys. 1982, 72, 337−347. (9) Cramer, C. J.; Dulles, F. J.; Storer, J. W.; Worthington, S. E. Chem. Phys. Lett. 1994, 218, 387−394. (10) Standard, J.; Steidl, R. J.; Beecher, C.; Quandt, R. J. Phys. Chem. A 2011, 115, 1243−1249. (11) González-Luque, R.; Merchan, M.; Roos, B. O. Z. Phys. D 1996, 36, 311−316. (12) Standard, J. M.; Quandt, R. W. J. Phys. Chem. A 2003, 107, 6877−6881. (13) Bobrowski, M.; Liwo, A.; Oldziej, S.; Jeziorek, D.; Ossowski, T. J. Am. Chem. Soc. 2000, 122, 8112−8119. (14) Garcia, V. M.; Castell, O.; Reguero, M.; Caballol, R. Mol. Phys. 1996, 87, 1395−1404. (15) Li, X. Z.; Paldus, J. Z. Phys. D 1996, 36, 311−316. (16) Paldus, J. Relativistic and Electron Correlation Effects in Molecules and Solids; Malli, G. L., Ed.; Springer: New York, 1994; pp 207−282. (17) Parr, R. G.; Yang, W. Density Functional Theory of Atoms and Molecules; Oxford University Press: Oxford, U.K., 1989. (18) Dreizler, D. M.; Gross, E. Density Functional Theory. An Approach to the Quantum Many Body Problem; Springer: Berlin, Germany, 1990. (19) Lim, M. H.; Worthington, S. E.; Dulles, F. J.; Cramer, C. J. ACS Symp. Ser. 1996, 629, 402−422. (20) Wang, J.; Becke, A. D.; Smith, U. H. J. Chem. Phys. 1995, 102, 3477−3480. (21) de Visser, S. P.; Filatov, M.; Shaik, S. Phys. Chem. Chem. Phys. 2000, 22, 5046−5048. (22) Grafenstein, J.; Cremer, D.; Kraka, E. Chem. Phys. Lett. 1998, 288, 593−602. (23) Grafenstein, J.; Cremer, D. Chem. Phys. Lett. 2000, 316, 569− 577. (24) Garza, A. J.; Jiménez-Hoyos, C. A.; Scuseria, G. E. J. Chem. Phys. 2013, 138, 134102. (25) Tsuchimochi, T.; Scuseria, G. E. J. Chem. Phys. 2011, 134, 064101. (26) Saito, T.; Nishihara, S.; Yamanaka, S.; Kitagawa, Y.; Kawakami, T.; Yamada, S.; Isobe, H.; Okumura, M.; Yamaguchi, K. Theor. Chem. Acc. 2011, 130, 739−748. 8080

dx.doi.org/10.1021/jp405755z | J. Phys. Chem. A 2013, 117, 8073−8080