resub take 2

Manuscript/CBD.bib

@@ -1,13 +1,24 @@
 %% This BibTeX bibliography file was created using BibDesk.
 %% https://bibdesk.sourceforge.io/
 
-%% Created for Pierre-Francois Loos at 2022-06-08 18:49:54 +0200 
+%% Created for Pierre-Francois Loos at 2022-06-30 10:43:59 +0200 
 
 
 %% Saved with string encoding Unicode (UTF-8) 
 
 
 
+@article{Mulliken_1955,
+	author = {R. S. Mulliken},
+	date-added = {2022-06-30 10:42:34 +0200},
+	date-modified = {2022-06-30 10:43:54 +0200},
+	doi = {10.1063/1.1740655},
+	journal = {J. Chem. Phys.},
+	pages = {1997--2011},
+	title = {Report on Notation for the Spectra of Polyatomic Molecules},
+	volume = {23},
+	year = {1955}}
+
 @article{Knowles_1988,
 	abstract = {A new method for evaluating one-particle coupling coefficients in a general configuration interaction calculation is presented. Through repeated application and use of resolutions of the identity, two-, three- and four-body coupling coefficients and density matrices may be built in a simple and efficient way. The method is therefore of use in both multiconfiguration SCF (MC SCF) and multireference configuration interaction (MRCI) calculations. Examples show that the approach is efficient for both these applications.},
 	author = {Peter J. Knowles and Hans-Joachim Werner},

Manuscript/CBD.tex

@@ -13,12 +13,13 @@
 \DeclareSIUnit\kcal{\kilo\cal}
 \newcommand{\kcalmol}{\si{\kcal\per\mole}}
 
-\usepackage[
-	colorlinks=true,
-    citecolor=blue,
-    breaklinks=true
-	]{hyperref}
-\urlstyle{same}
+%\usepackage[
+%	colorlinks=false,
+%    citecolor=black,
+%    breaklinks=false
+%	]{hyperref}
+%\urlstyle{same}
+\renewcommand{\thetable}{\arabic{table}} 
 
 \newcommand{\ie}{\textit{i.e.}}
 \newcommand{\eg}{\textit{e.g.}}
@@ -33,7 +34,7 @@
 \newcommand{\fnm}{\footnotemark}
 \newcommand{\fnt}{\footnotetext}
 \newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
-\newcommand{\SupInf}{\textcolor{blue}{supporting information}}
+\newcommand{\SupInf}{\textcolor{black}{supporting information}}
 \newcommand{\QP}{\textsc{quantum package}}
 \newcommand{\T}[1]{#1^{\intercal}}
 
@@ -86,16 +87,11 @@
 \begin{abstract}
 \paragraph*{Abstract:}	
 Cyclobutadiene 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 \alert{standard} adiabatic time-dependent density-functional theory (TD-DFT) or \alert{standard} equation-of-motion coupled cluster (EOM-CC), are notoriously known to struggle in such situations.
+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 standard adiabatic time-dependent density-functional theory (TD-DFT) or standard equation-of-motion coupled cluster (EOM-CC), are notoriously known to struggle in such situations.
 In this work, using a large panel of methods and basis sets, we provide an extensive computational study of the automerization barrier (defined as the difference between the square and rectangular ground-state energies) and the vertical excitation energies at $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 qualitatively correct description of these diradical states, is also tested within TD-DFT (combined with numerous exchange-correlation functionals) and the algebraic diagrammatic construction [ADC(2)-s, ADC(2)-x, and ADC(3)] schemes. 
 A theoretical best estimate is defined for the automerization barrier and for each vertical transition energy.
-\bigskip
-\begin{center}
-        \boxed{\includegraphics[width=0.4\linewidth]{TOC}}
-\end{center}
-\bigskip
 \end{abstract}
 
 \maketitle
@@ -117,31 +113,31 @@ This was confirmed by several experimental studies by Pettis and co-workers \cit
 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 is a diradical that has two degenerate singly occupied frontier orbitals.
 Therefore, one must take into account, at least, two electronic configurations to properly model this multi-configurational scenario. 
-Of course, \alert{standard} single-reference methods are naturally unable to describe such situations. 
+Of course, standard single-reference methods are naturally unable to describe such situations. 
 Interestingly, the {\sBoneg} ground state 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 in the {\Dfour} arrangement. 
 The automerization barrier (AB) is thus defined as the difference between the square and rectangular ground-state energies.
 The energy of this barrier is estimated, experimentally, in the range of \SIrange{1.6}{10}{\kcalmol}, \cite{Whitman_1982} while previous state-of-the-art \textit{ab initio} calculations yield values in the \SIrange{7}{9}{\kcalmol} range. \cite{Eckert-Maksic_2006,Li_2009,Shen_2012,Zhang_2019}
 
 The lowest-energy excited states of CBD in both symmetries are represented in Fig.~\ref{fig:CBD}, where we have reported the {\oneAg} and {\tBoneg} states for the rectangular geometry and the {\sBoneg} and {\Atwog} states for the square one. 
 Due to the energy scale, the higher-energy states ({\sBoneg} and {\twoAg} for {\Dtwo} and {\Aoneg} and {\Btwog} for {\Dfour}) are not shown. 
-Interestingly, the {\twoAg} and {\Aoneg} states have a strong contribution from doubly-excited configurations and these so-called double excitations \cite{Loos_2019} are known to be inaccessible with \alert{standard} 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 \alert{remain challenging for standard hierarchy of EOM-CC methods that are using ground-state Hartree-Fock reference}. \cite{Kucharski_1991,Kallay_2004,Hirata_2000,Hirata_2004}
+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 inaccessible with standard 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 remain challenging for standard hierarchy of EOM-CC methods that are using ground-state Hartree-Fock reference. \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 approaches. 
-The most evident way is to rely on \alert{multi-reference} methods, which are naturally designed to address such scenarios.
+The most evident way is to rely on multi-reference 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}
 
 Another way to deal with double excitations and multi-reference situations is to use high level truncation of the EOM formalism \cite{Rowe_1968,Stanton_1993} of CC theory. \cite{Kucharski_1991,Kallay_2003,Kallay_2004,Hirata_2000,Hirata_2004}
 However, to provide a correct description of these situations, one has 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}
 Although multi-reference CC methods have been designed, \cite{Jeziorski_1981,Mahapatra_1998,Mahapatra_1999,Lyakh_2012,Kohn_2013} they are computationally demanding and remain far from being black-box.
 
