CBD/Manuscript/CBD.tex
2022-04-04 09:56:20 +02:00

1099 lines
91 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: Autoisomerization 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 autoisomerization 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 autoisomerization 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_2002,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 autoisomerization 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 multi-reference calculations yield an energy barrier in the range of \SIrange{6}{7}{\kcalmol}. \cite{Eckert-Maksic_2006}
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 autoisomerization barrier and the low-lying excited states of CBD at the {\Dtwo} and {\Dfour} ground-state geometries.
Computational details are reported in Section \ref{sec:compmet}.
Section \ref{sec:res} is devoted to the discussion of our results.
Finally, our conclusions are drawn in Section \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 autoisomerization 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 (16.7\% 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}
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 autoisomerisation 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 Subsection \ref{sec:SCI}), these are mostly used as an additional safety net to further check the convergence of the CCSDTQ estimates.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 autoisomerization barrier and the excited states of the CBD molecule.
First, because the autoisomerization 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 with the 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 all 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{Autoisomerization barrier}
\label{sec:auto}
%%% TABLE I %%%
\begin{squeezetable}
\begin{table}
\caption{Autoisomerization 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.}
\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-M11 & $11.03$ & $10.25$ & $11.22$ & $11.12$ \\
SF-TD-LC-$\omega$PBE08 & $19.05$ & $18.98$ & $19.74$ & $19.71$ \\[0.1cm]
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]& $[ 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]{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{Autoisomerization barrier (in \si{\kcalmol}) of CBD at various levels of theory using the aug-cc-pVTZ basis.}
\label{fig:AB}
\end{figure*}
%%% %%% %%% %%%
The results concerning the autoisomerization 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} for a given basis set between the different functionals.
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 autoisomerization 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.0}{\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 autoisomerization 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 smaller (4e,4o) active space, the TBEs are bracketed by the CASPT2 and NEVPT2 values that differ by approximately \SI{1.5}{\kcalmol} for all the bases.
In this case, the NEVPT2 values are fairly accurate with differences below half a \si{\kcalmol} with 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).
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 with the increase in computational cost of $\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{Excited States}
\label{sec:states}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{{\Dtwo} symmetry}
\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_tddft_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.}
\label{tab:D2h}
\begin{ruledtabular}
\begin{tabular}{llrrr}
& & \mc{3}{c}{Excitation energies (eV)} \\
\cline{3-5}
Method & Basis & {\tBoneg} & {\sBoneg} & {\twoAg} \\
\hline
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 & & $\left[3.128\right]$\fnm[1] & $\left[4.143\right]$\fnm[1]\\[0.1cm]
CCSDTQ &6-31+G(d)& & $3.340$ & $4.073$ \\
& aug-cc-pVDZ & & $\left[3.161\right]$\fnm[2]& $\left[4.047\right]$\fnm[2] \\
& aug-cc-pVTZ & & $\left[3.125\right]$\fnm[3]& $\left[4.149\right]$\fnm[3]\\[0.1cm]
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 & $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]
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]{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.}
%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_tddft_D2h} reports, at the {\Dtwo} rectangular equilibrium geometry of the {\oneAg} ground state, the singlet and triplet 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 and CC 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.
\alert{First, let us discuss the basis set effects at the SF-TD-DFT level.
These are small and the results are basically converged to the basis set limit with the triple-$\zeta$ basis.
Again, we clearly see that the functional with the largest amount of short-range exact exchange are the most accurate ones.
However, their accuracy remains average especially for the {\sBoneg} state (due to spin contamination?).
Mention that the amount of exact exchange is key for SF-TD-DFT as this is the only non-vanishing term in the spin-flip block. \cite{Shao_2003}}
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.
Then we discuss the various ADC scheme (ADC(2)-s, ADC(2)-x and ADC(3)) results. For ADC(2) we have vertical energy differences of about 0.03 eV for the {\tBoneg} state and around 0.06 eV for {\twoAg} state throughout all bases.
However for the {\sBoneg} state we have an energy difference of about 0.2 eV.
For the ADC(2)-x and ADC(3) schemes the calculations with the aug-cc-pVQZ basis are not feasible with our resources.
With ADC(2)-x we have similar vertical energies for the triplet and the {\sBoneg} states but for the {\twoAg} state the energy difference between the ADC(2) and ADC(2)-x schemes is about \SIrange{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 \SIrange{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 {\tBoneg} and the {\sBoneg} states the CCSDT and the CC3 values are close with an energy difference of \SIrange{0.01}{0.02}{\eV} for all bases.
The energy difference is larger for the {\twoAg} state with around \SIrange{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 {\tBoneg} state can not be described with these methods.
\alert{Due to the lack of dynamic correlation, CASSCF poorly describes the {\sBoneg} ionic state.
Although the autoisomerization barrier is not bad with the (4e,4o) active space, the {\sBoneg} is poorly describe mainly due to the poor CASSCF treatement that the second-order correction cannot fix.}
Then we review the vertical energies obtained with multi-reference 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 {\tBoneg} and the {\twoAg} states but a larger variation for the {\sBoneg} state with around \SI{0.1}{\eV}.
We can observe that we have the inversion of the states compared to all methods discussed so far between the {\twoAg} and {\sBoneg} states with {\sBoneg} higher than {\twoAg} due to the lack of dynamical correlation in the CASSCF methods.
The {\sBoneg} 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 \SIrange{0.22}{0.25}{\eV} of energy difference for the triplet state for all bases and \SIrange{0.32}{0.36}{\eV} for the {\twoAg} state, the largest energy difference is for the {\sBoneg} state with \SIrange{1.5}{1.6}{\eV}.
For the XMS-CASPT2(4,4) only the {\twoAg} 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 {\sBoneg} and the {\twoAg} states with approximatively \SIrange{0.002}{0.003}{\eV} and \SIrange{0.02}{0.03}{\eV} of energy difference for all bases, respectively.
The energy difference for the {\sBoneg} state is slightly larger with \SI{0.05}{\eV} for all bases.
Note that for this state the vertical energy varies of \SI{0.23}{eV} from the 6-31+G(d) basis to the aug-cc-pVDZ one.
Then we use a larger active space with twelve electrons in twelve orbitals, the CASSF(12,12) values are close to the CASSCF(4,4) value for the triplet state with 0.01-0.02 eV of energy differences.
For the {\twoAg} state we have an energy difference of about \SI{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 \SI{0.7}{\eV} for the {\sBoneg} 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 \SIrange{0.17}{0.2}{\eV} for the {\tBoneg} and the {\twoAg} states and for the different bases.
Again, the energy difference for the {\twoAg} state is larger with \SIrange{0.5}{0.7}{\eV} depending on the basis.
In a similar way than with XMS-CASPT2(4,4), XMS-CASPT(12,12) only describes the {\twoAg} 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 {\twoAg} states the energies are similar with an energy difference between the SC-NEVPT2(12,12) and the PC-NEVPT2(12,12) values of about \SIrange{0.03}{0.04}{\eV} and \SIrange{0.02}{0.03}{\eV} respectively.
Figure \ref{fig:D2h} shows the vertical energies of the studied excited states described in Tables \ref{tab:sf_tddft_D2h} and \ref{tab:D2h}.
We see that the SF-TD-DFT vertical energies are rather close to the TBE one for the doubly excited state {\twoAg}.
Indeed, as presented later in subsection \ref{sec:TBE}, SF methods are able to describe effectively states with strong multi-configurational character.
The same observation can be done for the SF-ADC values but with much better results for the two other states.
Again we can notice that the SF-ADC(2)-x and SF-ADC(3) schemes do not improve the values.
For the multi-reference methods we see that using the smallest active space do not provide a good description of the {\sBoneg} state and in the case of CASSCF(4,4), as previously said, we even have {\sBoneg} state above the {\twoAg} state.
For the CC methods taking into account the quadruple excitations is necessary to have a good description of the doubly excited state {\twoAg}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{{\Dfour} symmetry}
\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.}
\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
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 & & $\left[1.511 \right]$\fnm[1] &$\left[2.021 \right]$\fnm[1] \\[0.1cm]
CCSDTQ & 6-31+G(d) & $0.205$ & $1.593$ & $2.134$ \\
& aug-cc-pVDZ & $\left[0.160\right]$\fnm[2] & $\left[1.528 \right]$\fnm[4]&$\left[1.947\right]$\fnm[4] \\
& aug-cc-pVTZ & $\left[0.144\right]$\fnm[3] & $\left[1.500 \right]$\fnm[5]&$\left[2.034\right]$\fnm[5] \\ [0.1cm]
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]
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]{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]{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.}
%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*}
%%% %%% %%% %%%
\alert{In the {\Dtwo} symmetry we could not compute $\%T_1$ values associated with the {\Atwog}, {\Aoneg}, and {\Btwog} excited states.
However, it is clear from the inspection of the wave function that}
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 {\Atwog} and the ground state {\sBoneg}.
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 \SIrange{0.004}{0.007}{\eV} for the triplet state {\Atwog}.
We have \SIrange{0.015}{0.021}{\eV} of energy difference for the {\Aoneg} state through all bases, we can notice that this state is around \SI{0.13}{\eV} (considering all bases) higher with the PBE0 functional.
We can make the same observation for the {\Btwog} state for the B3LYP and the PBE0 functionals, indeed we have small energy differences for all bases and the state is around \SIrange{0.14}{0.15}{\eV} for the PBE0 functional.
For the BH\&HLYP functional the {\Aoneg} and {\Btwog} states are higher in energy than for the two other hybrid functionals with about \SIrange{0.65}{0.69}{\eV} higher for the {\Aoneg} state and \SIrange{0.75}{0.77}{\eV} for the {\Btwog} 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 {\Aoneg} and {\Btwog} states with a maximum energy difference of \SIrange{0.01}{0.02}{\eV} for the{\Aoneg} state and \SIrange{0.005}{0.009}{\eV} for the {\Btwog} state considering all bases.
The maximum energy difference for the triplet state is larger with \SIrange{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 {\Atwog} and the ground state {\sBoneg}.
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 {\Aoneg} state are close to the BH\&HLYP energies for the {\Btwog} 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 \SIrange{0.16}{0.17}{\eV} for the {\Aoneg} state and \SIrange{0.17}{0.18}{\eV} for the {\Btwog} state considering all bases.
For the triplet state {\Atwog} the energy differences are smaller with \SIrange{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 \SI{0.003}{\eV} considering all bases, and are closer to the BH\&HLYP results for the two other states with \SIrange{0.06}{0.07}{\eV} and \SIrange{0.07}{0.08}{\eV} of energy difference for the {\Aoneg} and {\Btwog} states, respectively.
Then we discuss the various ADC scheme results, note that we were not able to obtain the vertical energies with the aug-cc-pVQZ 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 \SI{0.09}{\eV} for the triplet state whereas we have \SI{0.15}{\eV} and \SI{0.25}{eV} for the {\Aoneg} and {\Btwog} 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 {\Aoneg} state between ADC(2)-s and ADC(2)-x with around \SIrange{0.52}{0.58}{\eV} through all bases.
The ADC(3) vertical energies are very similar to the ADC(2) ones for the {\Btwog} state with an energy difference of \SIrange{0.01}{0.02}{\eV} for all bases, whereas we have an energy difference of \SIrange{0.04}{0.11}{\eV} and \SIrange{0.17}{0.22}{\eV} for the {\Aoneg} and {\Btwog} states, respectively.
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 aug-cc-pVTZ 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.
Note that for CC we started from a restricted Hartree-Fock (RHF) reference and in that case the ground state is the {\Aoneg}state, then the {\sBoneg} state is the single deexcitation and the {\Btwog} state is the double excitation from our ground state.
For the CC3 method we do not have the vertical energies for the triplet state {\Atwog}.
Considering all bases for the {\Aoneg}and {\Btwog} states we have an energy difference of about \SI{0.15}{\eV} and \SI{0.12}{\eV}, respectively.
The CCSDT energies are close to the CC3 ones for the {\Aoneg} state with an energy difference of around \SIrange{0.03}{0.06}{\eV} considering all bases.
For the{\Btwog} state the energy difference between the CC3 and the CCSDT values is larger with \SIrange{0.18}{0.27}{\eV}.
We can make a similar observation between the CC4 and the CCSDTQ values, for the {\Aoneg} state we have an energy difference of about \SI{0.01}{\eV} and this time we have smaller energy difference for the {\Btwog} with \SI{0.01}{\eV}.
Then we discuss the multi-reference results and this time we were able to reach the aug-cc-pVQZ 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 {\Aoneg} and {\Btwog} states, for the {\Aoneg} state we have an energy difference of about \SIrange{0.67}{0.74}{\eV} and \SIrange{1.65}{1.81}{\eV} for the {\Btwog} state.
The energy difference is smaller for the triplet state with \SIrange{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 {\Aoneg} and {\Btwog} states with {\Btwog} higher in energy than {\Aoneg} 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 \SI{0.06}{\eV} for all bases but larger energy difference for the {\Aoneg} state with around \SIrange{0.28}{0.29}{\eV} and \SIrange{0.79}{0.81}{\eV} for the {\Btwog} 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 {\Aoneg} state, considering all bases, with an energy difference of around \SIrange{0.05}{0.06}{\eV} and \SIrange{0.02}{0.05}{\eV} respectively.
The energy difference is larger for the {\Btwog} state with about \SIrange{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 {\Aoneg} and {\Btwog} states.
Again, every vertical energy is higher in the NEVPT2(12,12) case than for the NEVPT2(4,4) one.
Figure \ref{fig:D4h} shows the vertical energies of the studied excited states described in Tables \ref{tab:sf_D4h} and \ref{tab:D4h}.
We see that all methods exhibit larger variations of the vertical energies due to the high symmetry of the {\Dfour} structure.
For the triplet state most of the methods used are able to give a vertical energy close to the TBE one, we can although see that CASSCF, with the two different active spaces, show larger energy error than most of the DFT functionals used.
Then for the, strongly multi-configurational character, {\Aoneg} state we have a good description by the CC and multi-reference methods with the largest active space, except for CASSCF.
The SF-ADC(2) and SF-ADC(3) schemes are also able to provide a good description of the {\Aoneg} state and even for the {\Btwog} we see that SF-ADC(2)-x give a worst result than SF-ADC(2).
Multi-configurational methods with the smallest active space do not demonstrate a good description of the two singlet excited states especially for the CASSCF method and for the the NEVPT2 methods that give the vertical energy of the {\Btwog} state below the {\Aoneg} one.
Note that CASPT2 improve a lot the description of all the states compared to CASSCF.
The various TD-DFT functionals are not able to describe correctly the two singlet excited states.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Theoretical best estimates}
\label{sec:TBE}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{squeezetable}
\begin{table*}
\caption{Energy differences between the various methods and the reference TBE values.
Note that AB stands for the autoisomerization barrier and is reported in \si{\kcalmol}.
The numbers reported in parenthesis are the percentage of single excitations involved in the transition ($\%T_1$) calculated at the CC3/aug-cc-pVTZ level.
The values between square brackets have been obtained by extrapolation via the procedure described in the corresponding footnote.}
\label{tab:TBE}
\begin{ruledtabular}
\begin{tabular}{lrrrrrrr}
%\begin{tabular}{*{1}{*{8}{l}}}
& &\mc{3}{c}{{\Dtwo} excitation energies (eV)} & \mc{3}{c}{{\Dfour} excitation energies (eV)} \\
\cline{3-5} \cline{6-8}
Method & AB & {\tBoneg}(99\%) & {\sBoneg}(95\%)& {\twoAg}(1\%) & {\Atwog} & {\Aoneg} & {\Btwog} \\
\hline
SF-TD-B3LYP & $10.41$ & $0.241$ & $-0.926$ & $-0.161$ & $-0.164$ & $-1.028$ & $-1.501$ \\
SF-TD-PBE0 & $8.95$ & $0.220$ & $-0.829$ & $-0.068$ & $-0.163$ & $-0.903$ & $-1.357$ \\
SF-TD-BHHLYP & $3.79$ & $0.078$ & $-0.393$ & $0.343$ & $-0.099$ & $-0.251$ & $-0.603$ \\
SF-TD-M06-2X & $1.42$ & $0.000$ & $-0.354$ & $0.208$ & $-0.066$ & $-0.097$ & $-0.432$ \\
SF-TD-CAM-B3LYP & $9.90$ & $0.280$ & $-0.807$ & $-0.011$ & $-0.134$ & $-0.920$ & $-1.370$ \\
SF-TD-$\omega $B97X-V & $10.01$ & $0.335$ & $-0.774$ & $0.064$ & $-0.118$ & $-0.928$ & $-1.372$ \\
SF-TD-M11 & $2.29$ & $0.097$ & $-0.474$ & $0.151$ & $-0.063$ & $-0.312$ & $-0.675$ \\
SF-TD-LC-$\omega $PBE08 & $10.81$ & $0.435$ & $-0.710$ & $0.199$ & $-0.086$ & $-0.939$ & $-1.376$ \\[0.1cm]
SF-ADC(2)-s & $-0.30$ & $0.069$ & $-0.026$ & $-0.018$ & $0.112$ & $0.112$ & $-0.190$ \\
SF-ADC(2)-x & $1.44$ & $0.077$ & $-0.094$ & $-0.446$ & $0.068$ & $-0.409$ & $-0.303$ \\
SF-ADC(2.5) & $0.18$ & $0.013$ & $0.006$ & $0.029$ & $0.024$ & $0.094$ & $-0.185$ \\
SF-ADC(3) & $0.65$ & $-0.043$ & $0.037$ & $0.075$ & $-0.065$ & $0.075$ & $-0.181$ \\[0.1cm]
CCSD & $0.95$ & & & & $-0.059$ & $0.100$ & \\
CC3 & $-1.05$ & $-0.060$ & $-0.006$ & $0.628$ & & $0.162$ & $0.686$ \\
CCSDT & $-0.25$ & $-0.051$ & $0.014$ & $0.280$ & $0.005$ & $0.131$ & $0.503$ \\
CC4 & $-0.11$ & & $0.003$ & $-0.006$ & & $0.011$ & $-0.013$ \\
CCSDTQ & $0.00$ & & $0.000$ & $0.000$ & $0.000$ & $0.000$ & $0.000$ \\[0.1cm]
CASSCF(4,4) & $-1.55$ & $0.208$ & $1.421$ & $0.292$ & $0.290$ & $0.734$ & $1.390$ \\
CASPT2(4,4) & $-1.16$ & $-0.050$ & $-0.202$ & $-0.077$ & $-0.016$ & $0.006$ & $-0.399$ \\
%XMS-CASPT2(4,4) & & & & $-0.035$ & & & \\
SC-NEVPT2(4,4) & $0.30$ & $-0.083$ & $-0.703$ & $-0.041$ & $-0.120$ & $-0.072$ & $-0.979$ \\
PC-NEVPT2(4,4) & $0.31$ & $-0.080$ & $-0.757$ & $-0.066$ & $-0.118$ & $-0.097$ & $-1.031$ \\
%MRCI(4,4) & & $0.106$ & $0.553$ & $0.121$ & $0.127$ & $0.324$ & $0.381$ \\[0.1cm]
CASSCF(12,12) & $2.66$ & $0.224$ & $0.719$ & $0.068$ & $0.226$ & $0.443$ & $0.600$ \\
CASPT2(12,12) & $-0.42$& $0.018$ & $0.058$ & $-0.106$ & $0.039$ & $0.038$ & $-0.108$ \\
%XMS-CASPT2(12,12) & & & & $-0.090$ & & & \\
SC-NEVPT2(12,12) & $-0.64$ & $0.039$ & $0.063$ & $-0.063$ & $0.021$ & $0.046$ & $-0.142$ \\
PC-NEVPT2(12,12) & $-0.65$ & $0.000$ & $-0.062$ & $-0.093$ & $-0.013$ & $-0.024$ & $-0.278$ \\[0.1cm]
%CIPSI & & $-0.001\pm 0.030$ & $0.017\pm 0.035$ & $-0.120\pm 0.090$ & $0.025\pm 0.029$ & $0.130\pm 0.050$ & \\
\bf{TBE} & $[\bf{8.93}]$\fnm[1] & $[\bf{1.462}]$\fnm[2] & $[\bf{3.125}]$\fnm[3] & $[\bf{4.149}]$\fnm[3] & $[\bf{0.144}]$\fnm[4] & $[\bf{1.500}]$\fnm[3] & $[\bf{2.034}]$\fnm[3] \\[0.1cm]
Literature & $8.53$\fnm[5] & $1.573$\fnm[5] & $3.208$\fnm[5] & $4.247$\fnm[5] & $0.266$\fnm[5] & $1.664$\fnm[5] & $1.910$\fnm[5] \\
& $10.35$\fnm[6] & $1.576$\fnm[6] & $3.141$\fnm[6] & $3.796$\fnm[6] & $0.217$\fnm[6] & $1.123$\fnm[6] & $1.799$\fnm[6]\\
& $9.58$ \fnm[7]& $1.456$\fnm[7] & $3.285$\fnm[7] & $4.334$\fnm[7] & $0.083$\fnm[7] & $1.621$\fnm[7] & $1.930$\fnm[7] \\
& $7.48$\fnm[8]& $1.654$\fnm[8] & $3.416$\fnm[8] & $4.360$\fnm[8] & $0.369$\fnm[8] & $1.824$\fnm[8] & $2.143$\fnm[8] \\
\end{tabular}
\end{ruledtabular}
\fnt[1]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
\fnt[2]{Value obtained using the NEVPT2(12,12) one.}
\fnt[3]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
\fnt[4]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
\fnt[5]{Value obtained from Ref.~\onlinecite{Lefrancois_2015} at the SF-ADC(2)-s/cc-pVTZ level with the geometry obtained at the CCSD(T)/cc-pVTZ level.}
\fnt[6]{Value obtained from Ref.~\onlinecite{Lefrancois_2015} at the SF-ADC(2)-x/cc-pVTZ level with the geometry obtained at the CCSD(T)/cc-pVTZ level.}
\fnt[7]{Value obtained from Ref.~\onlinecite{Lefrancois_2015} at the SF-ADC(3)/cc-pVTZ level with the geometry obtained at the CCSD(T)/cc-pVTZ level.}
\fnt[8]{Value obtained from Ref.~\onlinecite{Manohar_2008} at the EOM-SF-CCSD/cc-pVTZ level with the geometry obtained at the CCSD(T)/cc-pVTZ level.}
\end{table*}
\end{squeezetable}
%================================================
Table \ref{tab:TBE} shows the energy differences, for the autoisomerization barrier (AB) and the different states considered, between the various methods treated and the Theoretical Best Estimate (TBE) at the aug-cc-pVTZ level for the AB and the states.
The percentage $\%T_1$ shown in parentheses for the excited states of the {\Dtwo} geometry is a metric that gives the percentage of single excitation calculated at the CC3/aug-cc-pVTZ level and it allows us to characterize the transition.
First, we look at the AB energy difference.
SF-TD-DFT shows large variations of the energy with errors of \SIrange{1.42}{10.81}{\kcalmol} compared to the TBE value.
SF-ADC schemes provide smaller errors with \SIrange{0.30}{1.44}{\kcalmol} where we have that the SF-ADC(2)-x gives a worse error than the SF-ADC(2)-s method.
CC methods also give small energy differences with \SIrange{0.11}{1.05}{\kcalmol} and where the CC4 provides an energy very close to the TBE one.
Then we look at the vertical energy errors for the {\Dtwo} structure.
First we consider the {\tBoneg} state and we look at the SF-TD-DFT results.
We see that increasing the amount of exact exchange in the functional give closer results to the TBE, indeed we have \SI{0.24}{\eV} and \SI{0.22}{\eV} of errors for the B3LYP and the PBE0 functionals, respectively whereas we have an error of \SI{0.08}{\eV} for the BH\&HLYP functional.
For the other functionals we have errors of \SIrange{0.10}{0.43}{\eV}, note that for this state the M06-2X functional gives the same result than the TBE.
We can also notice that all the functionals considered overestimate the vertical energies.
The ADC schemes give closer energies with errors of \SIrange{0.04}{0.08}{\eV}, note that ADC(2)-x does not improve the result compared to ADC(2)-s and that ADC(3) understimate the vertical energy whereas ADC(2)-s and ADC(2)-x overestimate the vertical energy.
The CC3 and CCSDT results provide energy errors of \SIrange{0.05}{0.06}{\eV} respectively.
Then we go through the multi-reference methods with the two different active spaces, four electrons in four orbitals and twelve electrons in twelve orbitals.
For the smaller active space we have errors of \SIrange{0.05}{0.21}{\eV}, the largest error comes from CASSCF(4,4) which is improved by CASPT2(4,4) that gives the smaller error.
Then for the largest active space multi-reference methods provide energy errors of \SIrange{0.02}{0.22}{\eV} with again the largest error coming from CASSCF(12,12) which is again improved by CASPT2(12,12) gives the smaller error.
For the {\sBoneg} state of the {\Dtwo} structure we see that all the xc-functional underestimate the vertical excitation energy with energy differences of about \SIrange{0.35}{0.93}{\eV}.
The ADC values are much closer to the TBE with energy differences around \SIrange{0.03}{0.09}{\eV}.
Obviously, the CC vertical energies are close to the TBE one with around or less than \SI{0.01}{\eV} of energy difference.
For the CASSCF(4,4) vertical energy we have a large difference of around \SI{1.42}{\eV} compared to the TBE value due to the lack of dynamical correlation in the CASSCF method.
As previously seen the CAPT2(4,4) method correct this and we obtain a value of \SI{0.20}{\eV}.
The others multi-reference methods in this active space give energy differences of around \SIrange{0.55}{0.76}{\eV} compared the the TBE reference.
For the largest active space with twelve electrons in twelve orbitals we have an improvement of the vertical energies with \SI{0.72}{eV} of energy difference for the CASSCF(12,12) method and around 0.06 eV for the others multi-configurational methods.
Then, for the {\twoAg} state we obtain closer results of the SF-TD-DFT methods to the TBE than in the case of the {\sBoneg} state.
Indeed, we have an energy difference of about \SIrange{0.01}{0.34}{\eV} for the {\twoAg} state whereas we have \SIrange{0.35}{0.93}{\eV} for the {\sBoneg} state.
The ADC schemes give the same error to the TBE value than for the other singlet state with \SI{0.02}{\eV} for the ADC(2) scheme and \SI{0.07}{\eV} for the ADC(3) one.
The ADC(2)-x scheme provides a larger error with \SI{0.45}{\eV} of energy difference.
Here, the CC methods manifest more variations with \SI{0.63}{\eV} for the CC3 value and \SI{0.28}{\eV} for the CCSDT compared to the TBE values.
The CC4 method provides a small error with less than 0.01 eV of energy difference.
The multi-configurational methods globally give smaller error than for the other singlet state with, for the two active spaces, \SIrange{0.03}{0.12}{\eV} compared to the TBE value.
We can notice that CC3 and CCSDT provide larger energy errors for the {\twoAg} state than for the {\sBoneg} state, this is due to the strong multi-configurational character of the {\twoAg} state whereas the {\sBoneg} state has a very weak multi-configurational character.
It is interesting to see that SF-TD-DFT and SF-ADC methods give better results compared to the TBE than CC3 and even CCSDT meaning that spin-flip methods are able to describe double excited states.
Note that multi-reference methods obviously give better results too for the {\twoAg} state.
Finally we look at the vertical energy errors for the \Dfour structure.
First, we consider the {\Atwog} state, the SF-TD-DFT methods give errors of about \SIrange{0.07}{1.6}{\eV} where the largest energy differences are provided by the hybrid functionals.
The ADC schemes give similar error with around \SIrange{0.06}{1.1}{\eV} of energy difference.
For the CC methods we have an energy error of \SI{0.06}{\eV} for CCSD and less than \SI{0.01}{\eV} for CCSDT.
Then for the multi-reference methods with the four by four active space we have for CASSCF(4,4) \SI{0.29}{\eV} of error and \SI{0.02}{\eV} for CASPT2(4,4), again CASPT2 demonstrates its improvement compared to CASSCF.
The other methods provide energy differences of about \SIrange{0.12}{0.13}{\eV}.
A larger active space shows again an improvement with \SI{0.23}{\eV} of error for CASSCF(12,12) and around \SIrange{0.01}{0.04}{\eV} for the other multi-reference methods.
CIPSI provides similar error with \SI{0.02}{\eV}.
Then, we look at the {\Aoneg} state where the SF-TD-DFT shows large variations of error depending on the functionals, the energy error is about \SIrange{0.10}{1.03}{\eV}.
The ADC schemes give better errors with around \SIrange{0.07}{0.41}{\eV} and where again the ADC(2)-x does not improve the ADC(2)-s energy but also gives worse results.
For the CC methods we have energy errors of about \SIrange{0.10}{0.16}{\eV} and CC4 provides really close energy to the TBE one with \SI{0.01}{\eV} of error.
For the multi-reference methods we globally have an improvement of the energies from the four by four to the twelve by twelve active space with errors of \SIrange{0.01}{0.73}{\eV} and \SIrange{0.02}{0.44}{\eV} respectively with the largest errors coming from the CASSCF method.
Lastly, we look at the {\Btwog} state where we have globally larger errors.
The SF-TD-DFT exhibits errors of \SIrange{0.43}{1.50}{\eV} whereas ADC schemes give errors of \SIrange{0.18}{0.30}{\eV}.
CC3 and CCSDT provide energy differences of \SIrange{0.50}{0.69}{\eV} and the CC4 shows again close energy to the CCSDTQ TBE energy with \SI{0.01}{\eV} of error.
The multi-reference methods give energy differences of \SIrange{0.38}{1.39}{\eV} and \SIrange{0.11}{0.60}{\eV} for the four by four and twelve by twelve active spaces respectively.
We can notice that we have larger variations for the vertical energies of the square CBD than for the rectangular one.
This can be explained by the fact that because of the degeneracy in the {\Dfour} structure it leads to strong multi-configurational character states where single reference methods are unreliable.
We can also see that for the CC methods we have a better description of the {\Aoneg} state than the {\Btwog} state, this can be related, as previously said in Subsubsection \ref{sec:D4h}, to the fact that {\Btwog} corresponds to a double excitation from the reference state.
To obtain an improved description of the {\Btwog} state we have to include quadruples.
At the end of Table \ref{tab:TBE} we show some literature results obtain from Ref.~\onlinecite{Lefrancois_2015,Manohar_2008} where the cc-pVTZ basis is used.
The SF-ADC(2)-s, SF-ADC(2)-x and SF-ADC(3)results are presented and are consistent with our results with the exact same schemes but with the aug-cc-pVTZ basis.
%Here again we can make the same comment for the $2\,{}^1A_{1g}$ and $1\,{}^1B_{2g}$ states of the square CBD than the $1\,{}^1B_{1g}$ and $2\,{}^1A_{g}$ states of the rectangular CBD. The first state ($2\,{}^1A_{1g}$) has a strong multi-configurational character
%%% TABLE I %%%
%\begin{squeezetable}
%\begin{table*}
% \caption{Energy differences between the various methods and the TBE considered. Note that AB stands for the autoisomerization barrier and that the energies are in \kcalmol while the vertical energies are given in eV. The number in parenthesis is the percentage $\%T_1$ calculated at the CC3/aug-cc-pVTZ level.}
%
% \label{tab:TBE}
% \begin{ruledtabular}
% \begin{tabular}{lrrrrrrr}
%%\begin{tabular}{*{1}{*{8}{l}}}
%&\mc{3}{r}{\Dtwo excitation energies (eV)} & \mc{3}{r}{\Dfour excitation energies (eV)} \\
% \cline{3-5} \cline{6-8}
%Method & AB & $1\,{}^3B_{1g} $~(98.7 \%) & $1\,{}^1B_{1g} $ (95.0 \%)& $2\,{}^1A_{g} $(0.84 \%) & $1\,{}^3A_{2g} $ & $2\,{}^1A_{1g} $ & $1\,{}^1B_{2g} $ \\
% \hline
%SF-TD-B3LYP & $10.41$ & $0.241$ & $-0.926$ & $-0.161$ & $-0.164$ & $-1.028$ & $-1.501$ \\
%SF-TD-PBE0 & $8.95$ & $0.220$ & $-0.829$ & $-0.068$ & $-0.163$ & $-0.903$ & $-1.357$ \\
%SF-TD-BHHLYP & $8.95$ & $0.078$ & $-0.393$ & $0.343$ & $-0.099$ & $-0.251$ & $-0.603$ \\
%SF-TD-M06-2X & $1.42$ & $0.000$ & $-0.354$ & $0.208$ & $-0.066$ & $-0.097$ & $-0.432$ \\
%SF-TD-CAM-B3LYP & $9.90$ & $0.280$ & $-0.807$ & $-0.011$ & $-0.134$ & $-0.920$ & $-1.370$ \\
%SF-TD-$\omega $B97X-V & $10.01$ & $0.335$ & $-0.774$ & $0.064$ & $-0.118$ & $-0.928$ & $-1.372$ \\
%SF-TD-M11 & $2.29$ & $0.097$ & $-0.474$ & $0.151$ & $-0.063$ & $-0.312$ & $-0.675$ \\
%SF-TD-LC-$\omega $PBE08 & $10.81$ & $0.435$ & $-0.710$ & $0.199$ & $-0.086$ & $-0.939$ & $-1.376$ \\[0.1cm]
%SF-ADC(2)-s & $-0.30$ & $0.069$ & $-0.026$ & $-0.018$ & $0.112$ & $0.112$ & $-0.190$ \\
%SF-ADC(2)-x & $1.44$ & $0.077$ & $-0.094$ & $-0.446$ & $0.068$ & $-0.409$ & $-0.303$ \\
%SF-ADC(3) & $0.65$ & $-0.043$ & $0.037$ & $0.075$ & $-0.065$ & $0.075$ & $-0.181$ \\[0.1cm]
%CCSD & $0.95$ & & & & $-0.059$ & $0.100$ & \\
%CC3 & $-1.05$ & $-0.060$ & $-0.006$ & $0.628$ & & $0.162$ & $0.686$ \\
%CCSDT & $-0.25$ & $-0.051$ & $0.014$ & $0.280$ & $0.005$ & $0.131$ & $0.503$ \\
%CC4 & $-0.11$ & & $0.003$ & $-0.006$ & & $0.011$ & $-0.013$ \\
%CCSDTQ & $0.00$ & & $0.000$ & $0.000$ & $0.000$ & $0.000$ & $0.000$ \\[0.1cm]
%CASSCF(4,4) & & $0.208$ & $1.421$ & $0.292$ & $0.290$ & $0.734$ & $1.390$ \\
%CASPT2(4,4) & & $-0.050$ & $-0.202$ & $-0.077$ & $-0.016$ & $0.006$ & $-0.399$ \\
%XMS-CASPT2(4,4) & & & & $-0.035$ & & & \\
%SC-NEVPT2(4,4) & & $-0.083$ & $-0.703$ & $-0.041$ & $-0.120$ & $-0.072$ & $-0.979$ \\
%PC-NEVPT2(4,4) & & $-0.080$ & $-0.757$ & $-0.066$ & $-0.118$ & $-0.097$ & $-1.031$ \\
%MRCI(4,4) & & $0.106$ & $0.553$ & $0.121$ & $0.127$ & $0.324$ & $0.381$ \\[0.1cm]
%CASSCF(12,12) & $2.66$ & $0.224$ & $0.719$ & $0.068$ & $0.226$ & $0.443$ & $0.600$ \\
%CASPT2(12,12) & & $0.018$ & $0.058$ & $-0.106$ & $0.039$ & $0.038$ & $-0.108$ \\
%XMS-CASPT2(12,12) & & & & $-0.090$ & & & \\
%SC-NEVPT2(12,12) & $-0.65$ & $0.039$ & $0.063$ & $-0.063$ & $0.021$ & $0.046$ & $-0.142$ \\
%PC-NEVPT2(12,12) & & $0.000$ & $-0.062$ & $-0.093$ & $-0.013$ & $-0.024$ & $-0.278$ \\[0.1cm]
%%CIPSI & & $-0.001\pm 0.030$ & $0.017\pm 0.035$ & $-0.120\pm 0.090$ & $0.025\pm 0.029$ & $0.130\pm 0.050$ & \\
%\bf{TBE} & $\left[\bf{8.93}\right]$\fnm[1] & $\left[\bf{1.462}\right]$\fnm[2] & $\left[\bf{3.125}\right]$\fnm[3] & $\left[\bf{4.149}\right]$\fnm[3] & $\left[\bf{0.144}\right]$\fnm[4] & $\left[\bf{1.500}\right]$\fnm[3] & $\left[\bf{2.034}\right]$\fnm[3] \\
%\end{tabular}
%
% \end{ruledtabular}
% \fnt[1]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
% \fnt[2]{Value obtained using the NEVPT2(12,12) one.}
% \fnt[3]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ.}
% \fnt[4]{Value obtained using CCSDTQ/aug-cc-pVDZ corrected by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ.}
%
%\end{table*}
%\end{squeezetable}
%%% %%% %%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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.
%Otherwise we got the CCSDTQ/aug-cc-pVTZ values by correcting the CCSDTQ/aug-cc-pVDZ values by the difference between CC4/aug-cc-pVTZ and CC4/aug-cc-pVDZ (Eq.~\eqref{eq:aug-cc-pVTZ}) and we obtain the CCSDTQ/aug-cc-pVDZ values by correcting the CCSDTQ/6-31G+(d) values by the difference between CC4/aug-cc-pVDZ and CC4/6-31G+(d) (Eq.~\eqref{eq:aug-cc-pVDZ}).
%When the CC4/aug-cc-pVTZ values were not obtained we corrected the CC4/aug-cc-pVDZ values by the difference between CCSDT/aug-cc-pVTZ and CCSDT/aug-cc-pVDZ to obtain them (Eq.~\eqref{eq:CC4_aug-cc-pVTZ}).
%If the CC4 values have not been obtained then we used the same scheme that we just described but by using the CCSDT values.
%If neither the CC4 and CCSDTQ values were not available then we used the NEVPT2(12,12)/aug-cc-pVTZ values.
In order to provide a benchmark of the automerization barrier and vertical energies we used coupled-cluster (CC) methods with doubles (CCSD), with triples (CCSDT and CC3) and with quadruples (CCSTQ and CC4).
Due to the presence of multi-configurational states we used multi-reference methods (CASSCF, CASPT2 and NEVPT2) with two active spaces ((4,4) and (12,12)).
We also used spin-flip (SF-) within two frameworks, in TD-DFT with various global and range-separated hybrids functionals, and in ADC with the ADC(2)-s, ADC(2)-x and ADC(3) schemes.
The CC methods provide good results for the AB and vertical energies, however in the case of multi-configurational states CC with only triples is not sufficient and we have to include the quadruples to correctly describe these states.
multi-configurational methods also provide very solid results for the largest active space with second order correction (CASPT2 and NEVPT2).
With SF-TD-DFT the quality of the results are, of course, dependent on the functional but for the doubly excited states we have solid results.
In SF-ADC we have very good results compared to the TBEs even for the doubly excited states, nevertheless the ADC(2)-x scheme give almost systematically worse results than the ADC(2)-s ones and using the ADC(3) scheme does not always provide better values.
The description of the excited states of the {\Dtwo} structure give rise to good agreement between the single reference and multi-configurational methods due to the large $\%T_1$ percentage of the first two excited states.
When this percentage is much smaller as in the case of the doubly excited state {\twoAg} the spin-flip methods show very good results within the ADC framework but even at the TD-DFT level spin-flip display good results compared to the TBE value.
As said in the discussion, for the {\Dfour} geometry, the description of excited states is harder because of the strong multi-configurational character where SF-TD-DFT can present more than 1 eV of error compared to the TBE.
However, SF-ADC can show error of around \SIrange{0.1}{0.2}{\eV} which can be better than the multi-configurational methods results.
%%%%%%%%%%%%%%%%%%%%%%%%
\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}