9664
J. Phys. Chem. 1996, 100, 9664-9670
Trimethylenemethane. Comparison of Multiconfiguration Self-Consistent Field and Density Functional Methods for a Non-Kekule´ Hydrocarbon Christopher J. Cramer* and Bradley A. Smith Department of Chemistry and Supercomputer Institute, UniVersity of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455-0431 ReceiVed: December 13, 1995X
Relative energies for different multiplets of trimethylenemethane and methylenecyclopropane are calculated at MCSCF, CASPT2N, and DFT levels of theory. Comparison to the experimentally measured heat of formation for the 3A2′ state and to the experimentally measured 3A2′-1A1 gap permits analysis of the relative importance of active space and basis set selection in the multiconfigurational methods. Such comparison also reveals that while DFT accurately treats triplet trimethylenemethane, there is a limitation in present DFT functionals with respect to accurately treating nondynamic correlation effects for closed-shell singlets in a molecule characterized by degenerate frontier molecular orbitals. Implications for calculations on larger systems are discussed.
1. Introduction Trimethylenemethane (TMM) is the simplest non-Kekule´ hydrocarbon.1-7 The designation non-Kekule´ implies that it is a conjugated system that cannot be represented by any resonance structure containing n double bonds derived from its 2n pi electrons. In 1948, Coulson assigned credit to Moffitt for having first noted the unusual electronic structure of TMM.8 LonguetHiggins later provided a detailed analysis of the non-Kekule´ class of molecules in 1950.9 Experimental identification of TMM was reported by Dowd in 1966,10 and verification of its triplet ground state in an 88 K matrix was ultimately accomplished in 1976.11 A zeroth-order analysis predicts a triplet ground state for TMM insofar as its electronic structure is characterized by two degenerate (pi) frontier molecular orbitals (FMO’s); i.e., Hund’s rule applies and this molecule is a diradical (Figure 1). However, double occupation of one of the two degenerate FMO’s to create a planar closed-shell singlet state can be stabilized by a Jahn-Teller12 distortion. Moreover, a different Jahn-Teller distortion stabilizes a planar open-shell singlet state. Finally, a 90° rotation of one methylene in TMM leads to a structure of C2V symmetry that presumably may be described by low-energy triplet and open-shell singlet wave functions (although a closed-shell singlet can be constructed, it is expected to be less energetically relevant). This situation is summarized in Figure 2. While Dowd unambiguously identified the planar triplet as the ground state of TMM in 1976,11 only estimates of the gaps separating this state from the various singlet states were available from subsequent experimental studies.13,14 Theoretical estimates of the multiplet splittings also began to appear following isolation of the triplet.15-29 However, properly accounting for nondynamic and dynamic electron correlation effects, and moreover doing so in a consistent fashion for multiple electronic states, is not a trivial undertaking (and this was especially true for the time period immediately following the experimental isolation). The best approaches tended to employ multiconfiguration self-consistent-field (MCSCF) methods followed by a multireference configuration interaction (MRCI) calculation. However, it remains impractical even today to include all of the valence electrons in the MCSCF X
Abstract published in AdVance ACS Abstracts, May 1, 1996.
S0022-3654(95)03697-5 CCC: $12.00
Figure 1. The four pi orbitals for planar TMM shown on the right. A 90° rotation of one methylene transforms the 2b1 orbital to the b2 orbital shown at lower left. The possible wave functions (with the spin functions ignored) that can be formed using the two highest energy electrons are listed. The 3B2 state carries the label 3A2′ in full D3h symmetry.
active space, raising questions about which space is most evenhanded in its treatment of the various multiplets. It obviously follows that full CI is similarly impractical. As a result, the quality of the calculated multiplet state energies is difficult to determine. This ambiguity has been recently eliminated, at least for the case of the planar closed-shell singlet. Photoelectron spectroscopy (PES) of the TMM radical anion has unambiguously determined the 1A1 state to lie 16.1 ( 0.1 kcal/mol above the 3A ′ ground state.30 Although the best theoretical predictions 2 place the open-shell singlet slightly lower in energy, it is not observed in the photoelectron spectrum. This is consistent with theoretical findings that the rotated 1B1 structure is lower in energy than the planar 1B2 structure, where the latter is presumably a (hyper)saddle point on the potential energy surface. As a result, Franck-Condon overlap is insufficient for spectroscopic observation. The PES experiment also establishes the heat of formation of 3A2′ TMM to be 70 ( 3 kcal/mol. Since the heat of formation of methylenecyclopropane (MCP) is known (47.9 ( 0.4 kcal/mol),31 the accuracy of different theoretical methods can be assessed by calculating the enthalpy difference between 3A2′ TMM and MCP. © 1996 American Chemical Society
Trimethylenemethane
J. Phys. Chem., Vol. 100, No. 23, 1996 9665
Figure 2. Gross structures of the relevant multiplets of TMM and MCP. Full geometric details can be found in Table 1.
The chief purpose of the present article is to explore the utility of density functional theory (DFT) for the prediction of TMM multiplet energies. We have demonstrated DFT’s high accuracy in the prediction of singlet-triplet gaps for a variety of singlecenter diradicals (e.g., carbenes, nitrenium ions, etc.)sin general (2 kcal/mol by comparison to either experiment or converged multiconfigurational calculations.32-37 Since DFT is intrinsically a single-configuration methodology, such accuracy implies that nondynamic correlation effects, which are especially critical in the evaluation of closed-shell singlet energies in single-center diradicals, are well accounted for by construction of the employed functional of the electron density. TMM thus represents a particular challenge for DFT, since for the 1A1 state the largest portion of the nondynamic correlation arises from the (nearly) degenerate FMO’s. In addition, it is not straightforward to employ DFT methodology for an intrinsically twoconfiguration system, like an open-shell singlet.35,36,38 We explore here the utility of the sum method38 for this purpose. In order to calibrate various aspects of our DFT work, we also examine MCSCF calculations followed by multireference second-order perturbation theory calculations (CASPT2N39-41) using active spaces of different sizes. These latter calculations parallel similar high-level work of Hrovat and Borden, who have additionally examined MRCI methods.42 Applications of the CASPT2 method to other non-Kekule´ systems have appeared.43,44 Section 2 describes the methodology, section 3 presents results from the electronic structure calculations, and section 4 provides a comparison of the different computational techniques in the context of experiment and earlier work. Finally, section 5 provides some general conclusions. 2. Computational Methods MCSCF calculations were of the complete active space variety. Active spaces were constructed from either the FMO’s (2,2), from the full set of molecular orbitals dominated by the carbon 2p orbitals (4,4), or from the latter set plus the six valence orbitals most closely associated with the C-C sigma framework (10,10). Additional correlation was accounted for at the CASPT2N level39-41 using the multireference wave functions. Density functional calculations employed Becke’s gradientcorrected expression for exchange energy45,46 and Perdew’s gradient-corrected expression for the correlation energy,47 which we have abbreviate BP86. We explored a number of different functionalssmodern (i.e., gradient-corrected) functionals all gave similar results to within about 2 kcal/mol in the relative energies, and for the sake of brevity we will not discuss them further. DFT triplet energies were calculated in straightforward fashion. Open-shell singlet energies were calculated using the sum method. That is, we double the difference in energy between orbital representations of the density having S2 expectation values of 2 (i.e., corresponding to a triplet for the wave function that could be formed from those orbitals) and of 1 (i.e.,
corresponding to a 50:50 singlet-triplet wave function) and add that value to the triplet energy. A more detailed presentation of this approach can be found elsewhere.35,36,38 Early applications of DFT spin-annihilation techniques35 to TMM were problematic, and were not further explored. Finally, for the closed-shell singlet, energies were calculated either with enforced double occupation of a single frontier Kohn-Sham orbital or with two electrons worth of paired spin density “smeared” equally over both frontier orbitals. The latter approach results in a significantly decreased energy for the 1A1 state. Most calculations employed either the cc-pVDZ or cc-pVTZ basis sets of Dunning.48 The presented DFT, MCSCF(10,10), and CASPT2N(10,10) relative energies are converged to within 1 kcal/mol based on analysis of the results obtained with these basis sets. Molecular geometries were optimized in the standard fashion for most cases. MCSCF geometries were optimized for all choices of active space; however, little change was observed on going from one active space to the next, so only the (10,10) geometries are reported below. Open-shell singlet geometries were not optimized at the DFT level since the electronic energies for these states were calculated indirectly (geometries for the 50:50 states were optimized at this level). Analytic frequency calculations at the DFT level indicate MCP and the 3A2′ state of TMM to be local minima and zero-point energies were calculated therefrom; furthermore, the 3B1 structure is calculated to be a transition state for methylene rotation. Analytic DFT second derivatives are not available for the other (pure) multiplets. Frequencies calculated numerically at the MCSCF level are associated with considerable uncertainty since, in the absence of a complete valence active space, they suffer from energetic anomalies associated with orbital symmetry changes accompanying geometry changes that break C2V symmetry. Not surprisingly, this effect is largest for those frequencies associated with soft modes like torsions. On the other hand, these same modes have little impact on calculated zero-point vibrational energies (ZPVE), so these calculations can still provide some useful information. MCSCF and CASPT2N calculations were performed with the MOLCAS version 3.0 program suite.49 MCSCF calculations were also carried out using GAMESS version 08/11/94.50 The two programs gave agreement on absolute energies to within 0.5 mhartree for MCSCF(2,2) and (4,4) calculations; however, the (10,10) energies differed by up to 25 mhartree. Since the relative energies calculated with MOLCAS at the MCSCF(10,10) level differed only marginally from the MCSCF(4,4) level but the GAMESS relative energies changed significantly, we deem the latter untrustworthy (perhaps reflecting a poor assignment of orbitals to the active space, although considerable effort was made to remedy any such problems). Density functional calculations were carried out both with GAUSSIAN92/DFT51 and ADF version 1.1.4.52 Only the latter program permits the smeared-density calculation for the 1A1 state; relative energies
9666 J. Phys. Chem., Vol. 100, No. 23, 1996
Cramer and Smith
TABLE 1: Geometrical Data for TMM and MCP Structuresa
3A
MCSCF r1 r2 r3 r4 r5
1.438 1.438 1.081 1.081 1.081
a1 a2 a3 a4 d1(HCCC)
2′
1A
DFT 1.424 1.424 1.100 1.100 1.100
MCSCF 1.370 1.496 1.082 1.080 1.080
1B 2
1
DFT
b
1.381 1.450 1.104 1.097 1.098
MCSCF
3B 1
DFTc
MCSCF
Bond Lengths, Å 1.539 1.459 1.411 1.409 1.079 1.100 1.081 1.100 1.081 1.100
1B 1
DFT
1.485 1.417 1.081 1.081 1.082
1.486 1.404 1.102 1.100 1.102
120.0 120.9 120.9 120.9
120.0 121.0 121.0 121.0
121.1 121.2 120.2 120.9
114.1 121.4 119.1 122.0
Valence Angles, deg 119.0 119.6 119.5 120.3 120.8 120.8 121.4 121.2 121.3 120.6 120.9 120.7
119.2 120.9 121.3 120.7
0.0
0.0
0.0
0.0
Dihedral Angle, deg 0.0 0.0 90.0
90.0
a
MCP
MCSCFd
DFTc
1.520 1.416 1.081 1.082 1.081
1.490 1.403 1.111 1.100 1.102
MCSCF
DFT
1.348 1.498 1.083 1.084
1.335 1.478 1.102 1.104
119.2 120.7 120.6 121.2
119.2 120.9 121.3 120.6
148.4 121.0 118.3
148.4 121.2 118.8
90.0
90.0
73.6
73.3
b
Geometries were optimized at the MCSCF(10,10)/cc-pVDZ and BP86/cc-pVDZ levels of theory. Refers to restricted calculation enforcing double occupation of HOMO. c Refers to 50:50 state. d At the MCSCF(4,4) and MCSCF(10,10) levels, this geometry relaxes very slightly to one of Cs symmetry (1A′′); see text.
for the other states from the two programs agreed to within 0.4 kcal/mol. Calculations with ADF employed the package’s polarized split-valence basis sets IV and V (which appear to be of similar quality to the Dunning basis sets on the basis of the observed agreement in relative multiplet energies).53 3. Results Table 1 provides geometrical data for all of the relevant multiplets at the MCSCF(10,10)/cc-pVDZ and BP86/cc-pVDZ levels. Geometries calculated with smaller MCSCF active spaces do not deviate much from those listed in Table 1. Nor, on the basis of a comparison to the MCSCF(4,4)/STO-3G results of Feller et al.,21 is there much of a dependence of geometries on basis set quality. Where direct comparison is possible, the DFT geometries agree reasonably closely with the MCSCF results except for the 1A1 state, where the DFT calculation enforcing a doubly occupied HOMO predicts the longer C-C bond to be 0.046 Å shorter and also predicts a smaller CCC valence angle by 7° (the smeared-density calculation, on the other hand, gives a geometry that is essentially D3hsvide infra). For the open-shell singlets 1B2 and 1B1, the DFT geometries refer to the 50:50 mixed state. Since the 3B1 and 1B1 geometries are similar at the MCSCF level, the mixed-state DFT geometry is a good match for the 1B1 geometry. For the planar B2 states, on the other hand, there is a large difference between the singlet and triplet geometries (note that in the C2V point group the 3A2′ state correlates with 3B2). The mixed-state DFT geometry is indeed about halfway in between the MCSCF singlet and triplet geometries for most degrees of freedom. The 50:50 geometry is clearly a less satisfactory approximation for the B2 state than for the B1. At the MCSCF(4,4) level, the numerically evaluated Hessian for the 1B1 state has one negative eigenvalue, with an associated frequency of 158i. The eigenvector that is associated with this frequency pyramidalizes the rotated methylene. However, under the constraints of Cs symmetry, the optimized pyramidal structure (no negative eigenvalues in the numerically evaluated
TABLE 2: Relative Multiplet Energies as a Function of Active Space and Basis Seta,b method MCSCF(2,2)/cc-pVDZ MCSCF(4,4)/cc-pVDZ MCSCF(10,10)/cc-pVDZ MCSCF(10,10)/cc-pVTZ CASPT2N(2,2)/cc-pVDZ CASPT2N(4,4)/cc-pVDZ CASPT2N(10,10)/cc-pVDZ CASPT2N(10,10)/cc-pVTZ zero-point energiese
3A
2′
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
c
1A
1
10.5 19.5 19.2 18.9 23.0 20.1 19.1 19.1 -1.6
1B 2
3B 1
1B d 1
MCP
16.7 19.9 19.5 19.3 19.7 19.8 19.1 19.1 -2.9
13.9 13.4 13.8 13.8 15.1 14.1 13.9 14.4 -1.3
15.1 14.7 15.9 15.8 16.9 15.8 15.7 16.1 -1.2
-29.5 -21.5 -18.9 -19.8 -17.8 -20.7 -20.9 -21.9 2.4
a Relative energies in kcal/mol. b MCSCF geometries were fully optimized with the cc-pVDZ basis set; CASPT2N calculations were at the corresponding MCSCF geometries. c Absolute energies (Eh) for this column are -154.879 71, -154.918 57, -154.995 09, -155.035 62, -155.424 57, -155.420 94, -155.425 36, -155.630 73, 0.077 30. d At the MCSCF(4,4) and MCSCF(10,10) levels, this geometry relaxes very slightly to one of Cs symmetry (1A′′), but the energetic effect of this relaxation is less than 0.1 kcal/mol. e Calculated at the MCSCF(4,4)/ cc-pVDZ level with all frequencies scaled by 0.92.
Hessian) has HCCC dihedral angles at the rotated methylene of (80°; i.e., the extent of pyramidalization is only 10° away from C2V. The remaining geometric degrees of freedom are essentially unperturbed from those in Table 1 (bond lengths change by less than 0.001 Å and angles by less than 0.2°). Moreover, the energetic consequence of pyramidalization is to lower the relative energy by less than 0.1 kcal/mol. The same results are obtained at the MCSCF(10,10) level, although numerical frequency analysis was not carried out at this level. In essence, then, even if there is truly a double-well potential corresponding to methylene inversion, the first vibrational level may lie above the barrier, and it is more meaningful to simply regard the inversion mode as being extremely soft. For ease of discussion, we will continue to employ the C2V geometry of the 1B1 state in the ensuing discussion. Table 2 illustrates the effect of active-space size on the relative multiplet energies calculated at the MCSCF and CASPT2N levels. Fully optimized geometries were used in every case, but the changes in geometry from one level to the
Trimethylenemethane
J. Phys. Chem., Vol. 100, No. 23, 1996 9667
TABLE 3: BP86 Relative Multiplet Energies for MCSCF and DFT Geometriesa geometry MCSCF(10,10)/cc-pVDZ BP86/cc-pVDZ
3A
2′
0.0 0.0
b
1A c 1
1B d 2
27.6 17.1 23.0 16.9
3B
1
1B d 1
MCP
13.4 14.9 -26.0 13.6 15.4 -26.1 (12.6)e (-23.6)e
a Relative energies in kcal/mol using the cc-pVTZ basis set. Absolute energies (Eh) for this column are -155.971 42, -155.972 21. c Energies refer to smeared-density calculations. If double occupation of the lower energy FMO (after Jahn-Teller distortion) is enforced the relative energies in this column are 39.0 and 32.7. d Energies calculated using the sum method; DFT geometries are from 50:50 state optimizations. e Includes BP86/cc-pVDZ zero-point energy correction relative to 3A2′ state. b
next are small, and the associated changes in relative energy are similarly small (typically less than 0.4 kcal/mol). Table 3 contains BP86/cc-pVTZ relative multiplet energies calculated for both the DFT and MCSCF geometries. BP86 relative energies calculated with the smaller cc-pVDZ basis set are within 5% of the valence-triple-ζ results. We have observed this rapid convergence with respect to basis set size in DFT relative multiplet energies for single-center diradicals as well.32-34,36,37 The open-shell singlet energies for DFT geometries were calculated using the 50:50 mixed-state geometry. This is obviously an approximation, as discussed further in the next section. In principle an optimized open-shell singlet geometry can be located at the DFT level, but analytic gradients are not available to facilitate this process. The ZPVE’s calculated at the BP86/cc-pVDZ level for MCP and 3A2′ TMM are 51.4 and 48.9 kcal/mol, respectively. At the corresponding restricted Hartree-Fock (RHF) and unrestricted Hartree-Fock (UHF) levels, after scaling the frequencies by the standard factor of 0.89,54 these ZPVE’s are calculated to be 50.5 and 46.9 kcal/mol, respectively. The UHF wave function for 3A2′ TMM is moderately spin-contaminated (〈S2〉 ) 2.199), and this may explain some of the sizable difference between BP86 and UHF for this state. ZPVE’s calculated at the MCSCF (4,4)/cc-pVDZ level were also obtained. The unscaled ZPVE’s for MCP and 3A2′ TMM are 55.3 and 52.7 kcal/mol, respectively. If we assume the ZPVE for MCP as calculated by either DFT or scaled RHF to be trustworthy, a good scale factor to apply to the MCSCF frequencies55 would appear to be 0.92; this gives scaled ZPVE’s for MCP and 3A2′ TMM of 50.9 and 48.5 kcal/mol, respectively. The difference between the two ZPVE’s is 2.4 kcal/mol, in good agreement with the 2.5 kcal/mol difference calculated at the BP86 level. This scaling factor was applied to obtain all of the ZPVE’s listed in Table 2. Finally, we take note of one practical aspect of these calculations, namely the resources required in terms of CPU time and disk usage. In general, with the same basis set, the DFT calculations take 1 order of magnitude less CPU time than the CASPT2N calculations. The MCSCF and CASPT2N calculations are moderately demanding in terms of disk resources (about 1 GB is typical), while direct DFT calculations use only about 1% of this storage. These differences would be expected to increase for larger molecules. 4. Discussion We begin with an analysis of various aspects of the different modeling approaches, then discuss the correspondence between the best theoretical and experimental results, and finally conclude with some suggestions with respect to modeling larger non-Kekule´ systems.
The calculated relative energies in Table 2 indicate that the minimum active-space size required for accuracy with TMM is (4,4). As might be expected, perturbation theory in the form of CASPT2N(2,2) calculations can be used to partially make up for the inadequacy of the MCSCF(2,2) space, but the results remain rather unsatisfactory compared to those obtained with larger active spaces. Thus, at the CASPT2N(2,2) level the errors in the MCSCF(2,2) predicted relative energies are reduced by about 50% compared to the most complete calculations, but they still remain as large as 4 kcal/mol. As for the (4,4) calculations, it is not entirely obvious that even this size active space should be adequate for comparison of planar and twisted geometries of TMM. For planar geometries, it is clear that the (4,4) space represents the complete valence pi space. These orbitals are orthogonal to the sigma framework and are unique and distinct. For the twisted geometries, on the other hand, there is no longer a clean separation of pi and sigma orbitals. Nevertheless, the TMM relative energies appear well converged with respect to increasing the active space size, as judged by comparison to the (10,10) results. The latter include the entire heavy atom sigma framework in the active space. MCP, however, is not well treated with only a (4,4) active space at the MCSCF level. When one correlates the TMM orbitals with MCP, the (4,4) space includes the pi orbitals of the exocyclic double bond and also the sigma orbitals of the basal C-C bond in the cyclopropane. Obviously, it is unbalanced to include only the basal bond of the alicyclic ring. When the remaining C-C sigma bonds are included in the (10,10) calculations, a more balanced comparison to the TMM structures results. The CASPT2N calculations using either the (4,4) or the (10,10) active spaces agree closely with each other for MCP. That is, in this case the CASPT2N method not only accounts for dynamic electron correlation but also makes up for the imbalance in the active space. On the whole, the additional correlation effects captured by the CASPT2N calculations are essentially constant across the range of TMM multiplets, especially when the active space is increased to (10,10) (which is to say the relative energies for the various multiplets change only very slightly). Again, however, this is not the case for MCP, where, compared to the MCSCF level, the CASPT2N calculations reduce the energy of MCP by 2 kcal/mol relative to the TMM states. This additional reduction may be associated with dynamic electron correlation, which might be expected to be more important for the small ring present in MCP, but other considerations (e.g., different C-H hybridizations compared to the TMM states) probably enter as well. Finally, increasing the basis set size from cc-pVDZ to ccpVTZ has negligible effect on the relative TMM energies. It does, however, lower the relative energies of the twisted B1 states by 0.5 kcal/mol and that of MCP by 1 kcal/mol. McKee has discussed the importance of polarization functions and expanded valence spaces for accurate energies in the cyclopropylcarbinyl system, and similar considerations appear to be operative here for MCP.56 As for the relative DFT energies, with the exception of the 1A state they appear insensitive to choice of geometry, i.e., 1 MCSCF vs BP86. This is unsurprising for those states where the geometries are similar (see Table 1) but less expected for the 1B2 state, where the MCSCF geometry differs markedly from the BP86 mixed-state B2 geometry. At least at the DFT level, the potential surface appears to be quite flat in directions associated with the Jahn-Teller distortion observed for this state. Evidently, the potential energy surface for the open-shell
9668 J. Phys. Chem., Vol. 100, No. 23, 1996 singlet is quite flat for seVeral internal coordinates. As noted above, the inversion potential for the rotated methylene in the 1B state is very soft. Moreover, after accounting for differential 1 ZPVE’s, the barrier to methylene rotation (i.e., the energy difference between the 1B1 and 1B2 states) is only 1.3 kcal/mol at the best CASPT2N level. The relative energy of the 1A1 state, on the other hand, is considerably more sensitive to geometry. The calculated absolute energies for this state are discussed in more detail below, but we note here that the relative energy calculated at the BP86 level using the smeared-density geometry (which is very nearly D3h) is 4.6 kcal/mol lower than that calculated by the same method at the MCSCF geometry. To assess possible geometric effects of electron correlation not accounted for at the MCSCF(10,10) level, we have carried out CASPT2N (10,10)/cc-pVTZ calculations at the BP86 geometry. The calculated energy of the 1A1 state at this geometry increases by 3.0 kcal/mol, suggesting that the DFT results with respect to geometry are not trustworthy for this multiplet. As for the magnitudes of the DFT relative energies in comparison to CASPT2N, we begin with a discussion of those states to which DFT may be unambiguously applied, i.e., the triplets and MCP. For the ZPVE-corrected rotation of one methylene in the triplet, DFT predicts a barrier of about 12.6 kcal/mol, while the CASPT2N prediction is 13.1 kcal/mol. These barriers are similar in magnitude, and it is not presently clear which may be the more accurate. A larger difference is observed for the prediction of the relative energy of MCP. In this case, BP86 predicts MCP to be lower in relative energy by about 4 kcal/mol compared to CASPT2N. For the latter comparison, we may take advantage of the experimentally determined heats of formation of MCP31 and 3A ′ TMM30 to assess the accuracy of the two theoretical levels. 2 As noted in the Introduction, the difference in the heats of formation for these two molecules is 22.1 ( 3.0. The calculated differences in ∆Hf,0 are 19.5 and 23.6 kcal/mol at the CASPT2N and BP86 levels, respectively. Both results lie within the experimental margin of error. However, CASPT2N has been established to selectively stabilize high-spin states over lowspin statessAndersson et al. note that the effect is usually from 3 to 5 kcal/mol per pair of electrons14,57 (in recent work36 we have found instances where this overstabilization can be as much as twice as large in magnitude). Thus, the discrepancy between the BP86 and CASPT2N results may be entirely associated with this effect. Andersson et al.41 have considered a different form for the zeroth-order Hamiltonian employed in the multireference perturbation theory that is designed to reduce this imbalance between high- and low-spin systems, and they refer to the resulting method as CASPT2G (of which there are two variations). Hrovat and Borden,42 however, have not observed CASPT2G calculations on the TMM system to predict relative energies improved over CASPT2N, at least based on comparison to the experimental 3A2′-1A1 splitting. For the open-shell singlets, BP86 predicts relative energies that are 1-2 kcal/mol lower than those found at the CASPT2N level. This agreement is close enough that it is not obvious which level (if either) is more accurate. In any case, the agreement we observe is yet another example of the general utility of the DFT sum method36,38,58,59 in predicting open-shell singlet energies. Finally, the 1A1 state of TMM is obviously handled poorly at the BP86 level.60 The energies predicted using a smeareddensity approach, where the highest-energy pair of electrons fractionally occupies both FMO’s, are 4 and 7 kcal/mol higher than the CASPT2N predictions and experiment, respectively.
Cramer and Smith
Figure 3. Energetics for C2 conrotatory closure of 1A1 TMM to MCP at the DFT and CASPT2N levels. The torsional coordinate is defined as the dihedral angle between the plane(s) bisecting the rotating methylene group(s) and the plane(s) containing the three unique carbon atoms. The solid line was obtained at the MCSCF(10,10)/cc-pVDZ level, the dotted line at the restricted BP86/cc-pVDZ level, and the dashed line at the smeared-density level equivalent to BP86/cc-pVDZ.
The situation is much worse if the highest energy pair of electrons is forced to occupy the HOMO exclusively, and this observation allows us to suggest the source of this deficiency. DFT, in principle, accounts for all electron correlation effects,61 both nondynamic and dynamic, by construction of the Kohn-Sham operator. However, the exact Kohn-Sham operator is not known for most systems, and instead approximations (like BP86) are employed. In previous work, we have observed most modern functionals of the density to accurately predict relative energies between singlet and triplet states of singlecenter diradicals.32-34,36 Since nondynamic correlation effects are larger for the singlets than the triplets in such cases, it is apparent that the functionals properly account for this factor. However, in none of these cases have the FMO’s been degeneratesrather, there has been a small to moderate HOMOLUMO gap. It is evident that the ability of DFT to account completely for nondynamic correlation begins to break down as the HOMO-LUMO gap approaches zero. Smearing the paired spin density over both orbitals provides some relief, but it is ultimately unable to completely account for the full correlation energy. To more quantitatively assess this contention, we have examined the conrotatory opening of MCP to form 1A1 TMM (Figure 3). Feller et al.21 have provided a complete twodimensional analysis of the torsional potential energy surface for the TMM/MCP system at a lower level of theory, and our CASPT2N conrotatory coordinate is consistent with their earlier results. The CASPT2N and the restricted (i.e., not employing a smeared density) BP86 coordinates are in good agreement until the torsion angle for the rotating C-C bonds reaches about 60°. At this point, the BP86 coordinate continues to rise while the CASPT2N drops to the 1A1 minimum. In the MCSCF reference wave function for the CASPT2N calculations, the occupation numbers of the formal HOMO and LUMO orbitals are 1.519 and 0.487, respectively. More work remains to be done in order to establish whether these occupation numbers may be considered to be a general borderline for DFT difficulties. Figure 3 illustrates as well the smeared-density results, which are greatly improved over the restricted approach near the planar geometry. Nevertheless, the closed-shell singlet energy remains unphysically high. Moreover, it is clarly
Trimethylenemethane unsatisfactory to be reduced to rather ad hoc decisions about when smeared-density calculations should be done instead of restricted-density calculations. Solutions to this nondynamic correlation problem include the development of improved functionals specifically designed to address this challenge, or alternatively the construction of a “multiconfiguration” formulation62,63 of density functional theory. Having established the trustworthiness of the various levels of theory with respect to their fabrication, we may now turn to a comparison with available experimental results. In particular, the experimental assignment of the adiabatic 3A2′-1A1 gap to be 16.1 ( 0.1 kcal/mol is well supported by the CASPT2N calculations. With our largest basis set and underlying active space, we predict a ZPVE-corrected gap of 17.5 kcal/mol. This modest overestimation may be associated with the observed tendency of the CASPT2N method to stabilize triplets over singlets, but the CASPT2G results of Hrovat and Borden would tend to mitigate against that interpretation. At all levels of theory the present calculations suggest that the open-shell singlet is the lowest-energy singlet state. An experimental measurement of this state’s relative energy is unavailable. On a separate point, Dowd and Chow13,14 followed the disappearance of the 3A2′ TMM electron spin resonance signal in a cold matrix, and suggested that the process quenching that signal occurred with an activation barrier on the order of 7 kcal/mol. It is clear from both experiment and theory that electronic excitation of planar TMM from the 3A2′ state to the 1A state cannot account for this observation. Early MCSCF 1 work by Feller et al.21 determined that the energy of the 1A1 state rises with the rotation of a single methylene group, so crossing to the closed-shell singlet surface cannot take place at lower energy via such a conformational change. Furthermore, we find at the CASPT2N(10,10)/cc-pVDZ level that the triplet surface crosses the closed-shell singlet surface along the reaction coordinate for conrotatory closure at 21.4 kcal/mol (relative to the 3A2′ minimum). The full torsional analysis of Feller et al.21 at a lower level indicates that crossing points significantly lower in energy do not exist. Finally, our calculated open-shell energies suggest that a crossing to this surface also will not occur with an activation barrier much below 14 kcal/mol (and, of course, the Franck-Condon factors for such a crossing are far from optimal). It seems clear then that some other process, possibly associated with a bimolecular reaction taking place in the matrix, must be invoked to rationalize the observed triplet quenching. 5. Conclusions It remains a challenge to accurately calculate multiplet splittings in non-Kekule´ systems. TMM proves to be particularly demanding as a test case. CASPT2N calculations must be interpreted with care because of the tendency for that method to stabilize electronic states having unpaired spin over those having less or none at all. Moreover, it is clear that for larger non-Kekule´ molecules the construction of adequate multireference wave functions becomes increasingly costly. Density functional methods, on the other hand, suffer from difficulties in accurately determining open-shell singlet geometries. Moreover, TMM illustrates that the ability of DFT to adequately account for nondynamic correlation effects in closedshell singlets breaks down as the limit of degenerate frontier molecular orbitals is approached. Of course, while they may be present in certain highly symmetric systems,64,65 such orbital degeneracies are not a necessary characteristic of non-Kekule´ molecules, and the potential for DFT to be useful in systems with a modest HOMO-LUMO gap remains.
J. Phys. Chem., Vol. 100, No. 23, 1996 9669 With that in mind, it may be that the most efficient way to perform calculations on large non-Kekule´ molecules that do not have frontier-orbital near degeneracies will be at the DFT// MCSCF level. That is, obtain reasonable geometries for all multiplets using MCSCF with a tractable (small) active space, and then obtain electronic energies employing DFT, where ambiguities associated with consistency in active space construction and accounting for electron correlation are reduced. We will examine this hypothesis in future work. Acknowledgment. It is a pleasure to thank Profs. Weston Borden, Paul Dowd, Mark Gordon, Ken Jordan, and Robert Squires for stimulating discussions. References and Notes (1) Dowd, P. Acc. Chem. Res. 1972, 5, 242. (2) Ovchinnikov, A. A. Theor. Chim. Acta 1978, 47, 297. (3) Berson, J. A. Acc. Chem. Res. 1978, 11, 446. (4) Borden, W. T.; Davidson, E. R. Acc. Chem. Res. 1981, 14, 69. (5) Borden, W. T. In Diradicals; Borden, W. T., Ed.; WileyInterscience: New York, 1982; p 1. (6) Dougherty, D. A. Acc. Chem. Res. 1991, 24, 88. (7) Borden, W. T.; Iwamura, H.; Berson, J. A. Acc. Chem. Res. 1994, 27, 109. (8) Coulson, C. A. J. Chim. Phys. 1948, 45, 243. (9) Longuet-Higgins, H. C. J. Chem. Phys. 1950, 18, 265. (10) Dowd, P. J. Am. Chem. Soc. 1966, 88, 2587. (11) Baseman, R. J.; Pratt, D. W.; Chow, M.; Dowd, P. J. Am. Chem. Soc. 1976, 98, 5726. (12) Jahn, H. A.; Teller, E. Proc. R. Soc. London 1937, A161, 220. (13) Dowd, P.; Chow, M. J. J. Am. Chem. Soc. 1977, 99, 6438. (14) Dowd, P.; Chow, M. J. Tetrahedron 1982, 38, 799. (15) Yarkony, D. R.; Schaefer, H. F., III. J. Am. Chem. Soc. 1974, 96, 3754. (16) Davidson, E. R.; Borden, W. T. J. Am. Chem. Soc. 1977, 99, 2053. (17) Davis, J. H.; Goddard, W. A. J. Am. Chem. Soc. 1977, 99, 4242. (18) Hood, D. M.; Schaefer, H. F., III; Pitzer, R. M. J. Am. Chem. Soc. 1978, 100, 8009. (19) Feller, D.; Borden, W. T.; Davidson, E. R. J. Chem. Phys. 1981, 74, 2256. (20) Dixon, D. A.; Foster, R.; Halgren, T. A.; Lipscomb, W. N. J. Am. Chem. Soc. 1981, 100, 1359. (21) Feller, D.; Tanaka, K.; Davidson, E. R.; Borden, W. T. J. Am. Chem. Soc. 1982, 104, 967. (22) Auster, S. B.; Pitzer, R. M.; Platz, M. S. J. Am. Chem. Soc. 1982, 104, 3812. (23) Borden, W. T.; Davidson, E. R.; Feller, D. Tetrahedron 1982, 38, 737. (24) Lahti, P. M.; Rossi, A. R.; Berson, J. A. J. Am. Chem. Soc. 1985, 107, 2273. (25) Skancke, A.; Schaad, L.; Hess, B. A. J. Am. Chem. Soc. 1988, 110, 5315. (26) Olivella, S.; Salvador, J. Int. J. Quantum Chem. 1990, 37, 713. (27) Radhakrishnan, T. P. Tetrahedron Lett. 1991, 32, 4601. (28) Borden, W. T. Mol. Cryst. Liq. Cryst. 1993, 232, 195. (29) Ichimura, A. S.; Koga, N.; Iwamura, H. J. Phys. Org. Chem. 1994, 7, 207. (30) Wenthold, P. G.; Hu, J.; Squires, R. R.; Lineberger, W. C. J. Am. Chem. Soc. 1996, 118, 475. (31) Cox, J. D.; Pilcher, G. Thermochemistry of Organic and Organometallic Compounds; Academic: New York, 1970. (32) Cramer, C. J.; Dulles, F. J.; Storer, J. W.; Worthington, S. E. Chem. Phys. Lett. 1994, 218, 387. (33) Cramer, C. J.; Dulles, F. J.; Falvey, D. E. J. Am. Chem. Soc. 1994, 116, 9787. (34) Cramer, C. J.; Worthington, S. E. J. Phys. Chem. 1995, 99, 1462. (35) Cramer, C. J.; Dulles, F. J.; Giesen, D. J.; Almlo¨f, J. Chem. Phys. Lett. 1995, 245, 165. (36) Lim, M. H.; Worthington, S. E.; Dulles, F. J.; Cramer, C. J. In Density-Functional Methods in Chemistry; Laird, B. B., Ziegler, T., Ross, R., Eds.; American Chemical Society: Washington, DC; in press. (37) Worthington, S. E.; Cramer, C. J. To be published. (38) Ziegler, T.; Rauk, A.; Baerends, E. J. Theor. Chim. Acta 1977, 43, 261. (39) Andersson, K.; Malmqvist, P.-A° .; Roos, B. O.; Sadlej, A. J.; Wolinski, K. J. Phys. Chem. 1990, 94, 5483. (40) Andersson, K.; Malmqvist, P.-A° .; Roos, B. O. J. Chem. Phys. 1992, 96, 1218. (41) Andersson, K. Theor. Chim. Acta 1995, 91, 31. (42) Hrovat, D. A.; Borden, W. T. To be published.
9670 J. Phys. Chem., Vol. 100, No. 23, 1996 (43) Hrovat, D. A.; Borden, W. T. J. Am. Chem. Soc. 1994, 116, 6327. (44) Yamanaka, S.; Kawakami, T.; Okumura, M.; Yamaguchi, K. Chem. Phys. Lett. 1995, 233, 257. (45) Becke, A. J. Chem. Phys. 1986, 84, 4524. (46) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (47) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (48) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (49) Andersson, K.; Blomberg, M. R. A.; Fu¨lscher, M. P.; Karlstro¨m, G.; Kello¨, V.; Lindh, R.; Malmqvist, P.-A° .; Noga, J.; Olsen, J.; Roos, B. O.; Sadlej, A. J.; Siegbahn, P. E. M.; Urban, M.; Widmark, P.-O. MOLCAS3; University of Lund: Lund, Sweden, 1994. (50) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347. (51) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Wong, M. W.; Foresman, J. B.; Robb, M. A.; Head-Gordon, M.; Replogle, E. S.; Gomperts, R.; Andres, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian 92/DFT, ReVision G.1; Gaussian, Inc.: Pittsburgh, PA, 1993. (52) ADF, Version 1.1.4, Department of Theoretical Chemistry, Vrije Universiteit, Amsterdam, 1993. (53) Snijders, G. J.; Baerends, E. J.; Vernooijs, P. At. Data Nucl. Data Tables 1982, 26, 483.
Cramer and Smith (54) Pople, J. A.; Scott, A. P.; Wong, M. W.; Radom, L. Isr. J. Chem. 1993, 33, 345. (55) R. R. Squires (personal communication) has found that a scale factor of 0.88 must be applied to the ZPVE of benzene calculated at the MCSCF(6,6)/3-21G level to achieve satisfactory comparison with experiment. (56) McKee, M. L. J. Phys. Chem. 1986, 90, 4908. (57) Andersson, K.; Roos, B. O. Int. J. Quantum Chem. 1993, 45, 591. (58) Jones, R. O. J. Chem. Phys. 1985, 82, 325. (59) Murray, C. W.; Handy, N. C.; Amos, R. D. J. Chem. Phys. 1993, 98, 7145. (60) This situation is by no means unique to the BP86 functional. We investigated several modern DFT functionals, and all gave identical results for the TMM multiplets and MCP to within 1 or 2 kcal/mol in relative energies. Inclusion of some (exact) Hartree-Fock exchange in the functional led to larger quantitative disagreements in relative energies (and also gave results in poorer agreement with those obtained at the CASPT2N level). (61) Hohenberg, P.; Kohn, W. Phys. ReV. B 1964, 3, 864. (62) Jordan, K. Personal communication. (63) Gordon, M. S. Personal communication. (64) Dowd, P. Tetrahedron Lett. 1991, 32, 445. (65) Eldin, S.; Liebman, J. F.; Reynolds, L. D.; Dowd, P. Tetrahedron Lett. 1992, 33, 4525.
JP953697X