-In this context, an interesting alternative to \alert{multi-reference} 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 are able to provide near full CI (FCI) ground- and excited-state energies of small molecules. \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}
+In this context, an interesting alternative to multi-reference 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 are able to provide near full CI (FCI) ground- and excited-state energies of small molecules. \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 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) de-excitation and excitation from the lowest triplet state, respectively.
-\alert{One drawback of spin-flip methods is spin contamination} (\ie, the artificial mixing of electronic states with different spin multiplicities) due not only to the spin incompleteness in the spin-flip expansion but also to the potential spin contamination of the reference configuration. \cite{Casanova_2020} 
+One drawback of spin-flip methods is spin contamination (\ie, the artificial mixing of electronic states with different spin multiplicities) due not only to the spin incompleteness in the spin-flip expansion but also to the potential 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}
-\alert{Note that one can quantify the polyradical character associated to a given electronic state using Head-Gordon's index \cite{Head-Gordon_2003} that provides a measure of the number of unpaired electrons. \cite{Orms_2018}}
+Note that one can quantify the polyradical character associated to a given electronic state using Head-Gordon's index \cite{Head-Gordon_2003} that provides a measure of the number of unpaired electrons. \cite{Orms_2018}
 
 In the present work, we define highly-accurate reference values and investigate the accuracy of each family of computational methods mentioned above on the automerization barrier and the low-lying excited states of CBD at the {\Dtwo} and {\Dfour} ground-state geometries. 
 Computational details are reported in Sec.~\ref{sec:compmet}.
@@ -190,13 +186,13 @@ 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 (except for CCSD). 
 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 at the CCSDT level with MRCC. \cite{mrcc} 
 
-To avoid having to perform multi-reference CC calculations or 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} \alert{singlet state of {$A_g$} symmetry} as reference. 
+To avoid having to perform multi-reference CC calculations or 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 of {$A_g$} symmetry as reference. 
 Hence, the open-shell ground state, {\sBoneg}, and the {\Btwog} state appear as a de-excitation and an excitation, respectively.
