Article pubs.acs.org/JACS
Doubly Excited Character or Static Correlation of the Reference State in the Controversial 21Ag State of trans-Butadiene? Yinan Shu and Donald G. Truhlar* Department of Chemistry, Chemical Theory Center, and Minnesota Supercomputing Institute, University of Minnesota, Minneapolis, Minnesota 55455-0431, United States S Supporting Information *
ABSTRACT: Butadiene is the simplest polyene and has long served as a model system for many chemical and spectroscopic properties. However, this small molecule has presented significant challenges to theoretical chemistry. The 21Ag state, which is dark but photochemically important, is a prime source of this difficulty. Previous studies attributed the notorious difficulty in treating this state to strong double excitation character of the 21Ag state, which prevents the application of linear response (LR) methods. Therefore, one would require methods with much higher computational cost, especially for the analogues of this state in longer polyenes, and consequently studies of longer polyenes are very limited. In the present work, we argue that the difficulty stems more significantly from the inherently multiconfigurational character of the ground state. In addition, we validate the possibility of employing LR time-dependent density functional theory to investigate such a state with reasonable accuracy.
1. INTRODUCTION Experimental and theoretical studies of the photochemistry of all-trans polyenes have a long history dating back to the 1970s.1−10 This class of molecules is relevant to many areas of chemistry; for example, polyenes are model systems for carotenoids, which have important functions in photosynthesis,11 and the photodynamics of polyenes involves cis−trans isomerization, which serves as a representative example of light−mechanical energy conversion.12 These molecules provide the simplest examples of conjugated carbon bonds, but they provide significant challenges for theoretical chemistry, even for the shortest polyene, namely butadiene. There are several reasons: (1) Both of the low-lying excited states, 1Bu and 1Ag, of butadiene contribute significantly to its photodynamics,4,5,13 but the 1Ag excited state is a dark state, and hence only very weak absorption is observed experimentally. (2) The relative ordering of low-lying excited 1Bu and 1Ag states reverses when the length of the polyene changes.10,14−16 (3) The electronic structure of the 1Ag state is usually identified as an example of a state with significant double excitation character,3 which would mean that linear response methods cannot give the right answer for the right reason. (4) The 1Bu state includes some Rydberg character, while the 1Ag is known to be a valence state.17 The correlation effect for valence and Rydberg states are different; therefore it is not easy to find a balanced way to treat the two states. A recent study by Watson and Chan18 has summarized the widespread of the computed excitation energies the for S1 (11Bu) and S2 (21Ag) states of butadiene reported in the literature, and they found that chemical accuracy can be reached for the 21Ag state with the single-reference equation-of-motion coupled cluster (EOM© 2017 American Chemical Society
CC) method only when a high order of cluster operator and large basis sets are used. The electronic structure of the 11Bu state of butadiene is less complicated than that of the 21Ag state. The 11Bu state is dominated by a singly excited CSF (φHOMO → φLUMO), whereas the 21Ag state is conventionally represented by a mixing of singly and doubly excited CSFs (φHOMO−1 → φLUMO, φHOMO → φLUMO+1, and φHOMO2 → φLUMO2). This raises the question of whether linear response (LR) theory, for example, time-dependent Kohn−Sham density functional theory with full linear response (LR-TDDFT)19,20 or with the Tamm− Dancoff approximation (KS-TDA),21−24 is adequate for describing the 21Ag state. Note that to distinguish full LR and TDA, we denote the two methods as LR-TDDFT and KSTDA, where KS denotes Kohn−Sham. The matrix formalism of KS-TDA is similar to that of configuration interaction singles (CIS). Hence no doubly excited CSF is involved explicitly. On the other hand, it has been reported previously that KS-TDA predicts relatively accurate excitation energies with local functionals for the states that have similar character (1Bu state φHOMO → φLUMO and 1Ag state: φHOMO−1 → φLUMO, φHOMO → φLUMO+1). For example, KS-TDA/BLYP with the 6-31G** basis predicts vertical excitation energies of 6.53 and 6.56 eV for the 11Bu and 21Ag states,25 and KS-TDA/M06-L with 6311+G(2d,p) predicts 5.78 and 6.67 eV for the 11Bu and 21Ag states, respectively.26 For comparison, the experimental value is 5.92 eV for the 11Bu state,27 and a theoretical best estimate for the 21Ag state is 6.39 ± 0.07 eV.18 Received: June 16, 2017 Published: September 6, 2017 13770
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society Table 1. State-Averaged Natural Orbitals for the S0, S1 and S2 States Computed by SA-CASSCF and the Significantly Contributing CSFs with Their Corresponding Coefficients
As background we note that some molecules require a multiconfigurational wave function for a good description even of the ground state; such molecules are said to have high static correlation or to be strongly correlated. Yet another useful language28 is that they are said to have high multireference character. Here we are considering methods for calculating excitation energies, and in this context multireference methods are defined as those in which the orbitals are optimized with a multiconfigurational wave function, and single-reference methods are those with orbitals optimized with a single Slater determinant. In the current work, we compare the energies of low-lying excited states of butadiene computed using a multireference method with double excitations (MR) to those computed with single-reference methods that explicitly include only single excitations (SR-S) and a single-reference method that also explicitly includes double excitations (SR-D). By analyzing the excited 21Ag state of butadiene from an electron density point of view and comparing different methods in terms of energies, we examine the possibility of describing the 21Ag state with KS-TDA, which is an SR-S method. We show that the difficulty in describing the 21Ag state of butadiene may be better classified as a problem of static correlation of the reference state than as a consequence of needing double excitation character in the excited state, which is an unconventional view of this problem but which agrees with an analysis using other methods by Starcke et al.16
We performed equation-of-motion coupled cluster theory with single and double excitations (EOM-CCSD),35 state-specific complete active space self-consistent field36 (SS-CASSCF), state-averaged complete active space self-consistent field37,38 (SA-CASSCF), and multistate complete active space second-order perturbation theory39 (MS-CASPT2 based on SA-CASSCF orbitals) calculations by using Molpro.40 We generated state-specific natural orbitals (SSNOs) by diagonalizing the MS-CASPT2 first-order reduced density matrix, and we generated state-averaged natural orbitals (SANOs) from the SACASSCF first-order reduced density matrix. We also performed calculations by configuration interaction (CI) singles (S), by CIS with a doubles correction (CIS(D)),41,42 and by KS-TDA with functionals having various percentages of Hartree−Fock (HF) exchange by using the Q-Chem software package.43 The natural orbitals44 and effective number (nD or nU, depending on the method) of unpaired electrons45 were computed by using the MATROP function in Molpro. The natural transition orbitals46 of the KS-TDA calculations were computed by Q-Chem. All the single-point calculations are performed with the 6-31G** basis set because this basis has no diffuse functions and so avoids the problem that diffuse basis functions unphysically lower the Rydberg states in KS theory, which would cause them to couple unphysically to the valence states.16 (As mentioned above, the 2 1Ag state is known to be a valence state.) The effect of mixing with the Rydberg state will be investigated in section 5 by employing a larger basis set with diffuse functions, 6-311+G(2df,p). All density functional and wave function calculations in this paper are spin restricted, i.e., we do not employ different oribitals for different spins. Employing the broken-symmetry (unrestricted) formalism does not change the density functional results.
2. COMPUTATIONAL METHODS
3. THE WEIGHTED CONFIGURATION REPRESENTATION At the ground-state equilibrium geometry of butadiene, the MS-CASPT2 method with an active space of four electrons in four orbitals and mixing three states at the PT2 step predicts the first two singlet excited states to be at 6.69 and 6.48 eV for the 21Ag and 11Bu states, respectively. The excited states can be analyzed using the coefficients of the CSFs based on either the SANOs or the SSNOs; we used both methods. We label the key SANOs as HOMO−1, HOMO, LUMO, and LUMO+1. The assignments are made, as usual, such that
The ground-state geometry of butadiene is obtained from experiment, which yields a structure with C2h symmetry.29 We performed KS-TDA and dual-functional TDA (DF-TDA) calculations with the GAMESS and GAMESS+DF software packages.30−33 The DF-TDA method34 is an extension of KS-TDA that incorporates the capability of describing conical intersections of ground and excited states into KS-TDA. DF-TDA employs different exchange-correlation functionals for orbital optimization and for construction of the kernel used to evaluate excitations. The two functionals involved are called F1 and F2, respectively, with the notation DF-TDA/F1:F2. 13771
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
Table 2. Excitation Energies and Weighted Configuration Representations (WCRs) of the 11Bu and 21Ag States As Calculated by MS-CASPT2, EOM-CCSD, CIS(D), CIS, KS-TDA/M06, and DF-TDA/MN15:M06 with the 6-31G** Basis Set 21Ag
11Bu
method
type
orbitals for WCR
energy (eV)
WCR
energy (eV)
WCR
MS-CASPT2 EOM-CCSD CIS(D) CIS KS-TDA/M06 KS-TDA/M06-L DF-TDA/MN15:M06 Best Estimateb Experiment
MR SR-D SR-S(D) SR-S SR-S SR-S SR-S
SANO HF HF HF KS/M06 KS/M06-L KS/MN15
6.69 7.69 9.44 9.19 7.23 6.90 7.28 6.39 ± 0.07
0.93 1.13 0.34 h1.43 l l+1 −1 h 1.34 0.25 0.77 0.37 h−1 h l l+1 N/A 0.54 0.62 0.31 h1.24 l l+1 −1 h 0.53 0.53 0.47 h1.47 h l l+1 −1 0.53 0.54 0.46 h1.46 l l+1 −1 h 1.52 0.52 0.48 h1.48 l l+1 −1 h N/A
6.48 7.10 6.93 6.77 6.42 6.64 6.47 6.21 ± 0.02 5.92
0.94 0.94 0.00 h1.87 l l+1 −1 h 1.82 0.91 0.91 0.00 h−1 h l l+1 N/A 0.97 0.94 0.01 h1.90 l l+1 −1 h 0.95 0.94 0.01 h1.90 h l l+1 −1 0.95 0.91 0.04 h1.89 l l+1 −1 h 0.95 0.94 0.01 h1.90 l l+1 −1 h N/A
a
The MS-CASPT2 WCR is computed using SANOs from the SA-CASSCF calcualtions, the EOM-CISD, CIS(D), and CIS WCRs are based on Hartree−Fock orbitals, and the density functional WCRs are based on Kohn−Sham orbitals. bRef 18.
a double excitation, that perhaps this state will simply be missing in in a single-excitation theory, although it is known empirically48 that linear response TDDFT can give reasonably accurate excitation energies for doubly excited states. Our goal is to explain this apparent paradox. On general principles, it is clear that double excitation character contributes to low-lying states because doubly excited CSFs are coupled directly to both the ground state and singly excited CSFs. The MS-CASPT2 excited state wave function for 0.0 2.0 0.0 the 21Ag state is dominated by three CSFs, namely h2.0 −1h l l+1 , 1.0 2.0 1.0 0 2.0 1.0 0 1.0 h−1h l l+1, and h−1h l l+1 . Two of these CSFs are also present in the CIS, KS-TDA, and DF-TDA methods. The energies and WCRs computed by CIS, KS-TDA/M06,49 KSTDA/M06-L, and DF-TDA/MN15:M0650 are shown in Table 2. The KS/M06 and KS/M06-L orbitals employed in the WCR are very similar to the SANOs and HF orbitals. The CIS, KSTDA/M06, and DF-TDA/MN15:M06 calculations predict the 21Ag state at 9.19, 7.23, and 7.28 eV respectively. Although not fully summarized in Table 2, we found that when employing local functionals (i.e., those without HF exchange and hence with less static correlation error51 for the ground state), the computed results for the 21Ag state are greatly improved. For example, KS-TDA/M06-L, KS-TDA/BLYP52,53 and KS-TDA/ revPBE61 predict 6.90, 6.53, and 6.58 eV for the 21Ag state, respectively. By comparing the WCRs for the SR-D and SR-S methods that we have discussed above, we see that they are reasonably similar to each other. This is a key result of the present analysis. Since the EOM-CCSD excited state explicitly contains doubly excited CSFs, the 21Ag state computed at EOM-CCSD theoretical level should have a correct state character. Hence, since the HF and KS orbitals employed in the WCR are similar, and the EOM-CCSD and KS-TDA WCRs are similar, it is reasonable to argue that KS-TDA has the capability of describing the 21Ag state of butadiene with a qualitatively correct state character. The next question that arises naturally is why the 21Ag state predicted by KS-TDA has the correct character, although there are no doubly excited CSFs in this method. In fact, the doubly excited CSFs contribute significantly to the ground state. Consider a two-electron, two-orbital picture (let us label the two orbitals as a and b) where the ground state has a large coefficient for a2b0 and a small but significant coefficient for a0b2. The doubly excited state can be regarded as an alternative weighting of these two CSFs, a large coefficient of a0b2 and a small but significant coefficient of a2b0. In density functional
the HOMO−1, HOMO, LUMO, and LUMO+1 have the same character compared as the corresponding Hartree−Fock ground-state HOMO−1, HOMO, LUMO, and LUMO+1 orbitals. Table 1 shows the SANOs of butadiene and the coefficients of the CSFs in each of the states. For example, the ground state is dominated by a CSF with occupations 2.0 0 0 h2.0 −1h l l+1 and with coefficient 0.9. The contribution of this CSF to the ground-state wave function can be represented as 1.8 0 0 h1.8 −1h l l+1, where 1.8 is the product of the CSF coefficient and the natural orbital occupation. We call this the weighted configuration representation (WCR) in the following text. The total WCR is calculated by summing the contributions of all CSFs. In this way, calculating with SANOs, we find that the 0.93 1.13 0.34 l l+1 . The MS-CASPT2 WCR for the 21Ag state is h1.43 −1 h energies and the MS-CASPT2 WCRs of the 21Ag and 11Bu states are given in the first row of Table 2. We next consider the 21Ag state excitation energies computed by single-reference methods that include doubly excited CSFs in the wave function; in particular, we consider EOM-CCSD, which explicitly includes both single and double excitations, and CIS(D), which perturbatively includes the effect of double excitations on single excitations. Certainly the first of these methods and perhaps also the second should be able to describe states with double excitation character with reasonable accuracy. The energies and WCRs computed at both theoretical levels are shown in Table 2. EOM-CCSD/6-31G** predicts the 21Ag state at 7.69 eV, which is 1.30 eV higher than the best estimate (which is 6.39 ± 0.07 eV18). We know that this large error is not an effect of too small a basis set because it was previously shown18 that EOM-CCSD/cc-pV6Z predicts the 21Ag state at 7.25 eV. Table 2 shows that the WCR computed 1.25 0.77 0.37 l l+1 . The failure of by EOM-CCSD with 6-31G** is h1.34 −1 h EOM-CCSD was previously attributed to the fact that states with doubly excited CSFs require dynamic correlations from higher orders of excitations.18,47 However, inclusion the full triples in EOM-CC predicts the 21Ag state at 6.83 eV18, which is still not as accurate as MS-CASPT2. Table 2 shows that CIS(D)/6-31G** predicts the 21Ag state to be at 9.44 eV, which differs from the best estimate by 3.05 eV. Thus, neither the SR-D method nor the SR-S(D) method predicts an accurate excitation energy for the 21Ag state state. Next we consider calculations of the 21Ag state by the single reference methods without doubly excited CSFs, in particular CIS and KS-TDA. The first question that arises, given the history of the discussions of this state, is whether these methods can describe the 21Ag state at all. One can imagine, if it really is 13772
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
Table 3. State Specific Natural Orbitals for 21Ag State Computed at MS-CASPT2 and EOM-CCSD Theoretical Levels with Occupation Numbers That Strongly Deviate from 0 or 2a
a
In addition, the KS-TDA/M06 natural transition orbitals is shown in the table as well.
SSNO. This difference is also reflected in the big difference of the contribution to the WCR from the LUMO orbital. What causes this difference? One a priori difference between the MRD and SR-S/SR-D methods is the reference state. We pursue this issue in the next section.
theory, the ground state is already correlated despite being represented by a single determinant, and the doubly excited CSF information is implicitly contained in the determinant. This is not true for the CIS wave function, and the CIS computed energy has a much larger error (compared to the best estimates) for the 21Ag state than does KS-TDA. Compared to the best estimates, for the 21Ag state, the CIS error is 2.84 eV, the KS-TDA/M06-L error is 0.51 eV, and the KS-TDA/M06 error is 0.84 eV. As will be discussed in the next section, the larger error for M06 as compared to M06-L may be explained by the static correlation error51 brought in by the 27% Hartree−Fock exchange in the M06 functional. To further support the possibility of employing KS-TDA to describe the 21Ag state, we examine the state-specific natural orbitals (SSNOs) and the natural transition orbitals (NTOs). Table 3 shows the SSNOs for the 21Ag state computed by MSCASPT2 and EOM-CCSD, and the NTOs computed by KSTDA/M06. The similarities of EOM-CCSD SSNOs and the KS-TDA/M06 NTOs are clearly evident. It is important to note here that we are using the same isosurface values (0.02 and −0.02) for the orbital plots. From the above analysis, we find it is possible to employ KSTDA to describe the 21Ag state. However, there exists a clear difference between WCR derived from SR-S or SR-D and the WCR derived from MR. The WCR difference between MSCASPT2 and EOM-CCSD or KS-TDA is illustrated by the differences of the SSNOs and NTOs for the corresponding methods in Table 3. Especially for the LUMO orbital, we can observe a local π orbital is forming for the middle two carbons in the EOM-CCSD SSNO and KS-TDA NTO, but this local LUMO π orbital character is much weaker for the MS-CASPT2
4. THE MULTIREFERENCE CHARACTER OF THE REFERENCE STATE Table 1 shows that the ground state of butadiene has only a 90% contribution (i.e., 0.95 squared) from the dominant CSF, 2.0 0 0 namely h2.0 −1h l l+1. Hence it is strongly correlated. Another way to see this is by computing the number of effectively unpaired electrons for the ground state. The definition of the total number of effectively unpaired electrons is not unique, but it can be a useful tool for analysis.54 Because it is not unique, we will employ two variants. One definition is nD = Tr(D)
(1)
where the matrix D is a difference matrix,55,56 defined as D = 2P − P2
(2)
where P is the spinless one-particle reduced density matrix. This can be rewritten as N
nD =
∑
[2ni − ni2] (3)
(i = 1)
Another (more robust) definition is45 13773
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
Figure 1. Errors of the 11Bu and 21Ag states as a function of the percentage of the HF exchange employed in a functional; the remaining percentage of exchange and all the correlation taken from (a) the revPBE functional and (b) the BLYP functional. (c) The errors of the 11Bu and 21Ag states as functions of the percentage of HF exchange employed in the M06 family functionals. N
nU =
∑
diagnostic value greater than 0.10 represents large multireference character in the ground state.57 Martin and coworkers have shown that the M diagnostic, which is a direct and explicit measure of multireference character in the wave function, correlates very well with %TAE[T4+T5], which is a much more expensive but indirect diagnostic equal to the percentage of the molecular total atomization energy accounted for by connected quadruple and quintuple excitations; their analysis is consistent with these boundaries (0.05 and 0.10).58 For butadiene the computed value of the M diagnostic in the present work is M = 0.143. Therefore, the ground state has strong multireference character. To demonstrate the suitability of the M value as a trustworthy multireference diagnostic, we investigated the molecular size dependence of the M diagnostic. Here we compare the ground-state M values of three molecules, butadiene, ethylene, and ethylene dimer. (Both ethylene and ethylene dimer were optimized by KS-DFT with the M06-2X49 functional and the 6-311+G(2df,2p) basis set. The SSNOs are computed by the MS-CASPT2 method, averaging over two states with two active electrons in two active orbitals for ethylene and averaging over six states with four active electrons in four active orbitals for ethylene dimer.) The ethylene molecule was chosen for this illustration because ethylene can be treated as a “subunit” of ethylene dimer, and ethylene dimer was chosen because it has a molecular size similar to butadiene and the same number of π electrons as butadiene. Comparing the M value of butadiene to that of ethylene dimer is therefore another way to gauge the multireference character of the ground state of butadiene. The M values for ethylene and ethylene dimer are 0.091 and 0.097 respectively, whereas the M
[1 − |1 − ni|] (4)
(i = 1)
where N is the number of orbitals and ni is a natural orbital occupation number. If a state is dominated by a single closed-shell CSF, these numbers should be close to 0. The computed nD at MSCASPT2 theoretical level for ground-state butadiene is 1.62, and the computed nU for ground-state butadiene is 0.84. The effective numbers of unpaired electrons per double bond are nD = 0.81 and nU = 0.42, which deviate from 0 and indicate strong multireference character of the ground state. We further analyze the multireference character of the ground state and excited states by employing the M diagnostic analysis. The M diagnostic for a given state i is defined as57 1 M = [2 − ni(MCDONO) + 2 + ni(MCUNO)]
niSOMO
∑
|ni(j) − 1|
j
(5)
where ni(j) is the natural orbital occupation number of orbital j in state i, MCDOMO is the most correlated doubly occupied natural orbital in the dominant configuration for state i, MCUNO is the most correlated unoccupied orbital in the dominant configuration, and nSOMO is the number of singly i occupied molecular orbitals. The multireference character of a system is classified as small (M < 0.05), modest (0.05 ≤ M ≤ 0.10), or large (M > 0.10), depending on M. Previous studies indicate that an M diagnostic value greater than 0.05 represents modest multireference character in the ground state, and an M 13774
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
than the basis set used above. The EOM-CCSD method predicts the 21Ag state to be lowered to 7.18 eV above the ground state (which is still too high but now by “only” 0.79 eV); however, the state now has about 80% Rydberg character. The CIS and CIS(D) methods predict the 21Ag state to be at 7.48 and 7.49 eV, respectively; these excitation energies are much closer to the best estimate than those computed with the 6-31G** basis set (see Table 2 again), and the states are highly mixed with Rydberg character. The KS-TDA with M06-L, M06, M06-2X, revPBE and BLYP methods with the larger basis set predict the 21Ag state to be at 6.57, 6.34, 7.01, 6.14, and 6.07 eV respectively. The computed 21Ag states with the local functionals (M06-L, revPBE and BLYP) have a weaker Rydberg character compared to what is predicted by hybrid functionals. A comparison of the computed 21Ag state accuracy between 631G** and 6-311+G(2df,p) basis sets is shown in Figure 2 for
value of butadiene is 47% larger than that of ethylene dimer. The comparison of the M values of ethylene and ethylene dimer shows that the M value does not have a strong molecular size dependence, and it further supports our conclusion that the ground state of butadiene has large multiconfigurational character. The inclusion of ethylene dimer in this comparison is important because it shows that the M diagnostic does not have the undesirable feature of another widely used diagnostic, i.e., the percentage weight of the dominant CSF (which is given by the square of the leading coefficient in the CI expansion). This weight is 96.6% for ethylene, 92.4%, for ethylene dimer, and 90.3% for butadiene, and comparison of the values for ethylene (3.4% less than 100%) and ethylene dimer (7.6% less than 100%) demonstrates that the dominant-weight diagnostic does have a significant dependence on system size, which makes it unsatisfactory for comparing the multiconfigurational characters of different molecules. The above results show that the reference state of the EOMCC and KS-TDA methods may not have captured enough static correlation for the ground state; this affects the excitation energy and the excited state WCR in that the wrong groundstate description leads to an incorrect response. We conclude that this is the reason why MS-CASPT2 can effectively describe the 21Ag state of butadiene, but single-reference methods cannot. This conclusion is further supported by the fact that the KS-TDA computed 21Ag state is much more accurate when employing a local exchange-correlation functional than a hybrid functional.59 Figure 1(a) shows the errors of the 11Bu and 21Ag states as a function of the percentage of HF exchange in the exchange-correlation functional, where the rest of the exchange and all the correlation are from the revPBE60,61 functional. The errors are computed as the difference between the computed excitation energies and the best estimates. The accuracy of the 21Ag state is an almost linear function of the percentage of HF exchange. In contrast, the accuracy of 11Bu state (although it does increase when the percentage of the HF exchange increases) does not change much when different amounts of HF exchange are employed. We found similar behavior for the other functionals, for example, the results for exchange-adjusted BLYP functionals are plotted in Figure 1(b) and those for the M06 family of functionals49 as plotted in Figure 1(c). The error of the 21Ag state increases significantly as a function of the percentage of HF exchange. Note that in the M06 family, the parameters were reoptimized for each percentage of HF exchange; hence the errors are relatively reduced when a higher percentage of the HF exchange is employed, but the trend is still apparent. We interpret this persistent trend as caused by the fact that Hartree−Fock exchange increases the static correlation error51 of the ground state. The above trends for the 11Bu state vis-à-vis the 21Ag state are attributed to the fact that 11Bu state has a much weaker multireference character than the other states considered here. The computed M diagnostic value for the 11Bu state is 0.056, which may be compared to 0.143 for the ground state and 0.416 for the 21Ag state. Hence, unlike the case for the 21Ag state, the inclusion of the ground state’s multireference character is not necessary to build the multireference character into the description of the 11Bu state by EOM or LR.
Figure 2. Errors of 21Ag state computed at EOM-CCSD, CIS(D), CIS, KS-TDA/M06, DF-TDA/MN15:M06, KS-TDA/M06-L, KS-TDA/ revPBE, KS-TDA/BLYP, DF-TDA/M06-L:BLYP methods with 631G** and 6-311+G(2df,p) basis sets. MS-CASPT2 with 6-31G** basis set computed error is shown at the first column for comparison.
the above-mentioned single-reference methods. In addition, the MS-CASPT2 error computed with the 6-31G** basis set is also shown in Figure 2 for comparison. By comparing these excitation energies computed with two different basis sets within the same method, we find using a large diffuse basis set gives much better accuracy compared with a smaller basis set, although with a wrong state character (highly mixed with Rydberg character). Hence the improved accuracy when using large diffuse basis set is most likely an artifact from error cancellation since density functional theory with standard functionals tends to underestimate62−64 Rydberg state energies.
5. THE EFFECT OF RYDBERG CHARACTER The mixing of Rydberg character into the excited states is investigated by performing single-point calculations at the same geometry with 6-311+G(2df,p) basis set, which is more diffuse 13775
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
6. DISCUSSION The current work focuses on the 21Ag state of butadiene. Previous studies usually assume that the difficulty of achieving an accurate description of the 21Ag state comes from the strong mixing of double excitation character into the 21Ag state, and it is clear from the MS-CASPT2 and EOM-CC calculations that the double excitation character is very strong. However, the question then arises of why MS-CASPT2 provides a much more accurate way to describe such a state than does EOM-CC. It is also interesting that TDDFT with local functionals can achieve much better accuracy compared with EOM-CC methods, despite the fact that the doubly excited CSF is not in the KS-TDA expansion at all. Hence we focused on understanding the following two questions: (1) Is it possible to use TDDFT to get accurate results for the 21Ag state of butadiene? (2) If it is possible, why does it work? To analyze these we employed two approaches: one, called WCR, developed in the present work, and another called number of unpaired electrons, developed in earlier literature. The WCR analysis compares the averaged electron occupations for the most important four valence orbitals among five electronic structure methods (MS-CASPT2, EOM-CCSD, CIS, KS-TDA, and DF-TDA). Doing this, we found the WCR of KS-TDA and DF-TDA are very close to that of EOM-CC. Since we know that the EOM-CC expansion contains the doubly excited CSF, this shows that a single-excitation method like KS-TDA can effectively represent the same kind of character with a linear combination of singly exctied configurations. Also, the following facts should be emphasized: (1) The low-lying 21Ag excited state for butadiene contains strong double excitation character because the doubly excited CSF φHOMO2 → φLUMO2, mixes strongly with singly excited CSFs φHOMO−1 → φLUMO and φHOMO → φLUMO+1. Hence the 21Ag state is not completely missing in the KS-TDA expansion; although the doubly excited CSFs are missing. (2) In density functional theory, the double excitation information is partly buried in the ground state, because a large fraction of the dynamic correlation is captured by the exchange-correlation functional. (3) There exist functionals, for example, BLYP, revPBE, and M06-L, that can achieve high accuracy for the 21Ag state of butadiene in terms of energies from a LR-TDDFT or KS-TDA calculation. With these facts, it is reasonable to argue that LR-TDDFT or KS-TDA should be able to describe the low-lying excited 21Ag state of butadiene. So what is the real difficulty of computing the 21Ag state? We believe that this is more likely that the ground state of the butadiene has a strong multireference character. This is supported by the following facts: (1) The MS-CASPT2 method describes the 21Ag state much more efficiently than EOM-CC. (2) The effective number of unpaired electrons for the ground state of butadiene is nD = 1.62 or nU = 0.84 depend on the definition; and the M diagnostic value is 0.143. (3) Local functionals employed in TDDFT behaves much better than hybrid functionals for the 21Ag state, and there exists a strong positive correlation between the error of the computed 21Ag state and the percentage of HF exchange employed in the functionals. One may extend the previous observation to conclude that the multireference character of the ground state promotes the contribution of double excited configurations in low-lying excited states.
7. CONCLUSIONS In the current work, we have investigated the low-lying excited 21Ag state of butadiene by a combination of methods. We validated the possibility of employing TDDFT for the 21Ag state, which has a strong nominal double excitation character. In addition, we found that the difficulty of treating the 21Ag state comes from the fact that the ground state of butadiene has strong multireference character. The 21Ag state is known to have valence character, but employing a diffuse basis set lowers the Rydberg states and causes artificial coupling between the valence and Rydberg states. The unphysical mixing between the 21Ag state and a Rydberg state causes KS-TDA and EOMCCSD to predict the excitation energies of the 21Ag state more accurately but with a wrong character. Hence we suggest the use of local functionals and a relative small basis set to study the 21Ag state of butadiene. Although further investigation is required, these conclusions may extend to the whole class of polyene molecules. Since the early days of quantum mechanics, the treatment of conjugated systems of π electrons has played a prominent role in the quantum theory of electronic structure in chemistry, starting with Hückel65 and continuing through the work of Pariser, Parr, Pople, Jug, and others.66,67 Two paradigmshattering events rocked this field though. One was the discovery of the large role of σ electrons in aromaticity.68−70 The other was the discovery of the important role of doubly excited and Rydberg configurations in the low-lying states of even the simplest conjugated system, namely butadiene.71 The role of doubly excited configurations continues to be controversial. In the present article we have argued that the problem cannot be understood unless one uses a multiconfigurational treatment of the ground state. The presence of doubly excited configurations in the ground state (sometimes labeled as static correlation, strong correlation, or multireference character) is inevitably accompanied by a different mixture of these configurations in low-lying excited states. The progress of quantum chemical understanding has now reached the stage where inherently multiconfigurational character in ground states must play a greater role not just in reference states of quantitative treatments but also in our zero-order orbital pictures of electronic excitation.
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/jacs.7b06283.
■
Equilibrium geometry of butadiene used for the calculations (PDF)
AUTHOR INFORMATION
Corresponding Author
*
[email protected] ORCID
Yinan Shu: 0000-0002-8371-0221 Donald G. Truhlar: 0000-0002-7742-7294 Notes
The authors declare no competing financial interest. 13776
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society
■
General Atomic and Molecular Electronic Structure System (GAMESS)). (33) https://comp.chem.umn.edu/gamess+df. (34) Shu, Y.; Parker, K. A.; Truhlar, D. G. J. Phys. Chem. Lett. 2017, 8, 2107−2112. (35) Stanton, J. F.; Bartlett, R. J. J. Chem. Phys. 1993, 98, 7029−7039. (36) Roos, B. O.; Taylor, P. R. Chem. Phys. 1980, 48, 157. (37) Werner, H. J.; Knowles, P. J. J. Chem. Phys. 1985, 82, 5053. (38) Li, Y.; Francisco, J. S. Spectrochim. Acta, Part A 1999, 55, 477− 485. (39) Finley, J.; Malmqvist, P. A.; Roos, B. O.; Serrano-Andres, L. Chem. Phys. Lett. 1998, 288, 288−299. (40) Werner, H.-J.; Knowles, P. J.; Knizia, G.; Manby, F. R.; Schütz, M.; Celani, P.; Korona, T.; Lindh, R.; Mitrushenkov, A.; Rauhut, G.; Shamasundar, K. R.; Adler, T. B.; Amos, R. D.; Bernhardsson, A.; Berning, A.; Cooper, D. L.; Deegan, M. J. O.; Dobbyn, A. J.; Eckert, F.; Goll, E.; Hampel, C.; Hesselmann, A.; Hetzer, G.; Hrenar, T.; Jansen, G.; K öppl, C.; Liu, Y.; Lloyd, A. W.; Mata, R. A.; May, A. J.; McNicholas, S. J.; Meyer, W.; Mura, M. E.; Nicklass, A.; O’Neill, D. P.; Palmieri, P.; Pflüger, K.; Pitzer, R.; Reiher, M.; Shiozaki, T.; Stoll, H.; Stone, A. J.; Tarroni, R.; Thorsteinsson, T.; Wang, M.; Wolf, A. Molpro, version 2010.1, a package of ab initio programs; see http:// www.molpro.net. (41) Head-Gordon, M.; Rico, R. J.; Oumi, M.; Lee, T. J. Chem. Phys. Lett. 1994, 219, 21−29. (42) Head-Gordon, M.; Maurice, D.; Oumi, M. Chem. Phys. Lett. 1995, 246, 114−121. (43) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T. B.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; Khaliullin, T. K.; Landau, A.; Liu, J.; Proynov, E. I.; Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J.; Woodcock, H. L.; III; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chan, C. M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M.; Distasio, R. A., Jr.; Do, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; Fusti-Molnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W. D.; Harbach, P. H. P.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T. C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S. P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Oana, C. M.; OlivaresAmaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Sharada, S. M.; Sharma, S.; Small, D. W.; Sodt, A.; Stein, T.; Stück, D.; Su, Y. C.; Thom, A. J. W.; Tsuchimochi, T.; Vanovschi, V.; Vogt, L.; Vydrov, O.; Watson, M. A.; Wenzel, J.; White, A.; Williams, C. F.; Yeganeh, S.; Yost, S. R.; You, Z. Q.; Zhang, I. Y.; Zhang, X.; Zhao, Y.; Brooks, B. R.; Chan, G. K. L.; Chipman, D. M.; Cramer, C. J.; Goddard, W. A.; III; Gordon, M. S.; Hehre, W. J.; Klamt, A.; Schaefer, H. F.; III; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xu, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J. D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C. P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W. Z.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Voorhis, T. V.; Herbert, J. M.; Krylov, A. I.; Gill, P. M. W.; Head-Gordon, M. Mol. Phys. 2015, 113, 184−215. (44) Löwdin, P.-O.; Shull, H. Phys. Rev. 1956, 101, 1730−1739. (45) Head-Gordon, M. Chem. Phys. Lett. 2003, 372, 508−511. (46) Martin, R. L. J. Chem. Phys. 2003, 118, 4775−4777. (47) Saha, B.; Ehara, M. J. Chem. Phys. 2006, 125, 014316. (48) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 302, 375− 382. (49) Zhao, Y.; Truhlar, D. G. Theor. Chem. Acc. 2008, 120, 215−241. (50) Yu, H. S.; He, X.; Li, S. L.; Truhlar, D. G. Chem. Sci. 2016, 7, 5032−5051.
ACKNOWLEDGMENTS This work was supported in part by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0015997.
■
REFERENCES
(1) Buenker, R. J.; Whitten, J. L. J. Chem. Phys. 1968, 49, 5381−5387. (2) Shih, S.; Buenker, R. J.; Peyerimhoff, S. D. Chem. Phys. Lett. 1972, 16, 244−251. (3) Schulten, K.; Karplus, M. Chem. Phys. Lett. 1972, 14, 305−309. (4) Hudson, B. S.; Kohler, B. Chem. Phys. Lett. 1972, 14, 299−304. (5) Hudson, B. S.; Kohler, B. J. Chem. Phys. 1973, 59, 4984−5002. (6) Mosher, O. A.; Flicker, W. M.; Kuppermann, A. Chem. Phys. Lett. 1973, 19, 332−333. (7) Mosher, O. A.; Flicker, W. M.; Kuppermann, A. J. Chem. Phys. 1973, 59, 6502−6511. (8) Hudson, B. S.; Kohler, B. Annu. Rev. Phys. Chem. 1974, 25, 437− 460. (9) Flicker, W.; Mosher, O. A.; Kuppermann, A. Chem. Phys. 1978, 30, 307−314. (10) Graham, R. L.; Freed, K. F. J. Chem. Phys. 1992, 96, 1304−1316. (11) Gust, D.; Moore, T. A.; Moore, A. L. Acc. Chem. Res. 2001, 34, 40−48. (12) Levine, B. G.; Martínez, T. J. Annu. Rev. Phys. Chem. 2007, 58, 613−634. (13) Levine, B. G.; Martínez, T. J. J. Phys. Chem. A 2009, 113, 12815−12824. (14) Cave, R. J.; Davidson, E. R. Chem. Phys. Lett. 1988, 148, 190− 196. (15) Cave, R. J.; Davidson, E. R. J. Phys. Chem. 1988, 92, 614−620. (16) Starcke, J. H.; Wormit, M.; Schirmer, J.; Dreuw, A. Chem. Phys. 2006, 329, 39−49. (17) Serrano-Andrés, L.; Merchán, M.; Nebot-Gill, I.; Lindh, R.; Roos, B. O. J. Chem. Phys. 1993, 98, 3151−3162. (18) Watson, M. A.; Chan, G. K.-L. J. Chem. Theory Comput. 2012, 8, 4013−4018. (19) Casida, M. E. In Recent Advances in Density Functional Methods, Part I; Chong, D. P., Ed.; World Scientific: Singapore, 1995; pp 155− 192. (20) Bauernschmitt, R.; Ahlrichs, R. Chem. Phys. Lett. 1996, 256, 454−464. (21) Hirata, S.; Head-Gordon, M. Chem. Phys. Lett. 1999, 314, 291− 299. (22) Fetter, A. L.; Walecka, J. D. Quantum Theory of Many-Particle Systems; Courier Corporation: New York, 2012. (23) Handy, N. C.; Tozer, D. J. J. Comput. Chem. 1999, 20, 106−113. (24) Hirata, S.; Lee, T. J.; Head-Gordon, M. J. Chem. Phys. 1999, 111, 8904−8912. (25) Hsu, C.-P.; Hirata, S.; Head-Gordon, M. J. Phys. Chem. A 2001, 105, 451−458. (26) Jacquemin, D.; Perpètem, E. A.; Ciofini, I.; Adamo, C.; Valero, R.; Zhao, Y.; Truhlar, D. G. J. Chem. Theory Comput. 2010, 6, 2071− 2085. (27) Doering, J. P.; McDiarmid, R. J. Chem. Phys. 1980, 73, 3617− 3624. (28) Truhlar, D. G. J. Comput. Chem. 2007, 28, 73−86. (29) Haugen, W.; Traetteberg, M. Acta Chem. Scand. 1966, 20, 1726−1726. (30) 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., Jr. J. Comput. Chem. 1993, 14, 1347−1363. (31) Gordon, M. S.; Michael, W. S. Theory and Applications of Computational Chemistry: The First Forty Years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 1167−1189. (32) Shu, Y.; Marenich, A. V.; Parker, K.; Truhlar, D. G. Gamess+DF; University of Minnesota: Minneapolis, MN, 2017 (based on the 13777
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778
Article
Journal of the American Chemical Society (51) Yu, H. S.; Li, S. L.; Truhlar, D. G. J. Chem. Phys. 2016, 145, 130901. (52) Becke, A. D. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098− 3100. (53) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789. (54) Staroverov, V. N.; Davidson, E. R. Int. J. Quantum Chem. 2000, 77, 316−323. (55) Takatsuka, K.; Fueno, T.; Yamaguchi, K. Theor. Chim. Acta 1978, 48, 175−183. (56) Lain, L.; Torre, A.; Bochicchio, R. C.; Ponec, R. Chem. Phys. Lett. 2001, 346, 283−287. (57) Tishchenko, O.; Zheng, J.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 1208−1219. (58) Fogueri, U. R.; Kozuch, S.; Karton, A.; Martin, J. M. L. Theor. Chem. Acc. 2013, 132, 1291. (59) Cave, R. J.; Zhang, F.; Maitra, N. T.; Burke, K. Chem. Phys. Lett. 2004, 389, 39−42. (60) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865−3868. (61) Zhang, Y.; Yang, W. Phys. Rev. Lett. 1998, 80, 890. (62) Casida, M. E.; Jamorski, C.; Casida, K. C.; Salahub, D. R. J. Chem. Phys. 1998, 108, 4439−4449. (63) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2006, 110, 13126− 13130. (64) Li, S. L.; Truhlar, D. G. J. Chem. Theory Comput. 2015, 11, 3123−3130. (65) Hückel, E. Eur. Phys. J. A 1931, 70, 204−286. (66) Pople, J. A. Trans. Faraday Soc. 1953, 49, 1375−1385. (67) Bredow, T.; Jug, K. Theor. Chem. Acc. 2005, 113, 1−14. (68) Shaik, S. S.; Hiberty, P. C. J. Am. Chem. Soc. 1985, 107, 3089− 3095. (69) Jug, K.; Köster, A. M. J. Am. Chem. Soc. 1990, 112, 6772−6777. (70) Hiberty, P. C.; Danovich, D.; Shurki, A.; Shaik, S. J. Am. Chem. Soc. 1995, 117, 7760−7769. (71) Hudson, B.; Kohler, B. Annu. Rev. Phys. Chem. 1974, 25, 437− 460.
13778
DOI: 10.1021/jacs.7b06283 J. Am. Chem. Soc. 2017, 139, 13770−13778