CBD/Response_Letter/Response_Letter.tex

\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 revised version of the manuscript entitled 
\begin{quote}
	\textit{``Reference Energies for Cyclobutadiene: Automerization and Excited States''}.
\end{quote}
We thank the reviewers for their constructive comments and to support publication of the present manuscript.
Our detailed responses to their comments can be found below.
For convenience, changes are highlighted in red in the revised version of the manuscript. 
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.}

\newpage

%%% REVIEWER 1 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#1}
 
{This article presents a survey of spin-flip TD-DFT, spin-flip ADC, multireference (CASSCF and MRPT), and equation-of-motion coupled cluster methods as applied to the automerization and vertical excitation energies of cyclobutadiene (CBD). As the smallest example of anti-aromaticity (and one of the smallest and most interesting exemplars of strong PJT distortion), CBD is an illuminating and challenging test case for these methods. (EOM-)CCSDTQ values, with a “pyramidal” basis set extrapolation scheme are used as the newly-proposed theoretical best estimates, and limited selected full CI (CIPSI) calculations confirm their excellent accuracy. The authors reach some interesting and useful conclusions concerning the tested methods.
}
\\
\alert{
Thank you for supporting publication of the present manuscript.
As detailed below, we have taken into account the comments and suggestions of the reviewers that we believe have overall improved the quality of the present paper.
We thank the reviewer for his/her careful analysis of our manuscript.}

{This work could be published as-is in JPC, but some suggestions for ways in which the manuscript could be improved follow:}
\\
\alert{
}

\begin{enumerate}

\item 
{The authors opt not to test SF-EOM-CC methods. 
A justification or rationalization would be helpful. 
Also, is it expected that these methods would improve on SF-ADC and/or EOM-CC?}
\\
\alert{The authors thanks the reviewer for this comment. 
Results for SF-EOM-CCSD, SF-EOM-CCSD(dT) and SF-EOM-CCSD(fT) have been added in the manuscript (and in the supporting information) and are discussed in the text.}

\item 
{The issue of reference symmetry frame is very important at the D4h geometry. 
The correlated calculation (and often the reference SCF calculation) are performed in a D2h subgroup, of which there are two distinct possibilities: one with the C2’ axes running through the carbon atoms and one with the C2’ axes bisecting them. 
It seems the former has been employed. 
The latter actually could potentially provide a faster convergence to the A1g state since it exhibits strong mixing between the two major determinants via T2 even at the CCSD level. 
However, this same property leads to a distinct inability to properly access the B1g ground state via a single excitation in EOM-CC. 
Some illuminating comments on this issue would be welcome.}
\\
\alert{We thank the reviewer for this interesting comment. 
Indeed, at the $D_{4h}$ $T_1$ optimized geometry, we have used the conventional standard orientation where two $C_2$ axes run through the carbon atoms. 
In this conventional orientation, the singlet ground state $1 ^1B_{1g}$ remains $1 ^1B_{1g}$ in the $D_{2h}$ point group and the singlet excited state $1 ^1A_{1g}$ becomes $1 ^1A_g$ in the $D_{2h}$ point group. 
As pointed out by the reviewer, upon rotating the molecular framework by 45 degrees in the ($xy$) plane, the two $C_2$ axes then bisect the carbon-carbon bonds. 
This induces a different orbital picture. 
The correlation between the orbitals and states in the new molecular framework are illustrated in the figure below at the CASSCF(4,4) level. 
In this new orientation, the two singlet states $1 ^1B_{1g}$ and $1 ^1A_{1g}$ become both $1 ^1A_{g}$ in the $D_{2h}$ point group. 
Because of the different orbital picture (the frontier orbitals are localized on two carbon atoms in the standard orientation and on four carbon atoms in the other orientation), the new CI coefficients resulting from this rotation bring also a different wave function representation. 
Whereas the $1 ^1B_{1g}$ ground state is described in a one-electron-excitation picture in the standard orientation (the $1 ^1A_{1g}$ excited state involves a double excitation), the corresponding $1 ^1B_{1g}$ ground state in the new orientation involves a two-electron-excitation picture (the $1 ^1A_{1g}$ excited state also involves a double excitation). 
Of course, these two representations are perfectly equivalent at the CASSCF level which describes single and double excitations on the same footing. 
This is obviously not the case in linear response theory, as pointed out by the reviewer.
As mentioned in our manuscript in Section II.B., for the $D_{4h}$ arrangement, we have considered the lowest closed-shell singlet state $^1A_{g}$ as reference. 
Because this state has a substantial double-excitation character, we expect a significant error at the CCSD level. 
The $1 ^1B_{1g}$ ground state is obtained as a singly excited state from that reference, while the $1 ^1B_{2g}$ excited state should also be a mixture involving a double excitation. 
In the other (non-standard) orientation, the lowest $^1A_g$ state correlates with the $1 ^1B_{1g}$ ground state, which in this orientation has a strong double-excitation character. 
Then, the $1 ^1 A_{1g}$ excited state has also a strong double-excitation character, while the $1 ^1B_{2g}$ excited state is obtained by one-electron excitation. 
Thus, whatever the orientation of the molecule, we will face the same problem for the reference state. 
Note that in the case of the SF formalism, these three singlet states should all be described correctly if one takes the $1 ^3A_{2g}$ state as a reference high spin state, whatever the orientation.
The interested reader will find these comments in the supporting information of the revised manuscript.}