-With respect to \alert{this closed-shell reference}, {\sBoneg} has a dominant double excitation character, while {\Btwog} has a dominant single excitation character, hence their contrasting convergence behaviors with respect to the order of the CC expansion (see below). 
+With respect to this closed-shell reference, {\sBoneg} has a dominant double excitation character, while {\Btwog} has a dominant single excitation character, hence their contrasting convergence behaviors with respect to the order of the CC expansion (see below). 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-\subsection{\alert{Multi-reference} calculations}
+\subsection{Multi-reference calculations}
 \label{sec:Multi}
 State-averaged CASSCF (SA-CASSCF) calculations are performed for vertical transition energies, whereas state-specific CASSCF is used for computing the automerization barrier. \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).
@@ -206,7 +202,7 @@ For ionic excited states, like the {\sBoneg} state of CBD, it is particularly im
 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 intruder 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.25}{\hartree} to avoid systematic underestimation of the vertical excitation energies. \cite{Ghigo_2004,Schapiro_2013,Zobel_2017,Sarkar_2022}
-\alert{For the sake of comparison and completeness, for the (4e,4o) active space, we also report (in the {\SupInf}) multi-reference CI calculations including Davidson correction (MRCI+Q). \cite{Knowles_1988,Werner_1988}}
+For the sake of comparison and completeness, for the (4e,4o) active space, we also report (in the {\SupInf}) multi-reference CI calculations including Davidson correction (MRCI+Q). \cite{Knowles_1988,Werner_1988}
 All these calculations are carried out with MOLPRO. \cite{Werner_2020} 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -222,16 +218,21 @@ Here, we explore the spin-flip version \cite{Lefrancois_2015} of the algebraic-d
 These calculations are performed using Q-CHEM 5.4.1. \cite{qchem} 
 The spin-flip version of our recently proposed composite approach, namely SF-ADC(2.5), \cite{Loos_2020d} where one simply averages the SF-ADC(2)-s and SF-ADC(3) energies, is also tested in the following.
 
-We have also carried out spin-flip calculations within the TD-DFT framework (SF-TD-DFT), \cite{Shao_2003} with the same Q-CHEM 5.2.1 code. \cite{qchem}
+We have also carried out spin-flip calculations within the TD-DFT framework (SF-TD-DFT). \cite{Shao_2003}
 The B3LYP, \cite{Becke_1988b,Lee_1988a,Becke_1993b} PBE0 \cite{Adamo_1999a,Ernzerhof_1999} and BH\&HLYP global hybrid GGA functionals are considered, which contain 20\%, 25\%, 50\% of exact exchange, respectively. 
 These calculations are labeled as SF-TD-B3LYP, SF-TD-PBE0, and SF-TD-BH\&HLYP in the following.
 Additionally, we have also computed SF-TD-DFT excitation energies using range-separated hybrid (RSH) functionals: CAM-B3LYP (19\% of short-range exact exchange and 65\% at long range), \cite{Yanai_2004a} LC-$\omega$PBE08 (0\% of short-range exact exchange and 100\% at long range), \cite{Weintraub_2009a} and $\omega$B97X-V (16.7\% of short-range exact exchange and 100\% at long range). \cite{Mardirossian_2014}
 Finally, the hybrid meta-GGA functional M06-2X (54\% of exact exchange) \cite{Zhao_2008} and the RSH meta-GGA functional M11 (42.8\% of short-range exact exchange and 100\% at long range) \cite{Peverati_2011} are also employed.
 Note that all SF-TD-DFT calculations are done within the Tamm-Dancoff approximation. \cite{Hirata_1999}
 
