Despite the fact that excited states are involved in ubiquitious processes such as photochemistry, catalysis or in solar cell technology, none of the many methods existing is the reference in providing accurate excitation energies. Indeed, each method has its own flaws and there are so many chemical scenario that can occur, so it is still one of the biggest challenge in theoretical chemistry. Speaking of difficult task, cyclobutadiene (CBD) molecule has been a real challenge for experimental and theoretical chemists for many decades \cite{bally_1980}. Due to his antiaromaticity \cite{AromaticityAntiaromaticityElectronic,} and his large angular strain \cite{baeyer_1885} the CBD molecule presents a high reactivity which made the synthesis of this molecule a particularly difficult exercise. Hückel molecular orbital theory gives a triplet state with square ($D_{4h}$) geometry for the ground state of the CBD,with the two singly occupied frontier orbitals that are degenerated by symmetry. This degeneracy is lifted by the Jahn-Teller effect, meaning by distortion of the molecule (lowering symmetry), and gives a singlet state with rectangular ($D_{2h}$) geometry for the ground state.
Indeed, synthetic work from Pettis and co-workers \cite{reeves_1969} gives a rectangular geometry to the singlet ground state of CBD and then was confirmed by experimental works \cite{irngartinger_1983,ermer_1983,kreile_1986}.
At the ground state structrure ($D_{2h}$), the ${}^1A_g$ state has a weak multi-configurational character because of the well separated frontier orbitals and can be described by single-reference methods. But at the square ($D_{4h}$) geometry, the singlet state ${}^1B_{1g}$ has two singly occupied frontier orbitals that are degenerated so has a two-configurational character and single-reference methods are unreliable to describe it. The singlet ($D_{4h}$) is a transition state in the automerization reaction between the two rectangular structures (see Fig.\ref{fig:CBD}). The autoisomerization barrier for the CBD molecule is defined as the energy difference between the singlet ground state of the square ($D_{4h}$) structure and the singlet ground state of the rectangular ($D_{2h}$) geometry. The energy of this barrier was predicted, experimentally, in the range of 1.6-10 kcal.mol$^{-1}$\cite{whitman_1982} and multi-reference calculations gave an energy barrier in the range of 6-7 kcal.mol$^{-1}$\cite{eckert-maksic_2006}. All the specificities of the CBD molecule make it a real playground for excited-states methods.
Excited states of the CBD molecule in both geometries are represented in Fig.\ref{fig:CBD}. Are represented ${}^1A_g$ and $1{}^3B_{1g}$ states for the rectangular geometry and ${}^1B_{1g}$and $1{}^3A_{2g}$ for the square one. Due to energy scaling doubly excited states $1{}^1B_{1g}$ and $2{}^1A_{1g}$ for the $D_{2h}$ and $D_{4h}$ structures, respectively, are not drawn. Doubly excited states are known to be challenging to represent for adiabatic time-dependent density functional theory (TD-DFT) and even for state-of-the-art methods like the approximate third-order coupled-cluster (CC3) \cite{christiansen_1995,koch_1997} or equation-of-motion coupled-cluster with singles, doubles and triples (EOM-CCSDT) \cite{kucharski_1991,kallay_2004,hirata_2000,hirata_2004}.
In order to tackle the problems of multi-configurational character and double excitations several ways are explored. The most evident way that one can think about to describe multiconfigurational and double excitations are multiconfigurational methods. Among these methods, one can find complete active space self-consistent field (CASSCF) \cite{roos_1996}, the second perturbation-corrected variant (CASPT2) \cite{andersson_1990} and the second-order $n$-electron valence state perturbation theory (NEVPT2) \cite{angeli_2001b,angeli_2001a,angeli_2002}. The exponential scaling of these methods with the size of the active space is the limitation to the application of these ones to big molecules.
Another way to deal with double excitations is to use high level truncation of the equation-of-motion (EOM) formalism of coupled-cluster (CC) theory. However, to provide a correct description of doubly excited states one have to take into account contributions from the triple excitations in the CC expansion. Again, due to the scaling of CC methods with the number of basis functions the applicability of these methods is limited to small molecules.
An alternative to multiconfigurational and CC methods is the use of selected CI (SCI) methods for the computation of transition energies for singly and doubly excited states that are known to reach near full CI energies for small molecules. These methods allow to avoid an exponential increase of the size of the CI expansion by retaining the most energetically relevant determinants only, using a second-order energetic criterion to select perturbatively determinants in the FCI space.
Finally, to describe doubly excited states, one can think of spin-flip formalism established by Krylov in 2001. To briefly introduce the spin-flip idea we can present it like: instead of taking the singlet ground state as reference, the reference configuration is taken as the lowest triplet state. So one can access the singlet ground state and the singlet doubly-excited state via a spin-flip deexcitation and excitation (respectively), the difference of these two excitation energies providing an estimate of the double excitation. Obviously spin-flip methods have their own flaws, especially the spin contamination \cite{casanova_2020} (i.e., an artificial mixing of electronic states of differents spin multiplicities) due to spin incompleteness of the spin-flip expansion and by spin contamination of the reference configuration. One can adress part of this problem by expansion of the excitation order but with an increase of the computational cost or by complementing the spin-incomplete configuration set with the missing configurations.
In the present work we investigate ${}^1A_g$, $1{}^3B_{1g}$, $1{}^1B_{1g}$, $2{}^1A_{g}$ and ${}^1B_{1g}$, $1{}^3A_{2g}$, $2{}^1A_{1g}$,$1{}^1B_{2g}$ excited states for the $D_{2h}$ and $D_{4h}$ geometries, respectively. Computational details are reported in Section \ref{sec:compmet} for SCI (Subsection \ref{sec:SCI}), EOM-CC (Subsection \ref{sec:CC}), multiconfigurational (Subsection \ref{sec:Multi}) and spin-flip (Subsection \ref{sec:sf}) methods. Section \ref{sec:res} is devoted to the discussion of our results, first we consider the ground state property studied which is the autoisomerization barrier (Subsection \ref{sec:auto}) and then we study the excited states (Subsection \ref{sec:states}) of the CBD molecule for both geometries.
%The system under investigation in this work is the cyclobutadiene (CBD) molecule, rectangular ($D_{2h}$) and square ($D_{4h}$) geometries are considered. The ($D_{2h}$) geometry is obtained at the CC3 level without frozen core using the aug-cc-pVTZ and the ($D_{4h}$) geometry is obtained at the RO-CCSD(T) level using aug-cc-pVTZ again without frozen core. All the calculations are performed using four basis set, the 6-31+G(d) basis and the aug-cc-pVXZ with X$=$D, T, Q. In the following we will use the notation AVXZ for the aug-cc-pVXZ basis, again with X$=$D, T, Q. The $\%T_1$ metric that gives the percentage of single excitation calculated at the CC3 level in \textcolor{red}{DALTON} allows to characterize the various states.Throughout all this work, spin-flip and spin-conserved calculations are performed with a UHF reference.
States energies and excitations energies calculations in the SCI framework are performed using QUANTUM PACKAGE where the CIPSI algorithm is implemented. The CIPSI algorithm allows to avoid the exponential increase of th CI expansion. To treat electronic states in the same way we use a state-averaged formalism meaning that the ground and excited states are represented with the same number and same set of determinants but using different CI coefficients. Then the SCI energy is the sum of two terms, the variational energy obtained by diagonalization of the CI matrix in the reference space and a second-order perturbative correction which estimates the contribution of the determinants not included in the CI space (estimate error in the truncation). It is possible to estimate the FCI limit for the total energies and compute the corresponding transition energies by extrapolating this second-order correction to zero. Extrapolation brings error and we can estime this one by energy difference between excitation energies obtained with the largest SCI wave function and the FCI extrapolated value. These errors provide a rough idea of the quality of the FCI extrapolation and cannot be seen as true bar error, they are reported in the following tables.
Different flavours of coupled-cluster (CC) calculations are performed using different codes. Indeed, CC theory provides a hierarchy of methods that provide increasingly accurate energies via the increase of the maximum excitation degree of the cluster operator. Without any truncation of the cluster operator one has the full CC (FCC) that is equivalent to the full configuration interaction (FCI) giving the exact energy and wave function of the system for a fixed atomic basis set. However, due to the computational exponential scaling with system size we have to use truncated CC methods. The CC with singles and doubles (CCSD), CC with singles, doubles and triples (CCSDT) calculations are achieved with \textcolor{red}{CFOUR}. The calculations in the context of CC response theory or ``approximate'' series (CC3,CC4) are performed with \textcolor{red}{DALTON}. The CC with singles, doubles, triples and quadruples (CCSDTQ) are done with the \textcolor{red}{CFOUR} code. The CC2, CC3 and CC4 methods can be seen as cheaper approximations of CCSD, CCSDT and CCSDTQ by skipping the most expensive terms and avoiding the storage of higher-excitations amplitudes.
State-averaged complete-active-space self-consistent field (SA-CASSCF) calculations are performed with \textcolor{red}{MOLPRO}. On top of those, NEVPT2 calculations, both partially contracted (PC) and strongly contracted (SC) scheme are considered. The PC-NEVPT2 is theoretically more accurate to the strongly contracted version due to the larger number of pertubers and greater flexibility. CASPT2 is performed and extended multistate (XMS) CASPT2 for strong mixing between states with same spin and spatial symmetries is also performed.
In both structures the CBD has a singlet ground state, for the spin-flip calculations we consider the lowest triplet state as reference. Spin-flip techniques are broadly accessible and here, among them, we explore equation-of-motion coupled-cluster singles and doubles (EOM-CCSD), configuration interaction singles (CIS), algebraic-diagrammatic construction (ADC) scheme and TD-DFT. The standard and extended spin-flip ADC(2) (SF-ADC(2)-s and SF-ADC(2)-x respectively) and SF-ADC(3) are performed using Q-CHEM 5.2.1. Spin-flip TD-DFT calculations are also performed using Q-CHEM 5.2.1. The B3LYP, PBE0 and BH\&HLYP Hybrid GGA functionals are considered, they contain 20\%, 25\%, 50\% of exact exchange and are labeled, respectively, as SF-BLYP, SF-B3LYP, SF-PBE0, SF-BH\&HLYP. We also have done spin-flip TD-DFT calculations using Range-Separated Hybrid (RSH) functionals as: CAM-B3LYP, LC-$\omega$PBE08 and $\omega$B97X-V. The main difference here between these RSH functionals is the amount of exact-exchange at long-range: 75$\%$ for CAM-B3LYP and 100$\%$ for LC-$\omega$PBE08 and $\omega$B97X-V. To complete the use of spin-flip TD-DFT we also considered the Hybrid meta-GGA functional M06-2X and the RSH meta-GGA functional M11. EOM-SF-CCSD and EOM-SF-CC(2,3) are also performed with Q-CHEM 5.2.1.
Two different sets of geometries obtained with different level of theory are considered for the ground state property and for the excited states of the CBD molecule. First, for the autoisomerization barrier because we consider an energy difference between two geometries it is needed to obtain these geometries at the same level of theory. Due to the fact that the square CBD is an open-shell molecule it is difficult to optimize the geometry so the most accurate method that we can use for both structures is the CASPT2(12,12) with the aug-cc-pVTZ (AVTZ) basis without frozen core. Then, for the excited states because we look at vertical energy transitions in one particular geometry we can use different methods for the different structures and use the most accurate method for each geometry. So in the case of the excited states of the CBD molecule we use CC3 without frozen core with the aug-cc-pVTZ basis for the rectangular ($D_{2h}$) geometry and we use RO-CCSD(T) with the aug-cc-pVTZ (AVTZ) basis again without frozen core for the square ($D_{4h}$) geometry. Table \ref{tab:geometries} shows the results on the geometry parameters obtained with the different methods.
The results for the calculation of the automerization barrier energy with and without spin-flip methods are shown in Table \ref{tab:auto_standard}. Two types of methods are used with spin-flip, SF-TD-DFT with several functionals and SF-ADC with the variants SF-ADC(2)-s, SF-ADC(2)-x and SF-ADC(3). First, one can see that there is large variations of the energy between the different functionals. Indeed if we look at the hybrid functionals B3LYP and BH\&HLYP we can see that there is a difference of around 7 \kcalmol throught all the bases, the difference in energy between the B3LYP and PBE0 functionals is much smaller with around 1.5 \kcalmol throught all the bases. We find a similar behaviour regarding the RSH functionals, we find a difference of about 8.5-9 \kcalmol between the M06-2X and the CAM-B3LYP functionals for all bases. The results between the CAM-B3LYP and the $\omega$B97X-V are very close in energy with a difference of around 0.15-0.25 \kcalmol . The energy difference between the M11 and the M06-2X functionals is larger with 0.7-0.8 \kcalmol for the AVXZ bases and with 1.79 \kcalmol for the 6-31+G(d) basis. For the SF-ADC methods the energy differences are smaller with 1.76-2 \kcalmol between the ADC(2)-s and the ADC(2)-x schemes, 0.94-1.61 \kcalmol between the ADC(2)-s and the ADC(3) schemes and 0.39-0.82 \kcalmol between the ADC(2)-x and the ADC(3) schemes.
Then we compare results for multireference methods, we can see a difference of about 2.91-3.22 \kcalmol through all the bases between the CASSCF(12,12) and the CASPT2(12,12) methods. These differences can be explained by the well known lack of dynamical correlation for the CASSCF method that is typically important for close-shell molecules which is the case for the ground states of the rectangular and square geometries of the CBD molecule. The results between the CASPT2(12,12) and the NEVPT2(12,12) are much closer with an energy difference of around 0.12-0.23 \kcalmol for all the bases. Finally the last results shown in Table \ref{tab:auto_standard} are the CC ones, for the autoisomerization barrier energy we consider the CCSD, CCSDT, CCSDTQ methods and the approximations of CCDT and of CCSDTQ, the CC3 and the CC4 methods, respectively. We can see that the CCSD values are higher than the other CC methods with an energy difference of around 1.05-1.24 \kcalmol between the CCSD and the CCSDT methods. The CCSDT and CCSDTQ autoisomerization barrier energies are closer with 0.25 \kcalmol of energy difference. The energy difference between the CCSDT and its approximation CC3 is about 0.67-0.8 \kcalmol for all the bases whereas the energy difference between the CCSDTQ and its approximate version CC4 is 0.11 \kcalmol.
All the calculations are performed using four basis set, the 6-31+G(d) basis and the aug-cc-pVXZ with X$=$D, T, Q. In the following we will use the notation AVXZ for the aug-cc-pVXZ basis, again with X$=$D, T, Q.
Table \ref{tab:sf_tddft_D2h} shows the results obtained for the vertical transition energies using spin-flip methods and Table \ref{tab:D2h} shows the results obtained with the standard methods. We discuss first the SF-TD-DFT values with hybrid functionals. For the B3LYP functional we can see that the energy differences for each state and throughout all bases are small with the largest one for the $1\,{}^3B_{1g}$ state with 0.012 eV. The same observation can be done for the PBE0 and BH\&HLYP functionals, we can also oberve that adding exact exchange to the functional (20\% of exact exchange for B3LYP, 25\% for PBE0 and 50\% for BH\&HLYP increase the energy difference between the $1\,{}^1B_{1g}$ and the $1\,{}^3B_{1g}$ states. Put another way, increasing the amount of exact exchange in the functional reduces the energy difference between the $1\,{}^3B_{1g}$ and the $\text{X}\,{}^1A_{g}$ states. We can also notice that the vertical energies of the different states do not vary in the same way when adding exact exchange, for instance the energy variation for the $1\,{}^3B_{1g}$ state from PBE0 to BH\&HLYP is around 0.1 eV whereas for the $1\,{}^1B_{1g}$ and the $1\,{}^1A_{1g}$ states this energy variation is about 0.4-0.5 eV and 0.34 eV respectively. For the RSH functionals we can not make the same observation, indeed we can see that for the CAM-B3LYP functional we have that the energy difference between the $1\,{}^1B_{1g}$ and the $1\,{}^3B_{1g}$ states is larger than for the $\omega$B97X-V and LC-$\omega$PBE08 functionals despite the fact that the latter ones have a bigger amount of exact exchange. However we can observe that we have a small energy difference for each state throughout the bases and for each RSH functional. The M06-2X functional is an Hybrid meta-GGA functional and contains 54\% of Hartree-Fock exchange, we can compare the energies with the BH\&HLYP functional and we can see that the energy differences are small with around 0.03-0.08 eV. We can notice that the upper bound of 0.08 eV in the energy differences is due to the $1\,{}^3B_{1g}$ state.
Then we discuss the various ADC scheme (ADC(2)-s, ADC(2)-x and ADC(3)) results. For ADC(2) we have vertical energy differencies of about 0.03 eV for the $1\,{}^3B_{1g}$ state and around 0.06 eV for $1\,{}^1A_{1g}$ state throughout all bases. However for the $1\,{}^1B_{1g}$ state we have an energy difference of about 0.2 eV. For the ADC(2)-x and ADC(3) schemes the calculations with the AVQZ basis are not feasible with our resources. With ADC(2)-x we have similar vertical energies for the triplet and the $1\,{}^1B_{1g}$ states but for the $1\,{}^1A_{1g}$ state the energy difference between the ADC(2) and ADC(2)-x schemes is about 0.4-0.5 eV. The ADC(3) values are closer from the ADC(2) than the ADC(2)-x except for the triplet state for which we have a energy difference of about 0.09-0.14 eV. Now, we look at Table \ref{tab:D2h} to discuss the results of standard methods. First we discuss the CC values, we have computed the vertical energies at the CCSDT and CCSDTQ level as well as their approximation the CC3 and CC4 methods, respectively. We can notice that for the $1\,{}^3B_{1g}$ and the $1\,{}^1B_{1g}$ states the CCSDT and the CC3 values are close with an energy difference of 0.009-0.02 eV for all bases. The energy difference is larger for the $1\,{}^1A_{1g}$ state with around 0.35-0.38 eV. The same observation can be done for CCSDTQ and CC4 with similar vertical energies for all bases. Note that the $1\,{}^3B_{1g}$ state can not be described with these methods.
Then we review the vertical energies obtained with multireference methods. The smallest active space considered is four electrons in four orbitals, for CASSCF(4,4) we have small energy variations throughout bases for the $1\,{}^3B_{1g}$ and the $1\,{}^1A_{1g}$ states but a larger variation for the $1\,{}^1B_{1g}$ state with around 0.1 eV. We can observe that we have the inversion of the states compared to all methods discussed so far between the $1\,{}^1A_{1g}$ and $1\,{}^1B_{1g}$ states with $1\,{}^1B_{1g}$ higher than $1\,{}^1A_{1g}$ due to the lack of dynamical correlation in the CASSCF methods. The $1\,{}^1B_{1g}$ state values in CASSCF(4,4) is much higher than for any of the other methods discussed so far. With CASPT2(4,4) we retrieve the right ordering between the states and we see large energy differences with the CASSCF values. Indeed, we have approximatively 0.22-0.25 eV of energy difference for the triplet state for all bases and 0.32-0.36 eV for the $1\,{}^1A_{1g}$ state, the largest energy difference is for the $1\,{}^1B_{1g}$ state with 1.5-1.6 eV.
For the XMS-CASPT2(4,4) only the $1\,{}^1A_{1g}$ state is described with values similar than for the CAPST2(4,4) method. For the NEVPT2(4,4) methods (SC-NEVPT2 and PC-NEVPT2) the vertical energies are similar for the $1\,{}^1B_{1g}$ and the $1\,{}^1A_{1g}$ states with approximatively 0.002-0.003 eV and 0.02-0.03 eV of energy difference for all bases, respectively. The energy difference for the $1\,{}^1B_{1g}$ state is slightly larger with 0.05 eV for all bases. Note that for this state the vertical energy varies of 0.23 eV from the 6-31+G(d) basis to the AVDZ one. Then we use a larger active space with twelve electrons in twelve orbitals, the CASSF(12,12) values are close to the the CASSCF(4,4) for the triplet state with 0.01-0.02 eV of energy differences. For the $1\,{}^1A_{1g}$ state we have an energy difference of about 0.2 eV between the CASSCF(4,4) and the CASSCF(12,12) values. We can notice that increasing the size of the active space gives the right ordering for the states and we have an energy difference of around 0.7 eV for the $1\,{}^1B_{1g}$ state between CASSCF(4,4) and the CASSCF(12,12) values. The CASPT2(12,12) method decreases the energy of all states compared to the CASSCF(12,12) method, again the decrease is not the same for all states. We have a diminution from CASSCF to CASPT2 of about 0.17 to 0.2 eV for the $1\,{}^3B_{1g}$ and the $1\,{}^1A_{1g}$ states and for the different bases. Again, the energy difference for the $1\,{}^1A_{1g}$ state is larger with 0.5-0.7 eV depending on the basis. In a similar way than with XMS-CASPT2(4,4), the XMS-CASPT(12,12) only describes the $1\,{}^1A_{1g}$ state and the vertical energies for this state are close to the CASPT(12,12) values. For the NEVPT2(12,12) schemes we see that for the triplet and the $1\,{}^1A_{1g}$ states the energies are similar with an energy difference between the SC-NEVPT2(12,12) and the PC-NEVPT2(12,12) values of about 0.03-0.04 eV and 0.02-0.03 eV respectively.
Spin-flip TD-DFT vertical excitation energies (with respect to the singlet $\text{X}\,{}^1A_{g}$ ground state) of the $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{g}$ states of CBD at the $D_{2h}$ rectangular equilibrium geometry of the $\text{X}\,{}^1 A_{g}$ ground state.
% Spin-flip CIS, ADC and CC vertical excitation energies (with respect to the singlet $\text{X}\,{}^1A_{g}$ ground state) of the $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{g}$ states of CBD at the $D_{2h}$ rectangular equilibrium geometry of the $\text{X}\,{}^1 A_{g}$ ground state.
Standard vertical excitation energies (with respect to the singlet $\text{X}\,{}^1A_{g}$ ground state) of the $1\,{}^3B_{1g}$, $1\,{}^1B_{1g}$, and $2\,{}^1A_{g}$ states of CBD at the $D_{2h}$ rectangular equilibrium geometry of the $\text{X}\,{}^1 A_{g}$ ground state.
Table \ref{tab:sf_D4h} shows vertical energies obtained with spin-flip methods and Table \ref{tab:D4h} vertical energies obtained with standard methods. As for the previous geometry we start with the SF-TD-DFT results with hybrid functionals. We can first notice that for the B3LYP and the PBE0 functionals we have a wrong ordering of the first triplet state $1\,{}^3A_{2g}$ and the ground state $\text{X}\,{}^1B_{1g}$. We retrieve the good ordering with the BH\&HLYP functional, so adding exact exchange to the functional allows us to have the right ordering between these two states. For the B3LYP and the PBE0 functionals we have that the energy differences for each states and for all bases are small with 0.004-0.007 eV for the triplet state $1\,{}^3A_{2g}$. We have 0.015-0.021 eV of energy difference for the $2\,{}^1A_{1g}$ state through all bases, we can notice that this state is around 0.13 eV (considering all bases) higher with the PBE0 functional. We can make the same observation for the $1\,{}^1B_{2g}$ state for the B3LYP and the PBE0 functionals, indeed we have small energy differences for all bases and the state is around 0.14-0.15 eV for the PBE0 functional. For the BH\&HLYP functional the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states are higher in energy than for the two other hybrid functionals with about 0.65-0.69 eV higher for the $2\,{}^1A_{1g}$ state and 0.75-0.77 eV for the $1\,{}^1B_{2g}$ state compared to the PBE0 functional. Then, we have the RSH functionals CAM-B3LYP, $\omega$B97X-V and LC-$\omega$PBE08. For these functionals the vertical energies are similar for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states with a maximum energy difference of 0.01-0.02 eV for the $2\,{}^1A_{1g}$ state and 0.005-0.009 eV for the $1\,{}^1B_{2g}$ state considering all bases. The maximum energy difference for the triplet state is larger with 0.047-0.057 eV for all bases. Note that the vertical energies obtained with the RSH functionals are close to the PBE0 ones except that we have the right ordering between the triplet state $1\,{}^3A_{2g}$ and the ground state $\text{X}\,{}^1B_{1g}$. The M06-2X results are closer to the BH\&HLYP ones but the M06-2X vertical energies are always higher than the BH\&HLYP ones. We can notice that the M06-2X energies for the $2\,{}^1A_{1g}$ state are close to the BH\&HLYP energies for the $1\,{}^1B_{2g}$ state. For these two states, when we compared the results obtained with the M06-2X and the BH\&HLYP functionals, we have an energy difference of 0.16-0.17 eV for the $2\,{}^1A_{1g}$ state and 0.17-0.18 eV for the $1\,{}^1B_{2g}$ state considering all bases. For the triplet state $1\,{}^3A_{2g}$ the energy differences are smaller with 0.03-0.04 eV for all bases. The M11 vertical energies are very close to the M06-2X ones for the triplet state with a maximum energy difference of 0.003 eV considering all bases, and are closer to the BH\&HLYP results for the two other states with 0.06-0.07 eV and 0.07-0.08 eV of energy difference for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states, respectively. Then we discuss the various ADC scheme results, note that we were not able to obtain the vertical energies with the AVQZ basis due to computational resources. For the ADC(2)-s scheme we can see that the energy difference for the triplet state are smaller than for the two other states, indeed we have an energy difference of 0.09 eV for the triplet state whereas we have 0.15 eV and 0.25 eV for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states, again when considering all bases. The energy difference for each state and through the bases are similar for the two other ADC schemes. We can notice a large variation of the vertical energies for the $2\,{}^1A_{1g}$ state between ADC(2)-s and ADC(2)-x with around 0.52-0.58 eV through all bases. The ADC(3) vertical energies are very similar to the ADC(2) ones for the $1\,{}^1B_{2g}$ state with an energy difference of 0.01-0.02 eV for all bases, whereas we have an energy difference of 0.04-0.11 eV and 0.17-0.22 eV for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states, respective
Now we look at Table \ref{tab:D4h} for vertical energies obtained with standard methods, we have first the various CC methods, for all the CC methods we were only able to reach the AVTZ basis.
%For the CCSD method we cannot obtain the vertical energy for the $1\,{}^1B_{2g}$ state, we can notice a small energy variation through the bases for the triplet state with around 0.06 eV. The energy variation is larger for the $2\,{}^1A_{1g}$ state with 0.20 eV.
Again, we have computed the vertical energies at the CCSDT and CCSDTQ level as well as their approximation the CC3 and CC4 methods, respectively. For the CC3 method we do not have the vertical energies for the triplet state $1\,{}^3A_{2g}$. Considering all bases for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states we have an energy difference of about 0.15 eV and 0.12 eV, respectively. The CCSDT energies are close to the CC3 ones for the $2\,{}^1A_{1g}$ state with an energy difference of around 0.03-0.06 eV considering all bases. For the $1\,{}^1B_{2g}$ state the energy difference between the CC3 and the CCSDT values is larger with 0.18-0.27 eV. We can make a similar observation between the CC4 and the CCSDTQ values, for the $2\,{}^1A_{1g}$ state we have an energy difference of about 0.01 eV and this time we have smaller energy difference for the $1\,{}^1B_{2g}$ with 0.01 eV. Then we discuss the multireference results and this time we were able to reach the AVQZ basis. Again we consider two different active spaces, four electrons in four orbitals and twelve electrons in twelve orbitals. If we compare the CASSCF(4,4) and CASPT2(4,4) values we can see large energy differences for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states, for the $2\,{}^1A_{1g}$ state we have an energy difference of about 0.67-0.74 eV and 1.65-1.81 eV for the $1\,{}^1B_{2g}$ state. The energy difference is smaller for the triplet state with 0.27-0.31 eV, we can notice that all the CASSCF(4,4) vertical energies are higher than the CASPT2(4,4) ones. Then we have the NEVPT2(4,4) methods (SC-NEVPT2(4,4) and PC-NEVPT2(4,4)) for which the vertical energies are quite similar, however we can notice that we have a different ordering of the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states with $1\,{}^1B_{2g}$ higher in energy than $2\,{}^1A_{1g}$ for the two NEVPT2(4,4) methods. Then we have the results for the same methods but with a larger active space. For the CASSCF(12,12) we have similar values for the triplet state with energy difference of about 0.06 eV for all bases but larger energy difference for the $2\,{}^1A_{1g}$ state with around 0.28-0.29 eV and 0.79-0.81 eV for the $1\,{}^1B_{2g}$ state. Note that from CASSCF(4,4) to CASSCF(12,12) every vertical energy is smaller. We can make the contrary observation for the CASPT2 method where from CASPT2(4,4) to CASPT2(12,12) every vertical energy is larger. For the CASPT2(12,12) method we have similar values for the triplet state and for the $2\,{}^1A_{1g}$ state, considering all bases, with an energy difference of around 0.05-0.06 eV and 0.02-0.05 eV respectively. The energy difference is larger for the $1\,{}^1B_{2g}$ state with about 0.27-0.29 eV. Next we have the NEVPT2 methods, the first observation that we can make is that by increasing the size of the active space we retrieve the right ordering of the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states. Again, every vertical energy is higher in the NEVPT2(12,12) case than for the NEVPT2(4,4) one.
%%% TABLE VI %%%
\begin{squeezetable}
\begin{table}
\caption{
Standard vertical excitation energies (with respect to the singlet $\text{X}\,{}^1B_{1g}$ ground state) of the $1\,{}^3A_{2g}$, $2\,{}^1A_{1g}$, and $1\,{}^1B_{2g}$ states of CBD at the $D_{4h}$ square-planar equilibrium geometry of the $1\,{}^3A_{2g}$ state.
Standard vertical excitation energies (with respect to the singlet $\text{X}\,{}^1B_{1g}$ ground state) of the $1\,{}^3A_{2g}$, $2\,{}^1A_{1g}$, and $1\,{}^1B_{2g}$ states of CBD at the $D_{4h}$ square-planar equilibrium geometry of the $1\,{}^3A_{2g}$ state.
EM, AS, and PFL acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).}