\includegraphics[width=\textwidth]{MOs}

\item{The authors note a significant improvement in the MRPT results as the active space is enlarged. 
However, it seems to me that the most appropriate active space (for the D4h geometry at least) is in fact (2e,2o) [i.e. $Eg^2$ at D4h]. 
Within this space, the CI coefficients become fixed at D4h, leading to an “exact” SCF reference, at least in terms of static correlation. 
Perhaps the major problem with the MRPT results is not active space insufficiency, then, but intruder states? 
Can the authors perform MRCI+Q  or MRAQCC calculations for comparison?}
\\
\alert{We agree with the reviewer that at the $D_{4h}$ geometry the (2e,2o) active space would be enough to describe the pure static correlation. 
However, to calculate the automerization barrier, we need to make the energy difference between the energy obtained for the ground state at the $D_{4h}$ geometry and that at the $D_{2h}$ geometry. 
At this last geometry, the correct description of the static correlation requires including (4e,4o) in the active space (i.e., all valence $\pi$ orbitals).
In addition, there are states with ionic character which required including the dynamic electron correlation (in particular the $\sigma$-$\pi$ polarization). 
Thus, the improvement of our results by including all $\sigma_{CC}$ is rather expected. 
We believe that the large differences observed between CASPT2 and NEVPT2 for the (4e,4o) active space is a consequence of the small active space.
As a matter of fact, when the active space is enlarged, all these issues disappear.
Note that we have minimized the intruder state problem by using an appropriate level shift and that this potential problem is not present at the NEVPT2 level.
As suggested by the reviewer, we have now added some results at the MRCI+Q level.
}