-\alert{There also exist spin-flip extensions of EOM-CC methods, \cite{Krylov_2001a,Levchenko_2004,Manohar_2008,Casanova_2009a,Dutta_2013} and we consider here the spin-flip version of EOM-CCSD, named SF-EOM-CCSD. \cite{Krylov_2001a}
+There also exist spin-flip extensions of EOM-CC methods, \cite{Krylov_2001a,Levchenko_2004,Manohar_2008,Casanova_2009a,Dutta_2013} and we consider here the spin-flip version of EOM-CCSD, named SF-EOM-CCSD. \cite{Krylov_2001a}
 Additionally, Manohar and Krylov introduced a non-iterative triples correction to EOM-CCSD and extended it to the spin-flip variant. \cite{Manohar_2008} 
-Two types of triples corrections were proposed: (i) EOM-CCSD(dT) that uses the diagonal elements of the similarity-transformed CCSD Hamiltonian, and (ii) EOM-CCSD(fT) where the Hartree-Fock orbital energies are considered instead.}
+Two types of triples corrections were proposed: (i) EOM-CCSD(dT) that uses the diagonal elements of the similarity-transformed CCSD Hamiltonian, and (ii) EOM-CCSD(fT) where the Hartree-Fock orbital energies are considered instead.
+
+All spin-flip calculations have been performed with Q-CHEM 5.4.1. \cite{qchem}
+Note that symmetry labels may vary as different packages use different standard orientations. 
+Here, we have consistently followed the so-called Mulliken conventions. \cite{Mulliken_1955}
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -355,7 +356,7 @@ 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$  \\
-\alert{SF-EOM-CCSD} & \alert{$5.86$} & \alert{$6.27$} & \alert{$7.40$} \\[0.1cm]
+SF-EOM-CCSD & $5.86$ & $6.27$ & $7.40$ \\[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$ \\
@@ -408,11 +409,11 @@ We note that SF-ADC(2)-x [which scales as $\order*{N^6}$] is probably not worth
 This behavior was previously reported by Dreuw's group. \cite{Wormit_2014,Harbach_2014,Dreuw_2015}
 Overall, even with the best exchange-correlation functional, SF-TD-DFT is clearly outperformed by the more expensive SF-ADC models.
 
-\alert{We observe that SF-EOM-CCSD/aug-cc-pVTZ tends to underestimate by about \SI{1.5}{\kcalmol} the energy barrier compared to the TBE, an observation in agreement with previous results by Manohar and Krylov. \cite{Manohar_2008}
+We observe that SF-EOM-CCSD/aug-cc-pVTZ tends to underestimate by about \SI{1.5}{\kcalmol} the energy barrier compared to the TBE, an observation in agreement with previous results by Manohar and Krylov. \cite{Manohar_2008}
 This can be alleviated by including the triples correction with SF-EOM-CCSD(fT) and SF-EOM-CCSD(dT) (see {\SupInf} where we have reported the data from Ref.~\onlinecite{Manohar_2008}). 
 We also note that the SF-EOM-CCSD values for the energy barrier are close to the ones obtained with the more expensive (standard) CC3 method, yet less accurate than values computed with the cheaper SF-ADC(2)-s formalism. 
 Note that, in contrast to a previous statement, \cite{Manohar_2008} the (fT) correction performs better than the (dT) correction for the energy barrier. 
-However, for the excited states, the situation is reversed (see below).}
+However, for the excited states, the situation is reversed (see below).
 
 Concerning the multi-reference approaches with the minimal (4e,4o) active space, the TBEs are bracketed by the CASPT2 and NEVPT2 values that differ by approximately \SI{1.5}{\kcalmol} for all bases.
 In this case, the NEVPT2 values are fairly accurate with differences below half a \si{\kcalmol} compared to the TBEs.
@@ -496,9 +497,9 @@ SF-ADC(2.5) & 6-31+G(d) & $1.496$ & $3.328$ & $4.219$ \\
 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$ \\[0.1cm]
- \alert{SF-EOM-CCSD} & \alert{6-31+G(d)} & \alert{$1.663$} & \alert{$3.515$} & \alert{$4.275$} \\
- & \alert{aug-cc-pVDZ} & \alert{$1.611$} & \alert{$3.315$} & \alert{$4.216$} \\
- & \alert{aug-cc-pVTZ} & \alert{$1.609$} & \alert{$3.293$} & \alert{$4.245$} \\[0.1cm]
+ SF-EOM-CCSD & 6-31+G(d) & $1.663$ & $3.515$ & $4.275$ \\
+ & aug-cc-pVDZ & $1.611$ & $3.315$ & $4.216$ \\
+ & aug-cc-pVTZ & $1.609$ & $3.293$ & $4.245$ \\[0.1cm]
 		\end{tabular}
 	\end{ruledtabular}
 \end{table}
@@ -567,7 +568,7 @@ CC4 &6-31+G(d)& & $3.343$ & $4.067$ \\
  & aug-cc-pVTZ & & $[3.128]$\fnm[1] & $[4.032]$\fnm[2]\\[0.1cm]
 CCSDTQ &6-31+G(d)&  $1.464$  & $3.340$ & $4.073$ \\
 & aug-cc-pVDZ & $[1.433]$\fnm[3]& $[3.161]$\fnm[4]&  $[4.046]$\fnm[4] \\
-& aug-cc-pVTZ & \alert{$[\bf 1.433]$\fnm[5]} & $[\bf 3.125]$\fnm[6]&  $[\bf 4.038]$\fnm[6]\\[0.1cm]
+& aug-cc-pVTZ & $[\bf 1.433]$\fnm[5] & $[\bf 3.125]$\fnm[6]&  $[\bf 4.038]$\fnm[6]\\[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$ \\
@@ -615,19 +616,19 @@ This further motivates the ``pyramidal'' extrapolation scheme that we have emplo
 Again, the extended version, SF-ADC(2)-x, does not seem to be relevant in the present context with much larger errors than the other schemes.
 Also, as reported previously, \cite{Loos_2020d} SF-ADC(2)-s and SF-ADC(3) have mirror error patterns making SF-ADC(2.5) particularly accurate except for the doubly-excited state {\twoAg} where the error with respect to the TBE (\SI{0.140}{\eV}) is larger than the SF-ADC(2)-s error (\SI{0.093}{\eV}). 
 
-\alert{Interestingly, we observe that the SF-EOM-CCSD excitation energies are systematically larger than the TBEs by approximately \SI{0.2}{\eV} with a nice consistency throughout the various (singly- and doubly-) excited states. 
+Interestingly, we observe that the SF-EOM-CCSD excitation energies are systematically larger than the TBEs by approximately \SI{0.2}{\eV} with a nice consistency throughout the various (singly- and doubly-) excited states. 
 Moreover, SF-EOM-CCSD excitation energies are somehow closer to their SF-ADC(2)-s analogs (with an energy difference of about \SI{0.1}{\eV}) than the other schemes as already noticed by LeFrançois and co-workers. \cite{Lefrancois_2015}
 We see that the SF-EOM-CCSD excitations energies for the triplet state are larger of about \SI{0.3}{\eV} compared to the CCSD ones, which was also pointed out in the study of Manohar and Krylov. \cite{Manohar_2008}
 Again, our SF-EOM-CCSD results are very similar to the ones obtained in previous studies \cite{Manohar_2008,Lefrancois_2015}. 
 We can logically expect similar trend for SF-EOM-CCSD(fT) and SF-EOM-CCSD(dT) that lower the excitation energies and tend to be in better agreement with respect to the TBE (see {\SupInf}). 
-Note that the (dT) correction slightly outperforms the (fT) correction as previously observed \cite{Manohar_2008} and theoretically expected.}
+Note that the (dT) correction slightly outperforms the (fT) correction as previously observed \cite{Manohar_2008} and theoretically expected.
 
 Let us now move to the discussion of the results obtained with standard wave function methods that are reported in Table \ref{tab:D2h}.
-Regarding the \alert{multi-reference} calculations, the most striking result is the poor description of the {\sBoneg} ionic state, especially with the (4e,4o) active space where CASSCF predicts this state higher in energy than the {\twoAg} state. 
+Regarding the multi-reference calculations, the most striking result is the poor description of the {\sBoneg} ionic state, especially with the (4e,4o) active space where CASSCF predicts this state higher in energy than the {\twoAg} state. 
 Of course, the PT2 correction is able to correct the state ordering problem but cannot provide quantitative excitation energies due to the poor zeroth-order treatment.
 Another ripple effect of the unreliability of the reference wave function is the large difference between CASPT2 and NEVPT2 that differ by half an \si{\eV}. 
 This feature is characteristic of the inadequacy of the active space to model such a state.
-\alert{Additional MRCI and MRCI+Q calculations (reported in the {\SupInf}) confirm this.}
+Additional MRCI and MRCI+Q calculations (reported in the {\SupInf}) confirm this.
 For the two other states, {\tBoneg} and {\twoAg}, the errors at the CASPT2(4,4) and NEVPT2(4,4) levels are much smaller (below \SI{0.1}{\eV}).
 Using a larger active space resolves most of these issues: CASSCF predicts the correct state ordering (though the ionic state is still badly described in term of energetics), CASPT2 and NEVPT2 excitation energies are much closer, and their accuracy is often improved (especially for the triplet and doubly-excited states) although it is difficult to reach chemical accuracy (\ie, an error below \SI{0.043}{\eV}) on a systematic basis.
 
@@ -701,9 +702,9 @@ SF-ADC(2.5) & 6-31+G(d) & $0.234$ & $1.705$ & $2.087$ \\
 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$ \\[0.1cm]
-\alert{SF-EOM-CCSD} & \alert{6-31+G(d)} & \alert{$0.446$} & \alert{$1.875$} & \alert{$2.326$} \\
-& \alert{aug-cc-pVDZ} & \alert{$0.375$} & \alert{$1.776$} & \alert{$2.102$} \\
-& \alert{aug-cc-pVTZ}& \alert{$0.354$} & \alert{$1.768$} & \alert{$2.060$} \\
+{SF-EOM-CCSD} & {6-31+G(d)} & {$0.446$} & {$1.875$} & {$2.326$} \\
+& {aug-cc-pVDZ} & {$0.375$} & {$1.776$} & {$2.102$} \\
+& {aug-cc-pVTZ}& {$0.354$} & {$1.768$} & {$2.060$} \\
 		\end{tabular}
 	\end{ruledtabular}
 	
@@ -819,7 +820,7 @@ Concerning the singlet-triplet gap, each scheme predicts it to be positive.
 Although it provides a decent singlet-triplet gap value, SF-ADC(2)-x seems to particularly struggle with the singlet excited states ({\Aoneg} and {\Btwog}), especially for the doubly-excited state {\Aoneg} where it underestimates the vertical excitation energy by \SI{0.4}{\eV}. 
 Again, averaging the SF-ADC(2)-s and SF-ADC(3) transition energies is beneficial in most cases at the exception of {\Aoneg}.
 Although the basis set effects are larger than at the SF-TD-DFT level, they remain quite moderate at the SF-ADC level, and this holds for wave function methods in general.
-\alert{Concerning the SF-EOM-CCSD excitation energies at the {\Dfour} square planar equilibrium geometry, very similar conclusions to the ones provided in the previous section dealing with the excitation energies at the {\Dtwo} rectangular equilibrium geometry can be drawn: (i) SF-EOM-CCSD systematically and consistently overestimates the TBEs by approximately \SI{0.2}{\eV} and is less accurate than SF-ADC(2)-s, (ii) the non-iterative triples corrections tend to give a better agreement with respect to the TBE (see {\SupInf}), and (iii) the (dT) correction performs better than the (fT) one.}
+Concerning the SF-EOM-CCSD excitation energies at the {\Dfour} square planar equilibrium geometry, very similar conclusions to the ones provided in the previous section dealing with the excitation energies at the {\Dtwo} rectangular equilibrium geometry can be drawn: (i) SF-EOM-CCSD systematically and consistently overestimates the TBEs by approximately \SI{0.2}{\eV} and is less accurate than SF-ADC(2)-s, (ii) the non-iterative triples corrections tend to give a better agreement with respect to the TBE (see {\SupInf}), and (iii) the (dT) correction performs better than the (fT) one.
 
 Let us turn to the multi-reference results (Table \ref{tab:D4h}).
 For both active spaces, expectedly, CASSCF does not provide a quantitive energetic description, although it is worth mentioning that the right state ordering is preserved.
@@ -827,7 +828,7 @@ This is, of course, magnified with the (4e,4o) active space for which the second
 In particular SC-NEVPT2(4,4)/aug-cc-pVTZ and PC-NEVPT2(4,4)/aug-cc-pVTZ underestimate the singlet-triplet gap by \SI{0.072}{} and \SI{0.097}{\eV} and, more importantly, flip the ordering of {\Aoneg} and {\Btwog}.
 Although {\Aoneg} is not badly described, the excitation energy of the ionic state {\Btwog} is off by almost \SI{1}{\eV}.
 Thanks to the IPEA shift in CASPT2(4,4), the singlet-triplet gap is accurate and the state ordering remains correct but the ionic state is still far from being well described.
-The (12e,12o) active space significantly alleviates these effects, and, as usual now, the agreement between CASPT2 and NEVPT2 is very much improved for each state, though the accuracy of \alert{multi-reference} approaches remains questionable for the ionic state with, \emph{e.g.,} an error up to \SI{-0.093}{\eV} at the PC-NEVPT2(12,12)/aug-cc-pVTZ level.
+The (12e,12o) active space significantly alleviates these effects, and, as usual now, the agreement between CASPT2 and NEVPT2 is very much improved for each state, though the accuracy of multi-reference approaches remains questionable for the ionic state with, \emph{e.g.,} an error up to \SI{-0.093}{\eV} at the PC-NEVPT2(12,12)/aug-cc-pVTZ level.
 
 Finally, let us analyze the excitation energies computed with various CC models that are gathered in Table \ref{tab:D4h}.
 As mentioned in Sec.~\ref{sec:CC}, we remind the reader that these calculations are performed by considering the {\Aoneg} state as reference, and that, therefore, 
@@ -855,14 +856,14 @@ This has been shown to be clearly beneficial for the automerization barrier and
 \item At the SF-ADC level, we have found that, as expected, the extended scheme, SF-ADC(2)-x, systematically worsen the results compared to the cheaper standard version, SF-ADC(2)-s. 
 Moreover, as previously reported, SF-ADC(2)-s and SF-ADC(3) have opposite error patterns which means that SF-ADC(2.5) emerges as an excellent compromise.
 
-\item \alert{SF-EOM-CCSD shows similar performance as the cheaper SF-ADC(2)-s formalism, especially for the excitation energies. 
-As previously reported, the two variants including non-iterative triples corrections, SF-EOM-CCSD(dT) and SF-EOM-CCSD(fT), improve the results, the (dT) correction performing slightly better for the vertical excitation energies computed at the {\Dtwo} and {\Dfour} equilibrium geometries.}
+\item SF-EOM-CCSD shows similar performance as the cheaper SF-ADC(2)-s formalism, especially for the excitation energies. 
+As previously reported, the two variants including non-iterative triples corrections, SF-EOM-CCSD(dT) and SF-EOM-CCSD(fT), improve the results, the (dT) correction performing slightly better for the vertical excitation energies computed at the {\Dtwo} and {\Dfour} equilibrium geometries.
 
 \item For the {\Dfour} square planar structure, a faithful energetic description of the excited states is harder to reach at the SF-TD-DFT level because of the strong multi-configurational character.
 In such scenario, the SF-TD-DFT excitation energies can exhibit errors of the order of \SI{1}{\eV} compared to the TBEs. 
 However, it was satisfying to see that the spin-flip version of ADC can lower these errors to \SIrange{0.1}{0.2}{\eV}.
 
-\item Concerning the \alert{multi-reference} methods, we have found that while NEVPT2 and CASPT2 can provide different excitation energies for the small (4e,4o) active space, the results become highly similar when the larger (12e,12o) active space is considered. 
+\item Concerning the multi-reference methods, we have found that while NEVPT2 and CASPT2 can provide different excitation energies for the small (4e,4o) active space, the results become highly similar when the larger (12e,12o) active space is considered. 
 From a more general perspective, a significant difference between NEVPT2 and CASPT2 is usually not a good omen and can be seen as a clear warning sign that the active space is too small or poorly chosen.
 The ionic states remain a struggle for both CASPT2 and NEVPT2, even with the (12e,12o) active space.
 
@@ -886,4 +887,15 @@ Included in the {\SupInf} are the raw data, additional calculations and geometri
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \bibliography{CBD}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+
+\begin{center}
+	\textbf{Table of Contents Graphic}
+	\\
+	\bigskip
+	\boxed{\includegraphics[width=\linewidth]{TOC}}
+\end{center}
+\bigskip
+
 \end{document}

Response_Letter/v2/Response_Letter.tex

@@ -0,0 +1,26 @@
+\documentclass[10pt]{letter}
+\usepackage{UPS_letterhead,xcolor,mhchem,ragged2e,hyperref}
+\newcommand{\alert}[1]{\textcolor{red}{#1}}
+\definecolor{darkgreen}{HTML}{009900}
+
+
+\begin{document}
+
+\begin{letter}%
+{To the Editors of the Journal of Chemical Physics,}
+
+\opening{Dear Editors,}
+
+\justifying
+Please find attached a second revised version of the manuscript entitled 
+\begin{quote}
+	\textit{``Reference Energies for Cyclobutadiene: Automerization and Excited States''}.
+\end{quote}
+We have taken two of the three optional minor suggestions of the reviewer into account as we prefer to stick with the more systematic EOM-SF-CCSD notations (point 2.).
+In addition, we have taken into account the non-scientific changes requested by the editorial team.
+We look forward to hearing from you.
+
+\closing{Sincerely, the authors.}
+
+\end{letter}
+\end{document}

Response_Letter/v2/UPS_letterhead.sty

@@ -0,0 +1,70 @@
+%ANU etterhead Yves
+%version 1.0	12/06/08
+%need to be improved
+
+
+\RequirePackage{graphicx}
+
+%%%%%%%%%%%%%%%%%%%%% DEFINE USER-SPECIFIC MACROS BELOW %%%%%%%%%%%%%%%%%%%%%
+\def\Who     {Pierre-Fran\c{c}ois Loos}
+\def\What    {Dr}
+\def\Where   {Universit\'e Paul Sabatier}
+\def\Address {Laboratoire de Chimie et Physique Quantiques}
+\def\CityZip {Toulouse, France}
+\def\Email   {loos@irsamc.ups-tlse.fr}
+\def\TEL     {+33 5 61 55 73 39}
+\def\URL     {} % NOTE: use $\sim$ for tilde
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MARGINS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\textwidth       6in
+\textheight      9.25in
+\oddsidemargin   0.25in
+\evensidemargin  0.25in
+\topmargin       -1.50in
+\longindentation 0.50\textwidth
+\parindent       5ex
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%% ADDRESS MACRO BELOW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\address{
+    \includegraphics[height=0.7in]{CNRS_logo.pdf}     \hspace*{\fill}\includegraphics[height=0.7in]{UPS_logo.pdf}
+	\\
+	\hrulefill
+	\\
+	{\small \What~\Who\hspace*{\fill} Telephone:\ \TEL
+	\\
+	\Where\hspace*{\fill}  Email:\ \Email
+	\\
+	\Address\hspace*{\fill}
+	\\
+	\CityZip\hspace*{\fill} \URL}
+	 }
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% OTHER MACROS BELOW %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%\signature{\What~\Who}
+
+\def\opening#1{\ifx\@empty\fromaddress
+ \thispagestyle{firstpage}
+ \hspace*{\longindendation}\today\par
+ \else \thispagestyle{empty}
+ {\centering\fromaddress \vspace{5\parskip} \\
+\today\hspace*{\fill}\par}
+ \fi
+ \vspace{3\parskip}
+ {\raggedright \toname \\ \toaddress \par}\vspace{3\parskip}
+ \noindent #1\par\raggedright\parindent 5ex\par
+ }
+
+%I do not know what does the macro below
+
+%\long\def\closing#1{\par\nobreak\vspace{\parskip}
+ %\stopbreaks
+ %\noindent
+ %\ifx\@empty\fromaddress\else
+ %\hspace*{\longindentation}\fi
+ %\parbox{\indentedwidth}{\raggedright
+ %\ignorespaces #1\vskip .65in
+ %\ifx\@empty\fromsig
+ %\else \fromsig \fi\strut}
+ %\vspace*{\fill}
+% \par}