CBD/Manuscript/CBD.tex
2022-04-06 14:26:14 +02:00

888 lines
68 KiB
TeX

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,longtable,wrapfig,txfonts}
\usepackage[version=4]{mhchem}
%\usepackage{natbib}
%\bibliographystyle{achemso}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{txfonts}
\usepackage{siunitx}
\DeclareSIUnit[number-unit-product = {\,}]
\cal{cal}
\DeclareSIUnit\kcal{\kilo\cal}
\newcommand{\kcalmol}{\si{\kcal\per\mole}}
\usepackage[
colorlinks=true,
citecolor=blue,
breaklinks=true
]{hyperref}
\urlstyle{same}
\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{red}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\trashXB}[1]{\textcolor{darkgreen}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\SupInf}{\textcolor{blue}{supporting information}}
\newcommand{\QP}{\textsc{quantum package}}
\newcommand{\T}[1]{#1^{\intercal}}
%% bold in Table
\newcommand{\bb}[1]{\textbf{#1}}
\newcommand{\rb}[1]{\textbf{\textcolor{red}{#1}}}
\newcommand{\gb}[1]{\textbf{\textcolor{darkgreen}{#1}}}
%geometries
\newcommand{\Dtwo}{$D_{2h}$}
\newcommand{\Dfour}{$D_{4h}$}
\sisetup{range-phrase=--}
\sisetup{range-units=single}
%states
%D2h states
\newcommand{\oneAg}{$1{}^1A_g$}
\newcommand{\tBoneg}{$1{}^3B_{1g}$}
\newcommand{\sBoneg}{$1{}^1B_{1g}$}
\newcommand{\twoAg}{$2{}^1A_g$}
%D4h states
%\newcommand{\oneBoneg}{$1{}^1B_{1g}$} same label as the D2h state
\newcommand{\Atwog}{$1{}^3A_{2g}$}
\newcommand{\Aoneg}{$1{}^1A_{1g}$}
\newcommand{\Btwog}{$1{}^1B_{2g}$}
\begin{document}
% addresses
\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}
\newcommand{\CEISAM}{Universit\'e de Nantes, CNRS, CEISAM UMR 6230, F-44000 Nantes, France}
\title{Reference Energies for Cyclobutadiene: Automerization and Excited States}
\author{Enzo \surname{Monino}}
\email{emonino@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Martial \surname{Boggio-Pasqua}}
\affiliation{\LCPQ}
\author{Anthony \surname{Scemama}}
\affiliation{\LCPQ}
\author{Denis \surname{Jacquemin}}
\affiliation{\CEISAM}
\author{Pierre-Fran\c{c}ois \surname{Loos}}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\begin{abstract}
The cyclobutadiene molecule is a well-known playground for theoretical chemists and is particularly suitable to test ground- and excited-state methods.
Indeed, due to its high spatial symmetry, especially at the $D_{4h}$ square geometry but also in the $D_{2h}$ rectangular arrangement, the ground and excited states of cyclobutadiene exhibit multi-configurational characters and single-reference methods, such as adiabatic time-dependent density-functional theory (TD-DFT) or equation-of-motion coupled cluster (EOM-CC), are notoriously known to struggle in such situations.
In this work, using a large panel of methods and basis sets, we provide an extensive computational study of the automerization barrier (defined as the difference between the square and rectangular ground-state energies) and the vertical excitation energies of the $D_{2h}$ and $D_{4h}$ equilibrium structures.
In particular, selected configuration interaction (SCI), multi-reference perturbation theory (CASSCF, CASPT2, and NEVPT2), and coupled-cluster (CCSD, CC3, CCSDT, CC4, and CCSDTQ) calculations are performed.
The spin-flip formalism, which is known to provide a correct description of states with multi-configurational character, is also tested within TD-DFT (where numerous exchange-correlation functionals are considered) and the algebraic diagrammatic construction [ADC(2)-s, ADC(2)-x, and ADC(3)].
A theoretical best estimate is defined for the automerization barrier and for each vertical transition energy.
\end{abstract}
\maketitle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Despite the fact that excited states are involved in ubiquitous processes such as photochemistry, \cite{Bernardi_1990,Bernardi_1996,Boggio-Pasqua_2007,Klessinger_1995,Olivucci_2010,Robb_2007,VanderLugt_1969} catalysis or in solar cell technology, \cite{Delgado_2010} none of the currently existing methods has shown to provide accurate excitation energies in all scenarios due to the complexity of the process, the size of the systems, environment effects and many other possible factors.
Indeed, each computational model has its own theoretical and/or technical limitations and the number of possible chemical scenarios is so vast that the design of new excited-state methodologies is still a very active field of theoretical quantum chemistry.\cite{Roos_1996,Piecuch_2002b,Dreuw_2005,Krylov_2006,Sneskov_2012,Gonzales_2012,Laurent_2013,Adamo_2013,Dreuw_2015,Ghosh_2018,Blase_2020,Loos_2020a,Hait_2021,Zobel_2021}
Speaking of difficult tasks, the cyclobutadiene (CBD) molecule has been a real challenge for experimental and theoretical chemistry for many decades. \cite{Bally_1980} Due to its antiaromaticity \cite{Minkin_1994} and large angular strain, \cite{Baeyer_1885} CBD presents a high reactivity which made its synthesis a particularly difficult exercise.
In the {\Dfour} symmetry, the simple H\"uckel molecular orbital theory (wrongly) predicts a triplet ground state (Hund's rule) with two singly-occupied frontier orbitals that are degenerate by symmetry, while state-of-the-art \textit{ab initio} methods (correctly) predict an open-shell singlet ground state.
This degeneracy is lifted by the so-called Jahn-Teller effect, \ie, by a descent in symmetry (from {\Dfour} to {\Dtwo} point group) via a geometrical distortion of the molecule, leading to a closed-shell singlet ground state (see below).
This was confirmed by several experimental studies by Pettis and co-workers \cite{Reeves_1969} and others. \cite{Irngartinger_1983,Ermer_1983,Kreile_1986}
In the {\Dtwo} symmetry, the {\oneAg} ground state has a weak multi-configurational character with well-separated frontier orbitals that can be described by single-reference methods.
However, in the {\Dfour} symmetry, the {\sBoneg} ground state has two singly occupied frontier orbitals that are degenerate.
Therefore, one must take into account, at least, two electronic configurations to properly model this multi-configurational scenario.
Of course, single-reference methods are naturally unable to describe such situations.
The singlet ground state, {\sBoneg}, of the square arrangement is a transition state in the automerization reaction between the two rectangular structures (see Fig.~\ref{fig:CBD}), while the lowest triplet state, {\Atwog}, is a minimum on the triplet potential energy surface.
Thus, the automerization barrier (AB) is defined as the difference between the square and rectangular ground-state energies.
The energy of this barrier is estimated, experimentally, in the range of \SIrange{1.6}{10}{\kcalmol}, \cite{Whitman_1982} while previous state-of-the-art \textit{ab initio} calculations yield an energy barrier in the range of \SIrange{7}{9}{\kcalmol}. \cite{Eckert-Maksic_2006,Li_2009,Shen_2012,Zhang_2019}
The lowest-energy excited states of CBD in both symmetries are represented in Fig.~\ref{fig:CBD}, where we have reported the {\oneAg} and {\tBoneg} states for the rectangular geometry and the {\sBoneg} and {\Atwog} states for the square one.
Due to the energy scale, the higher-energy states ({\sBoneg} and {\twoAg} for {\Dtwo} and {\Aoneg} and {\Btwog} for {\Dfour}) are not shown.
Interestingly, the {\twoAg} and {\Aoneg} states have a strong contribution from doubly-excited configurations and these so-called double excitations \cite{Loos_2019} are known to be an absolute nightmare for adiabatic time-dependent density-functional theory (TD-DFT) \cite{Runge_1984,Casida_1995,Tozer_2000,Maitra_2004,Cave_2004,Levine_2006,Elliott_2011,Maitra_2012,Maitra_2017} 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 problem of multi-configurational character and double excitations, we have explored several routes.
The most evident way is to rely on multi-configurational methods, which are naturally designed to address such scenarios.
Among these methods, one can mention the complete-active-space self-consistent field (CASSCF) method, \cite{Roos_1996} its second-order perturbatively-corrected variant (CASPT2) \cite{Andersson_1990,Andersson_1992,Roos_1995a} and the second-order $n$-electron valence state perturbation theory (NEVPT2) formalism. \cite{Angeli_2001,Angeli_2001a,Angeli_2002}
The exponential scaling of the computational cost (with respect to the size of the active space) associated with these methods is the principal limitation to their applicability to large molecules.
Another way to deal with double excitations and multi-reference situations is to use high level truncation of the equation-of-motion (EOM) formalism \cite{Rowe_1968,Stanton_1993} of coupled-cluster (CC) theory. \cite{Kucharski_1991,Kallay_2003,Kallay_2004,Hirata_2000,Hirata_2004}
However, to provide a correct description of these situations, one have to take into account, at the very least, contributions from the triple excitations in the CC expansion. \cite{Watson_2012,Loos_2018a,Loos_2019,Loos_2020b}
Again, due to the scaling of CC methods with the number of basis functions, their applicability is limited to small molecules.
Although multi-reference CC methods have been specifically designed for such purposes, \cite{Jeziorski_1981,Mahapatra_1998,Mahapatra_1999,Lyakh_2012,Kohn_2013} they are computationally demanding and still far from being black-box.
In this context, an interesting alternative to multi-configurational and CC methods is provided by selected configuration interaction (SCI) methods, \cite{Bender_1969,Whitten_1969,Huron_1973,Giner_2013,Evangelista_2014,Giner_2015,Caffarel_2016b,Holmes_2016,Tubman_2016,Liu_2016,Ohtsuka_2017,Zimmerman_2017,Coe_2018,Garniron_2018} which have proven, recently, to be able to provide near full CI (FCI) energies for small molecules for both ground- and excited-state energies. \cite{Caffarel_2014,Caffarel_2016a,Scemama_2016,Holmes_2017,Li_2018,Scemama_2018,Scemama_2018b,Li_2020,Loos_2018a,Chien_2018,Loos_2019,Loos_2020b,Loos_2020c,Loos_2020e,Garniron_2019,Eriksen_2020,Yao_2020,Williams_2020,Veril_2021,Loos_2021,Damour_2021}
For example, the \textit{Configuration Interaction using a Perturbative Selection made Iteratively} (CIPSI) method limits the 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. \cite{Huron_1973,Giner_2013,Giner_2015,Garniron_2017,Garniron_2018,Garniron_2019}
Nonetheless, SCI methods remain very expensive and can be applied to a limited number of situations.
Finally, another option to deal with these chemical scenarios is to rely on the cheaper spin-flip formalism, established by Krylov in 2001, \cite{Krylov_2001a,Krylov_2001b,Krylov_2002,Casanova_2020} where one accesses the ground and doubly-excited states via a single (spin-flip) excitation and deexcitation from the lowest triplet state, respectively.
Obviously, spin-flip methods have their own flaws, especially spin contamination (\ie, the artificial mixing of electronic states with different spin multiplicities) due principally to spin incompleteness in the spin-flip expansion but also to the spin contamination of the reference configuration. \cite{Casanova_2020}
One can address part of this issue by increasing the excitation order or by complementing the spin-incomplete configuration set with the missing configurations, \cite{Sears_2003,Casanova_2008,Huix-Rotllant_2010,Li_2010,Li_2011a,Li_2011b,Zhang_2015,Lee_2018}
both solutions being associated with an additional computational cost.
In the present work, we investigate the accuracy of each family of computational methods mentioned above on the automerization barrier and the low-lying excited states of CBD at the {\Dtwo} and {\Dfour} ground-state geometries.
Computational details are reported in Sec.~\ref{sec:compmet}.
Section \ref{sec:res} is devoted to the discussion of our results.
Finally, our conclusions are drawn in Sec.~\ref{sec:conclusion}.
%%% FIGURE 1 %%%
\begin{figure}
\includegraphics[width=\linewidth]{figure1}
\caption{Pictorial representation of the ground and excited states of CBD and the properties under investigation.
The singlet ground state (S) and triplet (T) properties are represented in black and red, respectively.
The automerization barrier (AB) is also represented.}
\label{fig:CBD}
\end{figure}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compmet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Selected configuration interaction calculations}
\label{sec:SCI}
For the SCI calculations, we rely on the CIPSI algorithm implemented in QUANTUM PACKAGE, \cite{Garniron_2019} which iteratively select determinants in the FCI space.
To treat electronic states on equal footing, we use a state-averaged formalism where the ground and excited states are expanded with the same set of determinants but with different CI coefficients.
Note that the determinant selection for these states are performed simultaneously via the protocol described in Refs.~\onlinecite{Scemama_2019,Garniron_2019}.
For a given size of the variational wave function, the CIPSI energy is the sum of two terms: the variational energy obtained by diagonalization of the CI matrix in the reference space $E_\text{var}$ and a second-order perturbative correction $E_\text{PT2}$ which estimates the contribution of the external determinants that are not included in the variational space at a given iteration.
The sum of these two energies is, for large enough wave functions, an estimate of the FCI energy, \ie, $E_\text{FCI} \approx E_\text{var} + E_\text{PT2}$.
It is possible to estimate more precisely the FCI energy via an extrapolation procedure, where the variational energy is extrapolated to $E_\text{PT2} = 0$. \cite{Holmes_2017}
Excitation energies are then computed as differences of extrapolated total energies. \cite{Chien_2018,Loos_2018a,Loos_2019,Loos_2020b,Loos_2020c}
Additionally, an error bar can be provided via a recent method based on Gaussian random variables as described in Ref.~\onlinecite{Veril_2021}.
This type of extrapolation procedures is now routine in SCI methods as well as other state-of-the-art techniques. \cite{Eriksen_2020,Loos_2020e,Eriksen_2021}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Coupled-cluster calculations}
\label{sec:CC}
Coupled-cluster theory provides a hierarchy of methods that yields increasingly accurate ground state energies by ramping up the maximum excitation degree of the cluster operator: \cite{Cizek_1966,Paldus_1972,Crawford_2000,Piecuch_2002b,Bartlett_2007,Shavitt_2009} CC with singles and doubles (CCSD), \cite{Cizek_1966,Purvis_1982} CC with singles, doubles, and triples (CCSDT), \cite{Noga_1987a,Scuseria_1988} CC with singles, doubles, triples, and quadruples (CCSDTQ), \cite{Oliphant_1991,Kucharski_1991a,Kucharski_1992} etc.
As mentioned above, CC theory can be extended to excited states via the EOM formalism, \cite{Rowe_1968,Stanton_1993} where one diagonalizes the similarity-transformed Hamiltonian in a CI basis of excited determinants yielding the following systematically improvable family of methods for neutral excited states:\cite{Noga_1987a,Koch_1990b,Kucharski_1991,Stanton_1993,Christiansen_1998,Kucharski_2001,Kowalski_2001,Kallay_2003,Kallay_2004,Hirata_2000,Hirata_2004} EOM-CCSD, EOM-CCSDT, EOM-CCSDTQ, etc.
In the following, we will omit the prefix EOM for the sake of conciseness.
Alternatively to the ``complete'' CC models, one can also employed the CC2, \cite{Christiansen_1995,Hattig_2000} CC3 \cite{Christiansen_1995,Koch_1995} and CC4 \cite{Kallay_2005,Matthews_2020,Loos_2021} methods which 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.
Here, we have performed CC calculations using various codes.
Typically, CCSD, CCSDT, and CCSDTQ as well as CC3 and CC4 calculations are achieved with CFOUR, \cite{Matthews_2020} with which only singlet excited states can be computed.
In some cases, we have also computed (singlet and triplet) excitation energies and properties (such as the percentage of single excitations involved in a given transition, namely $\%T_1$) at the CC3 level with DALTON \cite{Aidas_2014} and CCSDT level with MRCC. \cite{mrcc}
To avoid having to perform multi-reference CC calculations and because one cannot perform high-level CC calculations in the restricted open-shell or unrestricted formalisms, it is worth mentioning that, for the {\Dfour} arrangement, we have considered the lowest \textit{closed-shell} singlet state {\Aoneg} as reference.
Hence, the open-shell ground state {\sBoneg} and the {\Btwog} state appear as a deexcitation and an excitation, respectively.
With respect to {\Aoneg}, {\sBoneg} has a dominant double excitation character, while {\Btwog} have a dominant single excitation character, hence their contrasting convergence behaviors with respect to the order of the CC expansion (see below).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Multi-configurational calculations}
\label{sec:Multi}
State-averaged CASSCF (SA-CASSCF) calculations are performed with MOLPRO. \cite{Werner_2020}
For each excited state, a set of state-averaged orbitals is computed by taking into account the excited state of interest as well as the ground state (even if it has a different symmetry).
Two active spaces have been considered: (i) a minimal (4e,4o) active space including valence $\pi$ orbitals, and (ii) an extended (12e,12o) active space where we have additionally included the $\sigma_\text{CC}$ and $\sigma_\text{CC}^*$ orbitals.
For ionic states, like the {\sBoneg} state of CBD, it is particularly important to take into account the $\sigma$-$\pi$ coupling. \cite{Davidson_1996,Angeli_2009,BenAmor_2020}
On top of this CASSCF treatment, CASPT2 calculations are performed within the RS2 contraction scheme, while the NEVPT2 energies are computed within both the partially contracted (PC) and strongly contracted (SC) schemes. \cite{Angeli_2001,Angeli_2001a,Angeli_2002}
Note that PC-NEVPT2 is theoretically more accurate than SC-NEVPT2 due to the larger number of external configurations and greater flexibility.
In order to avoid the intruders state problem in CASPT2, a real-valued level shift of \SI{0.3}{\hartree} is set, \cite{Roos_1995b,Roos_1996} with an additional ionization-potential-electron-affinity (IPEA) shift of \SI{0.3}{\hartree} to avoid systematic overestimation of the vertical excitation energies. \cite{Ghigo_2004,Schapiro_2013,Zobel_2017,Sarkar_2022}
%For the sake of comparison, in some cases, we have also performed multi-reference CI (MRCI) calculations.
All these calculations are also carried out with MOLPRO. \cite{Werner_2020}
%and extended multistate (XMS) CASPT2 calculations are also performed in the cas of strong mixing between states with same spin and spatial symmetries.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Spin-flip calculations}
\label{sec:sf}
Within the spin-flip formalism, one considers the lowest triplet state as reference instead of the singlet ground state.
Ground-state energies are then computed as sums of the triplet reference state energy and the corresponding deexcitation energy.
Likewise, excitation energies with respect to the singlet ground state are computed as differences of excitation energies with respect to the reference triplet state.
Nowadays, spin-flip techniques are broadly accessible thanks to intensive developments in the electronic structure community (see Ref.~\onlinecite{Casanova_2020} and references therein).
Here, we explore the spin-flip version \cite{Lefrancois_2015} of the algebraic-diagrammatic construction \cite{Schirmer_1982} (ADC) using the standard and extended second-order ADC schemes, SF-ADC(2)-s \cite{Trofimov_1997,Dreuw_2015} and SF-ADC(2)-x, \cite{Dreuw_2015} as well as its third-order version, SF-ADC(3). \cite{Dreuw_2015,Trofimov_2002,Harbach_2014}
These calculations are performed using Q-CHEM 5.2.1. \cite{qchem}
The spin-flip version of our recently proposed composite approach, namely SF-ADC(2.5), where one simply averages the SF-ADC(2)-s and SF-ADC(3) energies, is also tested in the following. \cite{Loos_2020d}
We have also carried out spin-flip calculations within the TD-DFT framework (SF-TD-DFT), \cite{Shao_2003} and these are also performed with Q-CHEM 5.2.1. \cite{qchem}
The B3LYP, \cite{Becke_1988b,Lee_1988a,Becke_1993b} PBE0 \cite{Adamo_1999a,Ernzerhof_1999} and BH\&HLYP hybrid functionals are considered, which contain 20\%, 25\%, 50\% of exact exchange, respectively.
These calculations are labeled as SF-TD-BLYP, SF-TD-B3LYP, SF-TD-PBE0, and SF-TD-BH\&HLYP in the following.
Additionally, we have also computed SF-TD-DFT excitation energies using range-separated hybrid (RSH) functionals: CAM-B3LYP (19\% of short-range exact exchange and 65\% at long range), \cite{Yanai_2004a} LC-$\omega$PBE08 (0\% of short-range exact exchange and 100\% at long range), \cite{Weintraub_2009a} and $\omega$B97X-V (16.7\% of short-range exact exchange and 100\% at long range). \cite{Mardirossian_2014}
Finally, the hybrid meta-GGA functional M06-2X (54\% of exact exchange) \cite{Zhao_2008} and the RSH meta-GGA functional M11 (42.8\% of short-range exact exchange and 100\% at long range) \cite{Peverati_2011} are also employed.
Note that all SF-TD-DFT calculations are done within the Tamm-Dancoff approximation. \cite{Hirata_1999}
%EOM-SF-CCSD and EOM-SF-CC(2,3) are also performed with Q-CHEM 5.2.1.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Theoretical best estimates}
\label{sec:TBE_compdet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
When technically possible, each level of theory is tested with four Gaussian basis sets, namely, 6-31+G(d) and aug-cc-pVXZ with X $=$ D, T, and Q. \cite{Dunning_1989}
This helps us to assess the convergence of each property with respect to the size of the basis set.
More importantly, for each studied quantity (i.e., the automerization barrier and the vertical excitation energies), we provide a theoretical best estimate (TBE) established in the aug-cc-pVTZ basis.
In most cases, these TBEs are defined using extrapolated CCSDTQ/aug-cc-pVTZ values or, in a single occasion, NEVPT2(12,12).
The extrapolation of the CCSDTQ/aug-cc-pVTZ values is done via a ``pyramidal'' scheme, where we employ systematically the most accurate level of theory and the largest basis set available.
For example, when CC4/aug-cc-pVTZ and CCSDTQ/aug-cc-pVDZ data are available, we proceed via the following basis set extrapolation:
\begin{equation}
\label{eq:aug-cc-pVTZ}
\Delta \Tilde{E}^{\text{CCSDTQ}}_{\text{aug-cc-pVTZ}}
= \Delta E^{\text{CCSDTQ}}_{\text{aug-cc-pVDZ}}
+ \qty[ \Delta E^{\text{CC4}}_{\text{aug-cc-pVTZ}} - \Delta E^{\text{CC4}}_{\text{aug-cc-pVDZ}} ],
\end{equation}
while, when only CCSDTQ/6-31G+(d) values are available, we further extrapolate the CCSDTQ/aug-cc-pVDZ value as follows:
\begin{equation}
\label{eq:aug-cc-pVDZ}
\Delta \Tilde{E}^{\text{CCSDTQ}}_{\text{aug-cc-pVDZ}}
= \Delta E^{\text{CCSDTQ}}_{6-31\text{G}+\text{(d)}}
+ \qty[ \Delta E^{\text{CC4}}_{\text{aug-cc-pVDZ}} - \Delta E^{\text{CC4}}_{6-31\text{G}+\text{(d)}} ].
\end{equation}
If we lack the CC4 data, we can follow the same philosophy and rely on CCSDT.
For example,
\begin{equation}
\label{eq:CC4_aug-cc-pVTZ}
\Delta \Tilde{E}^{\text{CC4}}_{\text{aug-cc-pVTZ}}
= \Delta E^{\text{CC4}}_{\text{aug-cc-pVDZ}}
+ \qty[ \Delta E^{\text{CCSDT}}_{\text{aug-cc-pVTZ}} - \Delta E^{\text{CCSDT}}_{\text{aug-cc-pVDZ}} ],
\end{equation}
and so on.
If neither CC4, nor CCSDT are feasible, then we rely on NEVPT2(12,12).
The procedure for each extrapolated value is explicitly mentioned as a footnote.
Note that, due to error bar inherently linked to the CIPSI calculations (see Sec.~\ref{sec:SCI}), these are mostly used as an additional safety net to further check the convergence of the CCSDTQ estimates.
Additional tables gathering these TBEs as well as literature data for the automerization barrier and the vertical excitation energies can be found in the {\SupInf}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%================================================
\subsection{Geometries}
\label{sec:geometries}
Two different sets of geometries obtained with different levels of theory are considered for the automerization barrier and the excited states of the CBD molecule.
First, because the automerization barrier is obtained as a difference of energies computed at distinct geometries, it is paramount to obtain these at the same level of theory.
However, due to the fact that the ground state of the square arrangement is a transition state of singlet open-shell nature, it is technically difficult to optimize the geometry with high-order CC methods.
Therefore, we rely on CASPT2(12,12)/aug-cc-pVTZ for both the {\Dtwo} and {\Dfour} ground-state structures.
(Note that these optimizations are done without IPEA shift.)
Second, because the vertical transition energies are computed for a particular equilibrium geometry, we can afford to use different methods for the rectangular and square structures.
Hence, we rely on CC3/aug-cc-pVTZ to compute the equilibrium geometry of the {\oneAg} state in the rectangular ({\Dtwo}) arrangement and the restricted open-shell (RO) version of CCSD(T)/aug-cc-pVTZ to obtain the equilibrium geometry of the {\Atwog} state in the square ({\Dfour}) arrangement.
These two geometries are the lowest-energy equilibrium structure of their respective spin manifold (see Fig.~\ref{fig:CBD}).
The cartesian coordinates of these geometries are provided in the {\SupInf}.
Table \ref{tab:geometries} reports the key geometrical parameters obtained at these levels of theory as well as previous geometries computed by Manohar and Krylov at the CCSD(T)/cc-pVTZ level.
%================================================
%%% TABLE I %%%
\begin{squeezetable}
\begin{table}
\caption{Optimized geometries associated with several states of CBD computed with various levels of theory.
Bond lengths are in \si{\angstrom} and angles ($\angle$) are in degree.}
\label{tab:geometries}
\begin{ruledtabular}
\begin{tabular}{lllrrr}
State & Method & \ce{C=C} & \ce{C-C} & \ce{C-H} & $\angle\,\ce{H-C=C}$ \\
\hline
{\Dtwo} ({\oneAg}) &
CASPT2(12,12)/aug-cc-pVTZ \fnm[1] & 1.354 & 1.566 & 1.077 & 134.99 \\
&CC3/aug-cc-pVTZ \fnm[1] & 1.344 & 1.565 & 1.076 & 135.08 \\
&CCSD(T)/cc-pVTZ \fnm[2] & 1.343 & 1.566 & 1.074 & 135.09\\
{\Dfour} ({\sBoneg}) &
CASPT2(12,12)/aug-cc-pVTZ \fnm[1] & 1.449 & 1.449 & 1.076 & 135.00 \\
{\Dfour} ({\Atwog}) &
CASPT2(12,12)/aug-cc-pVTZ \fnm[1] & 1.445 & 1.445 & 1.076 & 135.00 \\
&RO-CCSD(T)/aug-cc-pVTZ \fnm[1] & 1.439 & 1.439 & 1.075 & 135.00\\
&RO-CCSD(T)/cc-pVTZ \fnm[2] & 1.439 & 1.439 & 1.073 & 135.00\\
\end{tabular}
\end{ruledtabular}
\fnt[1]{This work.}
\fnt[2]{From Ref.~\onlinecite{Manohar_2008}.}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%================================================
\subsection{Automerization barrier}
\label{sec:auto}
%%% TABLE I %%%
\begin{squeezetable}
\begin{table}
\caption{Automerization barrier (in \kcalmol) of CBD computed with various computational methods and basis sets.
The values in square parenthesis have been obtained by extrapolation via the procedure described in the corresponding footnote.
The TBE/aug-cc-pVTZ value is highlighted in bold.}
\label{tab:auto_standard}
\begin{ruledtabular}
\begin{tabular}{lrrrr}
& \mc{4}{c}{Basis sets} \\
\cline{2-5}
Method & 6-31+G(d) & aug-cc-pVDZ & aug-cc-pVTZ & aug-cc-pVQZ\\
\hline
%SF-CIS & $2.64$ & $2.82$ & $3.43$ & $3.43$ \\
%SF-TD-BLYP & $23.57$ & $23.62$ & $24.23$ & $24.22$ \\
SF-TD-B3LYP & $18.59$ & $18.64$ & $19.34$ & $19.34$ \\
SF-TD-PBE0& $17.18$ & $17.19$ & $17.88$ & $17.88$ \\
SF-TD-BH\&HLYP & $11.90$ & $12.02$ & $12.72$ & $12.73$ \\
SF-TD-M06-2X & $9.32$ & $9.62$ & $10.35$ & $10.37$ \\
SF-TD-CAM-B3LYP& $18.05$ & $18.10$ & $18.83$ & $18.83$ \\
SF-TD-$\omega$B97X-V & $18.26$ & $18.24$ & $18.94$ & $18.92$ \\
SF-TD-LC-$\omega$PBE08 & $19.05$ & $18.98$ & $19.74$ & $19.71$ \\[0.1cm]
SF-TD-M11 & $11.03$ & $10.25$ & $11.22$ & $11.12$ \\
SF-ADC(2)-s & $6.69$ & $6.98$ & $8.63$ & \\
SF-ADC(2)-x & $8.63$ & $8.96$ &$10.37$ & \\
SF-ADC(2.5) & $7.36$ & $7.76$ & $9.11$ & \\
SF-ADC(3) & $8.03$ & $8.54$ & $9.58$ \\[0.1cm]
CASSCF(4,4) & $6.17$ & $6.59$ & $7.38$ & $7.41$ \\
CASPT2(4,4) & $6.56$ & $6.87$ & $7.77$ & $7.93$ \\
SC-NEVPT2(4,4) & $7.95$ & $8.31$ & $9.23$ & $9.42$ \\
PC-NEVPT2(4,4) & $7.95$ & $8.33$ & $9.24$ & $9.41$ \\
CASSCF(12,12) & $10.19$ & $10.75$ & $11.59$ & $11.62$ \\
CASPT2(12,12) & $7.24$ & $7.53$ & $8.51$ & $8.71$ \\
SC-NEVPT2(12,12) & $7.10$ & $7.32$ & $8.29$ & $8.51$ \\
PC-NEVPT2(12,12) & $7.12$ & $7.33$ & $8.28$ & $8.49$ \\[0.1cm]
CCSD & $8.31$ & $8.80$ & $9.88$ & $10.10$ \\
CC3 & $6.59$ & $6.89$ & $7.88$ & $8.06$ \\
CCSDT & $7.26$ & $7.64$ & $8.68$ &$[ 8.86]$\fnm[1] \\
CC4 & $7.40$ & $7.78$ & $[ 8.82]$\fnm[2] & $[ 9.00]$\fnm[3]\\
CCSDTQ & $7.51$ & $[ 7.89]$\fnm[4]& $[\bf 8.93]$\fnm[5]& $[ 9.11]$\fnm[6]\\
%\alert{CIPSI} & $7.91\pm 0.21$ & $8.58\pm 0.14$ & & \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Value obtained using CCSDT/aug-cc-pVTZ corrected by the difference between CC3/aug-cc-pVQZ and CC3/aug-cc-pVTZ.}
\fnt[2]{Value obtained using CC4/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
\fnt[3]{Value obtained using CC4/aug-cc-pVTZ corrected by the difference between CCSDT/aug-cc-pVQZ and CCSDT/aug-cc-pVTZ.}
\fnt[4]{Value obtained using CCSDTQ/6-31+G(d) corrected by the difference between CC4/aug-cc-pVDZ basis and CC4/6-31+G(d).}
\fnt[5]{TBE value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
\fnt[6]{Value obtained using CCSDTQ/aug-cc-pVTZ corrected by the difference between CC4/aug-cc-pVQZ and CC4/aug-cc-pVTZ.}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%%% FIGURE II %%%
\begin{figure*}
\includegraphics[width=\linewidth]{AB_AVTZ}
\caption{Automerization barrier (in \si{\kcalmol}) of CBD at various levels of theory using the aug-cc-pVTZ basis.
See {\SupInf} for the raw data.}
\label{fig:AB}
\end{figure*}
%%% %%% %%% %%%
The results concerning the automerization barrier are reported in Table \ref{tab:auto_standard} for various basis sets and shown in Fig.~\ref{fig:AB} for the aug-cc-pVTZ basis.
First, one can see large variations of the energy barrier at the SF-TD-DFT level, with differences as large as \SI{10}{\kcalmol} between the different functionals for a given basis set.
Nonetheless, it is clear that the performance of a given functional is directly linked to the amount of exact exchange at short range.
Indeed, hybrid functionals with a large fraction of short-range exact exchange (\eg, BH\&HLYP, M06-2X, and M11) perform significantly better than the functionals having a small fraction of short-range exact exchange (\eg, B3LYP, PBE0, CAM-B3LYP, $\omega$B97X-V, and LC-$\omega$PBE08).
However, they are still off by \SIrange{1}{4}{\kcalmol} from the TBE reference value.
For the RSH functionals, the automerization barrier is much less sensitive to the amount of longe-range exact exchange.
Another important feature of SF-TD-DFT is the fast convergence of the energy barrier with the size of the basis set. \cite{Loos_2019d}
With the augmented double-$\zeta$ basis, the SF-TD-DFT results are basically converged to sub-\kcalmol accuracy, which is a drastic improvement compared to wave function approaches where this type of convergence is reached with the augmented triple-$\zeta$ basis.
For the SF-ADC family of methods, the energy differences are much smaller with a maximum deviation of \SI{2}{\kcalmol} between different versions.
In particular, we observe that SF-ADC(2)-s and SF-ADC(3), which scale as $\order*{N^5}$ and $\order*{N^6}$ respectively (where $N$ is the number of basis functions), under- and overestimate the automerization barrier, making SF-ADC(2.5) a good compromise with an error of only \SI{0.18}{\kcalmol} with the aug-cc-pVTZ basis.
Nonetheless, at a $\order*{N^5}$ computational scaling, SF-ADC(2)-s is particularly accurate, even compared to high-order CC methods (see below).
We note that SF-ADC(2)-x is probably not worth its extra cost [as compared to SF-ADC(2)-s] as it overestimates the energy barrier even more than SF-ADC(3).
Overall, even with the best exchange-correlation functional, SF-TD-DFT is clearly outperformed by the more expensive SF-ADC methods.
Concerning the multi-reference approaches with the minimal (4e,4o) active space, the TBEs are bracketed by the CASPT2 and NEVPT2 values that differ by approximately \SI{1.5}{\kcalmol} for all bases.
In this case, the NEVPT2 values are fairly accurate with differences below half a \si{\kcalmol} compared to the TBEs.
The CASSCF results predict an even lower barrier than CASPT2 due to the well known lack of dynamical correlation at the CASSCF level.
For the larger (12e,12o) active space, we see larger differences of the order of \SI{3}{\kcalmol} through all the bases between CASSCF and the second-order variants (CASPT2 and NEVPT2).
Thought, the deviations between CASPT2(12,12) and NEVPT2(12,12) are much smaller with an energy difference of around \SIrange{0.1}{0.2}{\kcalmol} for all bases, CASPT2 being slightly more accurate than NEVPT2 in this case.
For each basis set, both CASPT2(12,12) and NEVPT2(12,12) are less than a \si{\kcalmol} away from the TBEs.
For the two active spaces that we have considered here, the PC- and SC-NEVPT2 schemes provide nearly identical barriers independently of the size of the one-electron basis.
Finally, for the CC family of methods, we observe the usual systematic improvement following the series CCSD $<$ CC3 $<$ CCSDT $<$ CC4 $<$ CCSDTQ, which is also linked to their increase in computational cost: $\order*{N^6}$, $\order*{N^7}$, $\order*{N^8}$, $\order*{N^9}$, and $\order*{N^{10}}$, respectively.
Note that the introduction of the triple excitations is clearly mandatory to have an accuracy beyond SF-TD-DFT, while it is also clear that the iterative triples and quadruples can be included approximately via the CC3 and CC4 methods.
%================================================
%================================================
\subsection{Vertical excitation energies}
\label{sec:states}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{{\Dtwo} rectangular geometry}
\label{sec:D2h}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% TABLE II %%%
\begin{squeezetable}
\begin{table}
\caption{
Spin-flip TD-DFT and ADC vertical excitation energies (with respect to the singlet {\oneAg} ground state) of the {\tBoneg}, {\sBoneg}, and {\twoAg} states of CBD at the {\Dtwo} rectangular equilibrium geometry of the {\oneAg} ground state.}
\label{tab:sf_D2h}
\begin{ruledtabular}
\begin{tabular}{llrrr}
& \mc{4}{r}{Excitation energies (eV)} \hspace{0.5cm}\\
\cline{3-5}
Method & Basis & {\tBoneg} & {\sBoneg} & {\twoAg} \\
\hline
% SF-TD-BLYP & 6-31+G(d) & $1.829$ & $1.926$ & $3.755$ \\
% & aug-cc-pVDZ & $1.828$ & $1.927$ & $3.586$ \\
% & aug-cc-pVTZ & $1.825$ & $1.927$ & $3.546$ \\
% & aug-cc-pVQZ & $1.825$ & $1.927$ & $3.528$ \\[0.1cm]
SF-TD-B3LYP & 6-31+G(d) & $1.706$ & $2.211$ & $3.993$ \\
& aug-cc-pVDZ & $1.706$ & $2.204$ & $3.992$ \\
& aug-cc-pVTZ & $1.703$ & $2.199$ & $3.988$ \\
& aug-cc-pVQZ & $1.703$ & $2.199$ & $3.989$\\[0.1cm]
SF-TD-PBE0 & 6-31+G(d) & $1.687$ & $2.314$ & $4.089$ \\
& aug-cc-pVDZ & $1.684$ & $2.301$ & $4.085$ \\
& aug-cc-pVTZ & $1.682$ & $2.296$ & $4.081$ \\
& aug-cc-pVQZ & $1.682$ & $2.296$ & $4.079$\\[0.1cm]
SF-TD-BH\&HLYP & 6-31+G(d) & $1.552$ & $2.779$ & $4.428$ \\
& aug-cc-pVDZ & $1.546$ & $2.744$ & $4.422$ \\
& aug-cc-pVTZ & $1.540$ & $2.732$ & $4.492$ \\
& aug-cc-pVQZ & $1.540$ & $2.732$ & $4.415$ \\[0.1cm]
SF-TD-M06-2X & 6-31+G(d) & $1.477$ & $2.835$ & $4.378$ \\
& aug-cc-pVDZ & $1.467$ & $2.785$ & $4.360$ \\
& aug-cc-pVTZ & $1.462$ & $2.771$ & $4.357$ \\
& aug-cc-pVQZ & $1.458$ & $2.771$ & $4.352$ \\[0.1cm]
SF-TD-CAM-B3LYP & 6-31+G(d) & $1.750$ & $2.337$ & $4.140$ \\
& aug-cc-pVDZ & $1.745$ & $2.323$ & $4.140$ \\
& aug-cc-pVTZ & $1.742$ & $2.318$ & $4.138$ \\
& aug-cc-pVQZ & $1.743$ & $2.319$ & $4.138$ \\[0.1cm]
SF-TD-$\omega$B97X-V & 6-31+G(d) & $1.810$ & $2.377$ & $4.220$ \\
& aug-cc-pVDZ & $1.800$ & $2.356$ & $4.217$ \\
& aug-cc-pVTZ & $1.797$ & $2.351$ & $4.213$ \\
& aug-cc-pVQZ & $1.797$ & $2.351$ & $4.213$ \\[0.1cm]
SF-TD-LC-$\omega $PBE08 & 6-31+G(d) & $1.917$ & $2.445$ & $4.353$ \\
& aug-cc-pVDZ & $1.897$ & $2.415$ & $4.346$ \\
& aug-cc-pVTZ & $1.897$ & $2.415$ & $4.348$ \\
& aug-cc-pVQZ & $1.897$ & $2.415$ & $4.348$ \\[0.1cm]
SF-TD-M11 & 6-31+G(d) & $1.566$ & $2.687$ & $4.292$ \\
& aug-cc-pVDZ & $1.546$ & $2.640$ & $4.267$ \\
& aug-cc-pVTZ & $1.559$ & $2.651$ & $4.300$ \\
& aug-cc-pVQZ & $1.557$ & $2.650$ & $4.299$ \\[0.1cm]
%SF-CIS & 6-31+G(d) & $1.514$ & $3.854$ & $5.379$ \\
%& aug-cc-pVDZ & $1.487$ & $3.721$ & $5.348$ \\
%& aug-cc-pVTZ & $1.472$ & $3.701$ & $5.342$ \\
%& aug-cc-pVQZ & $1.471$ & $3.702$ & $5.342$ \\[0.1cm]
SF-ADC(2)-s & 6-31+G(d) & $1.577$ & $3.303$ & $4.196$ \\
& aug-cc-pVDZ & $1.513$ & $3.116$ & $4.114$ \\
& aug-cc-pVTZ & $1.531$ & $3.099$ & $4.131$ \\
& aug-cc-pVQZ & $1.544$ & $3.101$ & $4.140$ \\[0.1cm]
SF-ADC(2)-x & 6-31+G(d) & $1.557$ & $3.232$ & $3.728$ \\
& aug-cc-pVDZ & $1.524$ & $3.039$ & $3.681$ \\
& aug-cc-pVTZ & $1.539$ & $3.031$ & $3.703$ \\[0.1cm]
SF-ADC(2.5) & 6-31+G(d) & $1.496$ & $3.328$ & $4.219$ \\
& aug-cc-pVDZ & $1.468$ & $3.148$ & $4.161$ \\
& aug-cc-pVTZ & $1.475$ & $3.131$ & $4.178$ \\[0.1cm]
SF-ADC(3) & 6-31+G(d) & $1.435$ & $3.352$ & $4.242$ \\
& aug-cc-pVDZ & $1.422$ & $3.180$ & $4.208$ \\
& aug-cc-pVTZ & $1.419$ & $3.162$ & $4.224$ \\
% SF-EOM-CCSD & 6-31+G(d) & $1.663$ & $3.515$ & $4.275$ \\
% & aug-cc-pVDZ & $1.611$ & $3.315$ & $3.856$ \\
% & aug-cc-pVTZ & $1.609$ & $3.293$ & $4.245$ \\[0.1cm]
%SF-EOM-CC(2,3) & 6-31+G(d) & $1.490$ & $3.333$ & $4.061$ \\
%& aug-cc-pVDZ & $1.464$ & $3.156$ & $4.027$ \\
\end{tabular}
\end{ruledtabular}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%%% TABLE IV %%%
\begin{squeezetable}
\begin{table}
\caption{
Vertical excitation energies (with respect to the {\oneAg} ground state) of the {\tBoneg}, {\sBoneg}, and {\twoAg} states of CBD at the {\Dtwo} rectangular equilibrium geometry of the {\oneAg} ground state.
The values in square parenthesis have been obtained by extrapolation via the procedure described in the corresponding footnote.
The TBE/aug-cc-pVTZ values are highlighted in bold.}
\label{tab:D2h}
\begin{ruledtabular}
\begin{tabular}{llrrr}
& & \mc{3}{c}{Excitation energies (eV)} \\
\cline{3-5}
Method & Basis & {\tBoneg} & {\sBoneg} & {\twoAg} \\
\hline
CASSCF(4,4) &6-31+G(d)& $1.662$ & $4.657$ & $4.439$ \\
& aug-cc-pVDZ & $1.672$ & $4.563$ & $4.448$ \\
& aug-cc-pVTZ & $1.670$ & $4.546$ & $4.441$ \\
& aug-cc-pVQZ & $1.671$ & $4.549$ & $4.440$ \\[0.1cm]
CASPT2(4,4) &6-31+G(d)& $1.440$ & $3.162$ & $4.115$ \\
& aug-cc-pVDZ & $1.414$ & $2.971$ & $4.068$ \\
& aug-cc-pVTZ & $1.412$ & $2.923$ & $4.072$ \\
& aug-cc-pVQZ & $1.417$ & $2.911$ & $4.081$ \\[0.1cm]
%XMS-CASPT2(4,4) &6-31+G(d)& & & $4.151$ \\
%& aug-cc-pVDZ & & & $4.105$ \\
%& aug-cc-pVTZ & & & $4.114$ \\
%& aug-cc-pVQZ && & $4.125$ \\[0.1cm]
SC-NEVPT2(4,4) &6-31+G(d)& $1.407$ & $2.707$ & $4.145$ \\
& aug-cc-pVDZ & $1.381$ & $2.479$ & $4.109$ \\
& aug-cc-pVTZ & $1.379$ & $2.422$ & $4.108$ \\
& aug-cc-pVQZ & $1.384$ & $2.408$ & $4.125$ \\[0.1cm]
PC-NEVPT2(4,4) &6-31+G(d)& $1.409$ & $2.652$ & $4.120$ \\
& aug-cc-pVDZ & $1.384$ & $2.424$ & $4.084$ \\
& aug-cc-pVTZ & $1.382$ & $2.368$ & $4.083$ \\
& aug-cc-pVQZ & $1.387$ & $2.353$ & $4.091$ \\[0.1cm]
%MRCI(4,4) &6-31+G(d)& $1.564$ & $3.802$ & $4.265$ \\
%& aug-cc-pVDZ & $1.558$ & $3.670$ & $4.254$ \\
%& aug-cc-pVTZ & $1.568$ & $3.678$ & $4.270$ \\
%& aug-cc-pVQZ & $1.574$ & $3.681$ & $4.280$ \\[0.1cm]
CASSCF(12,12) &6-31+G(d)& $1.675$ & $3.924$ & $4.220$ \\
& aug-cc-pVDZ & $1.685$ & $3.856$ & $4.221$ \\
& aug-cc-pVTZ & $1.686$ & $3.844$ & $4.217$ \\
& aug-cc-pVQZ & $1.687$ & $3.846$ & $4.216$ \\[0.1cm]
CASPT2(12,12) &6-31+G(d)& $1.508$ & $3.407$ & $4.099$ \\
& aug-cc-pVDZ & $1.489$ & $3.256$ & $4.044$ \\
& aug-cc-pVTZ & $1.480$ & $3.183$ & $4.043$ \\
& aug-cc-pVQZ & $1.482$ & $3.163$ & $4.047$ \\[0.1cm]
%XMS-CASPT2(12,12) &6-31+G(d)& && $4.111$ \\
%& aug-cc-pVDZ & & & $4.056$ \\
%& aug-cc-pVTZ & & & $4.059$ \\
%& aug-cc-pVQZ & & & $4.065$ \\[0.1cm]
SC-NEVPT2(12,12) &6-31+G(d)& $1.522$ & $3.409$ & $4.130$ \\
& aug-cc-pVDZ & $1.511$ & $3.266$ & $4.093$ \\
& aug-cc-pVTZ & $1.501$ & $3.188$ & $4.086$ \\
& aug-cc-pVQZ & $1.503$ & $3.167$ & $4.088$ \\[0.1cm]
PC-NEVPT2(12,12) &6-31+G(d)& $1.487$ & $3.296$ & $4.103$ \\
& aug-cc-pVDZ & $1.472$ & $3.141$ & $4.064$ \\
& aug-cc-pVTZ & $\bf 1.462$ & $3.063$ & $4.056$ \\
& aug-cc-pVQZ & $1.464$ & $3.043$ & $4.059$ \\[0.1cm]
%MRCI(12,12) &6-31+G(d)& & & $4.125$ \\[0.1cm]
CC3 &6-31+G(d)& $1.420$ & $3.341$ & $4.658$ \\
& aug-cc-pVDZ & $1.396$ & $3.158$ & $4.711$ \\
& aug-cc-pVTZ & $1.402$ & $3.119$ & $4.777$ \\
& aug-cc-pVQZ & $1.409$ & $3.113$ & $4.774$ \\[0.1cm]
CCSDT &6-31+G(d)& $1.442$ & $3.357$ & $4.311$ \\
& aug-cc-pVDZ & $1.411$ & $3.175$ & $4.327$ \\
& aug-cc-pVTZ & $1.411$ & $3.139$ & $4.429$ \\[0.1cm]
CC4 &6-31+G(d)& & $3.343$ & $4.067$ \\
& aug-cc-pVDZ & & $3.164$ & $4.041$ \\
& aug-cc-pVTZ & & $[3.128]$\fnm[1] & $[4.143]$\fnm[1]\\[0.1cm]
CCSDTQ &6-31+G(d)& & $3.340$ & $4.073$ \\
& aug-cc-pVDZ & & $[3.161]$\fnm[2]& $[4.047]$\fnm[2] \\
& aug-cc-pVTZ & & $[\bf 3.125]$\fnm[3]& $[\bf 4.149]$\fnm[3]\\[0.1cm]
CIPSI &6-31+G(d)& $1.486\pm 0.005$ & $3.348\pm 0.024$ & $4.084\pm 0.012$ \\
& aug-cc-pVDZ & $1.458\pm 0.009$ & $3.187\pm 0.035$ & $4.04\pm 0.04$ \\
& aug-cc-pVTZ & $1.461\pm 0.030$ & $3.142\pm 0.035$ & $4.03\pm 0.09$ \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Value obtained using CC4/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
\fnt[2]{Value obtained using CCSDTQ/6-31+G(d) corrected by the difference between CC4/aug-cc-pVDZ and CC4/6-31+G(d).}
\fnt[3]{TBE value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%%% FIGURE III %%%
\begin{figure*}
\includegraphics[width=\linewidth]{D2h}
\caption{Vertical excitation energies of the {\tBoneg}, {\sBoneg}, and {\twoAg} states at the {\Dtwo} rectangular equilibrium geometry of the {\oneAg} ground state using the aug-cc-pVTZ basis.
See {\SupInf} for the raw data.}
%Purple, orange, green, blue and black lines correspond to the SF-TD-DFT, SF-ADC, multi-reference, CC, and TBE values, respectively.}
\label{fig:D2h}
\end{figure*}
%%% %%% %%% %%%
Table \ref{tab:sf_D2h} reports, at the {\Dtwo} rectangular equilibrium geometry of the {\oneAg} ground state, the vertical transition energies associated with the {\tBoneg}, {\sBoneg}, and {\twoAg} states obtained using the spin-flip formalism, while Table \ref{tab:D2h} gathers the same quantities obtained with the multi-reference, CC, and CIPSI methods.
Considering the aug-cc-pVTZ basis, the evolution of the vertical excitation energies with respect to the level of theory is illustrated in Fig.~\ref{fig:D2h}.
At the CC3/aug-cc-pVTZ level, the percentage of single excitation involved in the {\tBoneg}, {\sBoneg}, and {\twoAg} are 99\%, 95\%, and 1\%, respectively.
Therefore, the two formers are dominated by single excitations, while the latter state is a genuine double excitation.
First, let us discuss basis set effects at the SF-TD-DFT level (Table \ref{tab:sf_D2h}).
As expected, these are found to be small and the results are basically converged to the basis set limit with the triple-$\zeta$ basis, which is definitely not the case for the wave function methods. \cite{Giner_2019}
Regarding now the accuracy of the vertical excitation energies, again, we clearly see that, for each transition, the functionals with the largest amount of short-range exact exchange (\eg, BH\&HLYP, M06-2X, and M11) are the most accurate.
Functionals with large fraction of exact exchange are known to perform best in the SF-TD-DFT framework as the Hartree-Fock exchange term is the only non-vanishing term in the spin-flip block. \cite{Shao_2003}
However, their overall accuracy remains average especially for the singlet states, {\sBoneg} and {\twoAg}, with error of the order of \SIrange{0.2}{0.5}{\eV} compared to the TBEs.
The triplet state, {\tBoneg}, is much better described with errors below \SI{0.1}{\eV}.
Note that, as evidenced by the data reported in {\SupInf}, none of these states exhibit a strong spin contamination.
%First, we discuss 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 {\tBoneg} state with \SI{0.012}{\eV}.
%The same observation can be done for the PBE0 and BH\&HLYP functionals, we can also observe 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 {\sBoneg} and the {\tBoneg} states.
%Put another way, increasing the amount of exact exchange in the functional reduces the energy difference between the {\tBoneg} and the {\oneAg} 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 {\tBoneg} state from PBE0 to BH\&HLYP is around \SI{0.1}{\eV} whereas for the {\sBoneg} and the {\twoAg} states this energy variation is about \SIrange{0.4}{0.5}{\eV} and \SI{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 {\sBoneg} and the {\tBoneg} states is larger than for the $\omega$B97X-V and LC-$\omega$PBE08 functionals despite the fact that the latter ones have a larger 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 \SIrange{0.03}{0.08}{\eV}.
%We can notice that the upper bound of \SI{0.08}{\eV} in the energy differences is due to the {\tBoneg} state.
%The M11 vertical energies are close to the BH\&HLYP ones for the triplet and the first singlet excited states with \SIrange{0.01}{0.02}{\eV} and \SIrange{0.08}{0.10}{\eV} of energy difference, respectively.
%For the {\twoAg} state the M11 energies are closer to the $\omega$B97X-V ones with \SIrange{0.05}{0.09}{\eV} of energy difference.
Second, we discuss the various SF-ADC schemes (Table \ref{tab:sf_D2h}), \ie, SF-ADC(2)-s, SF-ADC(2)-x, and SF-ADC(3).
At the SF-ADC(2)-s level, going from the smallest 6-31+G(d) basis to the largest aug-cc-pVQZ basis, we see a small decrease in vertical excitation energies of about \SI{0.03}{\eV} for the {\tBoneg} state and around \SI{0.06}{\eV} for {\twoAg} state, while the transition energy of the {\sBoneg} state drops more significantly by about \SI{0.2}{\eV}.
[The SF-ADC(2)-x and SF-ADC(3) calculations with aug-cc-pVQZ were not feasible with our computational resources.]
These basis set effects are fairly transferable to the other wave function methods that we have considered here.
This further motivates the ``pyramidal'' extrapolation scheme that we have employed to produce the TBE values (see Sec.~\ref{sec:TBE_compdet}).
Again, the extended version, SF-ADC(2)-x, does not seem to be relevant in the present context with much larger errors than the other schemes.
Also, as reported previously, \cite{Loos_2020d} SF-ADC(2)-s and SF-ADC(3) have mirror error patterns which makes SF-ADC(2.5) particularly accurate with a maximum error of \SI{0.029}{\eV} for the doubly-excited state {\twoAg}.
Let us now move to the discussion of the results concerning standard wave function methods that are reported in Table \ref{tab:D2h}.
Regarding the multi-configurational calculations, the most striking result is the poor description of the {\sBoneg} ionic state, especially with the (4e,4o) active space where CASSCF predicts this state higher in energy than the {\twoAg} state.
Of course, the PT2 correction is able to correct the root-flipping problem but cannot provide quantitative excitation energies due to the poor zeroth-order treatment.
Another ripple effect of the unreliability of the reference wave function is the large difference between CASPT2 and NEVPT2 that differ by half an \si{\eV}.
This feature is characteristic of the inadequacy of the active space.
For the two other states, {\tBoneg} and {\twoAg}, the errors at the CASPT2(4,4) and NEVPT2(4,4) levels are much smaller and typically below \SI{0.1}{\eV}.
Using a larger active resolves most of these issues: CASSCF predicts the correct state ordering (though the ionic state is still badly described in term of energetics), CASPT2 and NEVPT2 excitation energies are much closer, and their accuracy is often improved (especially for the triplet state) although it is difficult to reach chemical accuracy (\ie, an error below \SI{0.043}{\eV}) on a systematic basis.
Finally, for the CC models (Table \ref{tab:D2h}), the two states with a large $\%T_1$ value, {\tBoneg} and {\sBoneg}, are already extremely accurate at the CC3 level, and systematically improved by CCSDT and CC4.
For the doubly-excited state, {\twoAg}, the convergence of the CC expansion is much slower but it is worth pointing out that the inclusion of approximate quadruples via CC4 is particularly effective in the present case.
The CCSDTQ excitation energies (which are used to define the TBEs) are systematically within the error bar of the CIPSI calculations, which confirms the outstanding performance of CC methods including quadruple excitations in the context of excited states.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{{\Dfour} square-planar geometry}
\label{sec:D4h}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% TABLE VI %%%
\begin{squeezetable}
\begin{table}
\caption{
Spin-flip TD-DFT and ADC vertical excitation energies (with respect to the singlet {\sBoneg} ground state) of the {\Atwog}, {\Aoneg}, and {\Btwog} states of CBD at the {\Dfour} square-planar equilibrium geometry of the {\Atwog} state.}
\label{tab:sf_D4h}
\begin{ruledtabular}
\begin{tabular}{llrrr}
& \mc{4}{r}{Excitation energies (eV)} \hspace{0.5cm}\\
\cline{3-5}
Method & Basis & {\Atwog} & {\Aoneg} & {\Btwog} \\
\hline
%SF-CIS & 6-31+G(d) & $0.355$ & $2.742$ & $3.101$ \\
%& aug-cc-pVDZ & $0.318$ & $2.593$ & $3.052$ \\
%& aug-cc-pVTZ & $0.305$ & $2.576$ & $3.053$ \\
%& aug-cc-pVQZ & $0.306$ & $2.577$ & $3.056$ \\[0.1cm]
SF-TD-B3LYP & 6-31+G(d) & $-0.016$ & $0.487$ & $0.542$ \\
& aug-cc-pVDZ & $-0.019$ & $0.477$ & $0.536$ \\
& aug-cc-pVTZ & $-0.020$ & $0.472$ & $0.533$ \\
& aug-cc-pVQZ & $-0.020$ & $0.473$ & $0.533$ \\[0.1cm]
SF-TD-PBE0 & 6-31+G(d) & $-0.012$ & $0.618$ & $0.689$ \\
& aug-cc-pVDZ & $-0.016$ & $0.602$ & $0.680$ \\
& aug-cc-pVTZ & $-0.019$ & $0.597$ & $0.677$ \\
& aug-cc-pVQZ & $-0.018$ & $0.597$ & $0.677$ \\[0.1cm]
SF-TD-BH\&HLYP & 6-31+G(d) & $0.064$ & $1.305$ & $1.458$ \\
& aug-cc-pVDZ & $0.051$ & $1.260$ & $1.437$ \\
& aug-cc-pVTZ & $0.045$ & $1.249$ & $1.431$ \\
& aug-cc-pVQZ & $0.046$ & $1.250$ & $1.432$ \\[0.1cm]
SF-TD-M06-2X & 6-31+G(d) & $0.102$ & $1.476$ & $1.640$ \\
& aug-cc-pVDZ & $0.086$ & $1.419$ & $1.611$ \\
& aug-cc-pVTZ & $0.078$ & $1.403$ & $1.602$ \\
& aug-cc-pVQZ & $0.079$ & $1.408$ & $1.607$ \\[0.1cm]
SF-TD-CAM-B3LYP & 6-31+G(d) & $0.021$ & $0.603$ & $0.672$ \\
& aug-cc-pVDZ & $0.012$ & $0.585$ & $0.666$ \\
& aug-cc-pVTZ & $0.010$ & $0.580$ & $0.664$ \\
& aug-cc-pVQZ & $0.010$ & $0.580$ & $0.664$ \\[0.1cm]
SF-TD-$\omega $B97X-V & 6-31+G(d) & $0.040$ & $0.600$ & $0.670$ \\
& aug-cc-pVDZ & $0.029$ & $0.576$ & $0.664$ \\
& aug-cc-pVTZ & $0.026$ & $0.572$ & $0.662$ \\
& aug-cc-pVQZ & $0.026$ & $0.572$ & $0.662$ \\[0.1cm]
SF-TD-LC-$\omega $PBE08 & 6-31+G(d) & $0.078$ & $0.593$ & $0.663$ \\
& aug-cc-pVDZ & $0.060$ & $0.563$ & $0.659$ \\
& aug-cc-pVTZ & $0.058$ & $0.561$ & $0.658$ \\
& aug-cc-pVQZ & $0.058$ & $0.561$ & $0.659$ \\[0.1cm]
SF-TD-M11 & 6-31+G(d) & $0.102$ & $1.236$ & $1.374$ \\
& aug-cc-pVDZ & $0.087$ & $1.196$ & $1.362$ \\
& aug-cc-pVTZ & $0.081$ & $1.188$ & $1.359$ \\
& aug-cc-pVQZ & $0.080$ & $1.185$ & $1.357$ \\[0.1cm]
SF-ADC(2)-s & 6-31+G(d) & $0.345$ & $1.760$ & $2.096$ \\
& aug-cc-pVDZ & $0.269$ & $1.656$ & $1.894$ \\
& aug-cc-pVTZ & $0.256$ & $1.612$ & $1.844$ \\[0.1cm]
SF-ADC(2)-x & 6-31+G(d) & $0.264$ & $1.181$ & $1.972$ \\
& aug-cc-pVDZ & $0.216$ & $1.107$ & $1.760$ \\
& aug-cc-pVTZ & $0.212$ & $1.091$ & $1.731$ \\[0.1cm]
SF-ADC(2.5) & 6-31+G(d) & $0.234$ & $1.705$ & $2.087$ \\
& aug-cc-pVDZ & $0.179$ & $1.614$ & $1.886$ \\
& aug-cc-pVTZ & $0.168$ & $1.594$ & $1.849$ \\[0.1cm]
SF-ADC(3) & 6-31+G(d) & $0.123$ & $1.650$ & $2.078$ \\
& aug-cc-pVDZ & $0.088$ & $1.571$ & $1.878$ \\
& aug-cc-pVTZ & $0.079$ & $1.575$ & $1.853$ \\
\end{tabular}
\end{ruledtabular}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%%% TABLE V %%%
\begin{squeezetable}
\begin{table}
\caption{
Vertical excitation energies (with respect to the {\sBoneg} ground state) of the {\Atwog}, {\Aoneg}, and {\Btwog} states of CBD at the {\Dfour} square-planar equilibrium geometry of the {\Atwog} state.
The values in square brackets have been obtained by extrapolation via the procedure described in the corresponding footnote.
The TBE/aug-cc-pVTZ values are highlighted in bold.}
\label{tab:D4h}
\begin{ruledtabular}
\begin{tabular}{llrrr}
& \mc{3}{r}{Excitation energies (eV)} \hspace{0.1cm}\\
\cline{3-5}
Method & Basis & {\Atwog} & {\Aoneg} & {\Btwog} \\
\hline
CASSCF(4,4) & 6-31+G(d) & $0.447$ & $2.257$ & $3.549$ \\
& aug-cc-pVDZ & $0.438$ & $2.240$ & $3.443$ \\
& aug-cc-pVTZ & $0.434$ & $2.234$ & $3.424$ \\
& aug-cc-pVQZ & $0.435$ & $2.235$ & $3.427$ \\[0.1cm]
CASPT2(4,4) & 6-31+G(d) & $0.176$ & $1.588$ & $1.899$ \\
& aug-cc-pVDZ & $0.137$ & $1.540$ & $1.708$ \\
& aug-cc-pVTZ & $0.128$ & $1.506$ & $1.635$ \\
& aug-cc-pVQZ & $0.128$ & $1.498$ & $1.612$ \\[0.1cm]
SC-NEVPT2(4,4) & 6-31+G(d) & $0.083$ & $1.520$ & $1.380$ \\
& aug-cc-pVDZ & $0.037$ & $1.465$ & $1.140$ \\
& aug-cc-pVTZ & $0.024$ & $1.428$ & $1.055$ \\
& aug-cc-pVQZ & $0.024$ & $1.420$ & $1.030$ \\[0.1cm]
PC-NEVPT2(4,4) & 6-31+G(d) & $0.085$ & $1.496$ & $1.329$ \\
& aug-cc-pVDZ & $0.039$ & $1.440$ & $1.088$ \\
& aug-cc-pVTZ & $0.026$ & $1.403$ & $1.003$ \\
& aug-cc-pVQZ & $0.026$ & $1.395$ & $0.977$ \\[0.1cm]
%MRCI(4,4) & 6-31+G(d) & $0.297$ & $1.861$ & $2.571$ \\
% & aug-cc-pVDZ & $0.273$ & $1.823$ & $2.419$ \\
% & aug-cc-pVTZ & $0.271$ & $1.824$ & $2.415$ \\
% & aug-cc-pVQZ & $0.273$ & $1.825$ & $2.413$ \\[0.1cm]
CASSCF(12,12) & 6-31+G(d) & $0.386$ & $1.974$ & $2.736$ \\
& aug-cc-pVDZ & $0.374$ & $1.947$ & $2.649$ \\
& aug-cc-pVTZ & $0.370$ & $1.943$ & $2.634$ \\
& aug-cc-pVQZ & $0.371$ & $1.945$ & $2.637$ \\[0.1cm]
CASPT2(12,12) & 6-31+G(d) & $0.235$ & $1.635$ & $2.170$ \\
& aug-cc-pVDZ & $0.203$ & $1.588$ & $2.015$ \\
& aug-cc-pVTZ & $0.183$ & $1.538$ & $1.926$ \\
& aug-cc-pVQZ & $0.179$ & $1.522$ & $1.898$ \\[0.1cm]
SC-NEVPT2(12,12) & 6-31+G(d) & $0.218$ & $1.644$ & $2.143$ \\
& aug-cc-pVDZ & $0.189$ & $1.600$ & $1.991$ \\
& aug-cc-pVTZ & $0.165$ & $1.546$ & $1.892$ \\
& aug-cc-pVQZ & $0.160$ & $1.529$ & $1.862$ \\[0.1cm]
PC-NEVPT2(12,12) & 6-31+G(d) & $0.189$ & $1.579$ & $2.020$ \\
& aug-cc-pVDZ & $0.156$ & $1.530$ & $1.854$ \\
& aug-cc-pVTZ & $0.131$ & $1.476$ & $1.756$ \\
& aug-cc-pVQZ & $0.126$ & $1.460$ & $1.727$ \\[0.1cm]
CCSD & 6-31+G(d) & $0.148$ & $1.788$ & \\
& aug-cc-pVDZ & $0.100$ & $1.650$ & \\
& aug-cc-pVTZ & $0.085$ & $1.600$ & \\
& aug-cc-pVQZ & $0.084$ & $1.588$ & \\[0.1cm]
CC3 & 6-31+G(d) & & $1.809$ & $2.836$ \\
& aug-cc-pVDZ & & $1.695$ & $2.646$ \\
& aug-cc-pVTZ & & $1.662$ & $2.720$ \\[0.1cm]
CCSDT & 6-31+G(d) & $0.210$ & $1.751$ & $2.565$ \\
& aug-cc-pVDZ & $0.165$ & $1.659$ & $2.450$ \\
& aug-cc-pVTZ & $0.149$ & $1.631$ & $2.537$ \\[0.1cm]
CC4 & 6-31+G(d) & & $1.604$ & $2.121$ \\
& aug-cc-pVDZ & & $1.539$ & $1.934$ \\
& aug-cc-pVTZ & & $[1.511 ]$\fnm[1] &$[2.021 ]$\fnm[1] \\[0.1cm]
CCSDTQ & 6-31+G(d) & $0.205$ & $1.593$ & $2.134$ \\
& aug-cc-pVDZ & $[0.160]$\fnm[2] & $[1.528 ]$\fnm[4]&$[1.947]$\fnm[4] \\
& aug-cc-pVTZ & $[\bf 0.144]$\fnm[3] & $[\bf 1.500 ]$\fnm[5]&$[\bf 2.034]$\fnm[5] \\ [0.1cm]
CIPSI & 6-31+G(d) & $0.201\pm 0.003$ & $1.602\pm 0.007$ & $2.13\pm 0.04$ \\
& aug-cc-pVDZ & $0.157\pm 0.003$ & $1.587\pm 0.005$ & $2.102\pm 0.027$ \\
& aug-cc-pVTZ & $0.169\pm 0.029$ & $1.63\pm 0.05$ & \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Value obtained using CC4/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
\fnt[2]{Value obtained using CCSDTQ/6-31+G(d) corrected by the difference between CCSDT/aug-cc-pVDZ and CCSDT/6-31+G(d).}
\fnt[3]{TBE value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
\fnt[4]{Value obtained using CCSDTQ/6-31+G(d) corrected by the difference between CC4/aug-cc-pVDZ and CC4/6-31+G(d).}
\fnt[5]{TBE value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
\end{table}
\end{squeezetable}
%%% %%% %%% %%%
%%% FIGURE IV %%%
\begin{figure*}
\includegraphics[width=\linewidth]{D4h}
\caption{Vertical excitation energies (in \si{\eV}) of the {\Atwog}, {\Aoneg}, and {\Btwog} states at the {\Dfour} square-planar equilibrium geometry of the {\Atwog} state using the aug-cc-pVTZ basis.
See {\SupInf} for the raw data.}
%Purple, orange, green, blue and black lines correspond to the SF-TD-DFT, SF-ADC, multi-reference, CC, and TBE values, respectively.}
\label{fig:D4h}
\end{figure*}
%%% %%% %%% %%%
In Table \ref{tab:sf_D4h}, we report, at the {\Dfour} square planar equilibrium geometry of the {\Atwog} state, the vertical transition energies associated with the {\Atwog}, {\Aoneg}, and {\Btwog} states obtained using the spin-flip formalism, while Table \ref{tab:D4h} gathers the same quantities obtained with the multi-reference, CC, and CIPSI methods.
The vertical excitation energies computed at various levels of theory are depicted in Fig.~\ref{fig:D4h} for the aug-cc-pVTZ basis.
Unfortunately, due to technical limitations, we could not compute $\%T_1$ values associated with the {\Atwog}, {\Aoneg}, and {\Btwog} excited states in the {\Dfour} symmetry.
However, it is clear from the inspection of the wave function that, with respect to the {\sBoneg} ground state, {\Atwog} and {\Btwog} are dominated by single excitations, while {\Aoneg} is strongly dominated by double excitations.
As for the previous geometry we start by discussing the SF-TD-DFT results (Table \ref{tab:sf_D4h}), and in particular the singlet-triplet gap, \ie, the energy difference between {\sBoneg} and {\Atwog}.
For all functionals, this gap is small (basically below \SI{0.1}{\eV} while the TBE value is \SI{0.144}{\eV}) but it is worth mentioning that B3LYP and PBE0 predict a negative singlet-triplet gap (hence a triplet ground state).
Increasing the fraction of exact exchange in hybrids or relying on RSHs (even with a small amount of short-range exact exchange) allows to recover a positive gap and a singlet ground state.
At the SF-TD-DFT level, the energy gap between the two singlet excited states, {\Aoneg} and {\Btwog}, is particularly small and grows moderately with the amount of short-range exact exchange.
The influence of the exact exchange on the singlet energies is quite significant with an energy difference of the order of \SI{1}{\eV} between the functional with the smallest amount of exact exchange (B3LYP) and the functional with the largest amount (M06-2X).
As for the excitation energies computed at the {\Dtwo} ground-state equilibrium structure and the automerization barrier, functionals with a large fraction of short-range exact exchange deliver much more accurate results.
Yet, the transition energy to {\Btwog} is off by more than half an \si{\eV} compared to the TBE, while the doubly-excited state is much closer to the reference value (errors of \SI{-0.251}{}, \SI{-0.097}{}, and \SI{-0.312}{\eV} for BH\&HLYP, M06-2X, and M11, respectively)
Again, for all the excited states, the basis set effects are extremely small at the SF-TD-DFT level.
Note again that the $\expval*{S^2}$ values reported in {\SupInf} show that there is no significant spin contamination in these excited states.
Next, we discuss the various ADC schemes (Table \ref{tab:sf_D4h}) where we were not able to compute the vertical energies with the aug-cc-pVQZ basis due to our limited computational resources.
Overall, we observe similar trends than the ones mentioned in Sec.~\ref{sec:D2h}.
Concerning the singlet-triplet gap, each scheme predicts it to be positive.
Although it provides a decent singlet-triplet gap value, SF-ADC(2)-x seems to particularly struggle with the singlet excited states ({\Aoneg} and {\Btwog}), especially for the doubly-excited state {\Aoneg} where it underestimates the vertical excitation energy by \SI{0.4}{\eV}.
Again, averaging the SF-ADC(2)-s and SF-ADC(3) transition energies is beneficial in most cases at the exception of {\Aoneg}.
Although the basis set effects are larger than at the SF-TD-DFT level, they remain quite moderate at the SF-ADC level, and for any wave function method in general.
Then, we discuss the multi-reference results (Table \ref{tab:D4h}).
For both active spaces, expectedly, CASSCF does not provide a quantitive energetic description of the excited states, although it is worth mentioning that the right state ordering is preserved.
This is, of course, magnified with the (4e,4o) active space for which the second-order perturbative treatment is unable to provide a satisfying description due to the restricted active space.
In particular SC-NEVPT2(4,4)/aug-cc-pVTZ and PC-NEVPT2(4,4)/aug-cc-pVTZ underestimate the singlet-triplet gap by \SI{0.072}{} and \SI{0.097}{\eV} and, more importantly, flip the ordering of {\Aoneg} and {\Btwog}.
Although {\Aoneg} is not badly described, the excitation energy of the ionic state {\Btwog} is off by \SI{1}{\eV}.
Thanks to the IPEA shift in CASPT2(4,4), the singlet-triplet gap is accurate and the state ordering remains correct but the ionic state is still far from being well described.
The (12e,12o) active space significantly alleviate these effects, and, as usually, the agreement between CASPT2 and NEVPT2 is very much improved for each state, though the accuracy of multi-configurational approaches remains questionable for the ionic state with an error up to \SI{-0.278}{\eV} at the PC-NEVPT2(12,12)/aug-cc-pVTZ level.
Finally, let us consider the excitation energies computed with various CC models and gathered in Table \ref{tab:D4h}.
As mentioned in Sec.~\ref{sec:CC}, we remind the reader that these calculations are performed by considering the {\Aoneg} state as reference, and that, therefore,
{\sBoneg} and {\Btwog} are obtained as a deexcitation and an excitation, respectively.
Consequently, with respect to {\Aoneg}, {\sBoneg} has a dominant double excitation character, while {\Btwog} have a dominant single excitation character.
This explains why one observes a slower convergence of the transition energies in the case of {\sBoneg} as shown in Fig.~\ref{fig:D4h}.
It is clear from the results of Table \ref{tab:D4h} that, if one wants to reach high accuracy, it is mandatory to include quadruple excitations.
Indeed, at the CCSDT/aug-cc-pVTZ level, the singlet-triplet gap is already very accurate (off by \SI{0.005}{\eV} only) while the excitation energies of the singlet states are still \SI{0.131}{} and \SI{0.503}{\eV} away from their respective TBE.
These deviations drop to \SI{0.011}{} and \SI{-0.013}{\eV} at the CC4/aug-cc-pVTZ level.
As a final comment, we can note that the agreement between our CCSDTQ-based TBEs and the CIPSI calculations are consistent if one takes into account the extrapolation error (see Sec.~\ref{sec:SCI}).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion}
\label{sec:conclusion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the present study, we have benchmarked a larger number of computational methods on the automerization barrier and the vertical energies of the cyclobutadiene (CBD) molecule in its square ({\Dfour}) and rectangular ({\Dtwo}) arrangements, for which we have defined theoretical best estimates (TBEs) based on extrapolated CCSDTQ/aug-cc-pVTZ data.
The main take-home messages of the present work can be summarized as follows:
\begin{itemize}
\item Within the SF-TD-DFT framework, we advice to use exchange-correlation (hybrids or range-separated hybrids) with a large fraction of short-range exact exchange.
This has been shown to be beneficial for the automerization barrier and the vertical excitation energies in the {\Dtwo} and {\Dfour} arrangements.
\item At the SF-ADC level, we have found that the extended scheme, SF-ADC(2)-x, systematically worsen the results compared to the cheaper standard version, SF-ADC(2)-s.
Moreover, as previously reported, SF-ADC(2)-s and SF-ADC(3) have opposite error patterns which means that SF-ADC(2.5) is an excellent compromise.
\item For the {\Dfour} square planar structure, a faithful energetic description of the excited states is harder to reach at the SF-TD-DFT level because of the strong multi-configurational character.
In such scenario, the SF-TD-DFT excitation energies can exhibit errors of the order of \SI{1}{\eV} compared to the TBEs.
However, it was satisfying to see that the spin-flip version of ADC can lower these errors to \SIrange{0.1}{0.2}{\eV}.
\item Concerning the multi-configurational methods, we have found that NEVPT2 and CASPT2 can provide different excitation energies for the small (4e,4o) active space, but they becomes very similar when the larger (12e,12o) active space is considered.
From a more general perspective, a significant difference between NEVPT2 and CASPT2 can be seen as a warning that the active space has been poorly chosen.
\item The ionic states remain a struggle for both CASPT2 and NEVPT2, even with the (12e,12o) active space.
\item In the context of CC methods, although the inclusion of triple excitations is very satisfactory in most cases, the inclusion of quadruples excitation (via CC4 or CCSDTQ) is mandatory to reach high accuracy (especially in the case of doubly-excited states).
We also point out that, within the error bar linked to the CIPSI extrapolation procedure, CCSDTQ and CIPSI yield similar excitation energies, hence confirming the outstanding accuracy of CCSDTQ in the context of molecular excited states.
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%%%
\acknowledgements{
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).
This work was performed using HPC ressources from CALMIP (Toulouse) under allocation 2022-18005.}
%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supporting Information Available}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Included in the {\SupInf} are the raw data, additional calculations and geometries, and the cartesian coordinates of the various optimized geometries.
%================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{CBD}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}