\item 
{It seems that extrapolated CCSDTQ/aQZ values are available for the automerization barrier. 
Why are the aTZ numbers used as the TBE instead?}
\\
\alert{For the sake of consistency with the excitation energies and comparison, we have defined all the TBEs of the manuscript at the aug-cc-pVTZ level.
We believe that aug-cc-pVTZ is an adequate basis in order to get accurate values for the automerization barrier and the vertical excitation energies.
Defining the TBE at the aug-cc-pVQZ level would make comparison with other methods quite expensive (and sometimes undoable for some of the most expensive methods.}

\end{enumerate}

%%% REVIEWER 2 %%%
\noindent \textbf{\large Authors' answer to Reviewer \#2}
 
{This is a useful addition to the literature, presenting extensive benchmarks on a popular system, cyclobutadiene or CB. I recommend it for the publication once the following issues are addressed.}
\\
\alert{Thank you for these positive comments and for supporting publication of our manuscript. 
Below, we address the points raised by Reviewer \#2.
}

\begin{enumerate}

\item 
{The results for EOM-SF-CCSD and EOM-SF-CCSD(fT/dT) must be included in the paper and in the
analysis/discussion of the results. Why to exclude the best-performing SF methods? Since this paper aspires to be a comprehensive benchmark on CB, I believe it is absolutely essential.  Moreover, some of these results are already available (e.g., Ref. 105 has the results for excitation energies obtained in the same basis -- aug-cc-pVTZ that is used in the paper). Even if one needs to redo the calculations, they are very quick and can be done on a laptop in a few minutes.}
\\
\alert{As mentioned in the response to Reviewer \#1, results for SF-EOM-CCSD, SF-EOM-CCSD(dT) and SF-EOM-CCSD(fT) have been added in the manuscript (and in the supporting information) and are discussed in the text.}

\item 
{The comparison of SF-ADC with EOM-SF-CCSD will be illuminating for the readers. 
For example, EOM-SF is more robust wrt reference spin-contamination compared to SF-ADC because of the CC ansatz. 
There could be other interesting differences to discuss.}
\\
\alert{See previous point.}

\item 
{I also recommend to include EOM-DEA-CCSD results -- this is another extension of EOM-CCSD, which can treat diradicals. 
It does not suffer from spin-contamination. 
The method is available in Q-Chem. 
See here for theory description and examples: J. Chem. Phys. 154, 114115 (2021). EOM-DIP is another method, which can deal wit this type of electronic structure, but it has difficulties with diffuse basis sets (e.g., J. Chem. Phys. 135, 084109 (2011)) -- so I am not asking to add the DIP numbers, but mentioning it would be appropriate.}
\\
\alert{EOM-DEA-CCSD results have been added to the supporting information of the revised manuscript.
EOM-DIP has also been mentioned as suggested by the reviewer.}

\item 
{The analysis would benefit greatly if the authors provide Head-Gordon's indices, which can be used to compare wave-functions computed by different methods in a meaningful way, as illustrated here:J. Chem. Theo. Comp. 14, 638 (2018). }
\\
\alert{The authors thank the referee for this valuable comment. 
Unfortunately, in order to obtain Head-Gordon's indices for the different spin-flip methods we would have to do all the calculations or at least for the aug-cc-pVTZ basis which will take too much time and resources. We have mentioned these indices (alongside relevant references) in the text and we will definitely use them in future works.}

\item 
{CAS-based methods are multi-reference (and also able to treat multi-configurational wfns). 
EOM-SF and EOM-EE are single-reference methods that are able to describe multi-configurational wfns. 
Please correct the section names and discussion appropriately.}
\\
\alert{We have modified the section names and discussion accordingly to the reviewer's suggestion.}

\item 
{Abstract and introduction:
please insert the word 'standard' before 'time-dependent density-functional theory (TD-DFT) or equation-of-motion ... are notoriously known to struggle in such situations.' 
The SF and DEA/DIP variants of these methods do not struggle.
}
\\
\alert{This has been corrected.}

\item
{Intro: "Of course, single-reference methods are naturally unable to describe such situations." 
This is incorrect -- see above (EOM-SF/DIP/DEA are single reference methods capable of describing multi-configurational wfns). 
Adding the word 'standard' might help. 
Below: "and remain tortuous for state-of-the-art methods ..." -- again, need to correct, e.g., consider 'remains challenging for standard hierarchy of EOM-CC methods that are using ground-state Hartree-Fock reference'.}
\\
\alert{We have performed these two corrections.}

\item{First page, last paragraph -- replace multi-configurational by multi-reference, as per above.}
\\
\alert{This has been corrected.}

\item{Last paragraph of intro -- here you can introduce the idea of single-reference approach to multi-reference wfns and describe SF/DEA/DIP methods.  
Please do not call SF 'cheaper' -- this does not make sense (the cost depends on correlation treatment). 
SF is more robust and more effective at each correlation level. 
It is also more black-box, as it does not require active-space selection. 
Also, EOM-SF can systematically converge to the exact FCI answer once sufficiently high excitations are included. 
'Obviously, spin-flip methods have their own flaws, especially spin' -- why 'obviously'? 
I suggest to reword, e.g., 'One drawback of SF methods is ...'}
\\
\alert{The reviewer is right. We have followed the suggestion of the reviewer and performed the required modifications.}

\item{Section IIC -- please rename.}
\\
\alert{Done.}

\end{enumerate}

\end{letter}
\end{document}