loos/CBD
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

corrections martial

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-04-11 13:33:38 +02:00
parent 45e72c53f7
commit 7258819f6f
1 changed files with 24 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
Manuscript/CBD.tex
View File

@ -86,9 +86,9 @@
\begin{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 adiabatic time-dependent density-functional theory (TD-DFT) or equation-of-motion coupled cluster (EOM-CC), are notoriously known to struggle in such situations.
In this work, using a large panel of methods and basis sets, we provide an extensive computational study of the automerization barrier (defined as the difference between the square and rectangular ground-state energies) and the vertical excitation energies of both the $D_{2h}$ and $D_{4h}$ equilibrium structures.
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 states with multi-configurational character, 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.
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.
\end{abstract}
@ -105,7 +105,7 @@ Indeed, each computational model has its own theoretical and/or technical limita
Speaking of difficult tasks, the cyclobutadiene (CBD) molecule has been a real challenge for both 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 making 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 in the rectangular geometry (see below).
This degeneracy is lifted by the so-called pseudo 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 in the rectangular geometry (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.
@ -135,7 +135,7 @@ For example, the \textit{Configuration Interaction using a Perturbative Selectio
%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) de-excitation and excitation 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 not only to the spin incompleteness in the spin-flip expansion but also to the spin contamination of the reference configuration. \cite{Casanova_2020}
Obviously, spin-flip methods have their own flaws, especially 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}
%both solutions being associated with an increased computational cost.
@ -194,16 +194,16 @@ With respect to {\Aoneg}, {\sBoneg} has a dominant double excitation character,
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Multi-configurational calculations}
\label{sec:Multi}
State-averaged CASSCF (SA-CASSCF) calculations are performed with MOLPRO. \cite{Werner_2020}
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).
Two active spaces have been considered: (i) a minimal (4e,4o) active space including the 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}
\titou{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.25}{\hartree} to avoid systematic overestimation of the vertical excitation energies. \cite{Ghigo_2004,Schapiro_2013,Zobel_2017,Sarkar_2022}
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}
%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}
All these calculations are 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -221,7 +221,7 @@ The spin-flip version of our recently proposed composite approach, namely SF-ADC
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}
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-BLYP, SF-TD-B3LYP, SF-TD-PBE0, and SF-TD-BH\&HLYP in the following.
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}
@ -284,7 +284,7 @@ Two different sets of geometries obtained with different levels of theory are co
First, because the automerization barrier is obtained as a difference of energies computed at distinct geometries, it is paramount to obtain these at the same level of theory.
However, due to the fact that the ground state of the square arrangement is a transition state of singlet open-shell nature, it is technically difficult to optimize the geometry with high-order CC methods.
Therefore, we rely on CASPT2(12,12)/aug-cc-pVTZ for both the {\Dtwo} and {\Dfour} ground-state structures.
\titou{(Note that these optimizations are done without IPEA shift.)}
(Note that these optimizations are done without IPEA shift but with a level shift and a state-specific reference CASSCF wave function.)
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}).
@ -347,8 +347,8 @@ SF-TD-BH\&HLYP & $11.90$ & $12.02$ & $12.72$ & $12.73$ \\
SF-TD-M06-2X & $9.32$ & $9.62$ & $10.35$ & $10.37$ \\
SF-TD-CAM-B3LYP& $18.05$ & $18.10$ & $18.83$ & $18.83$ \\
SF-TD-$\omega$B97X-V & $18.26$ & $18.24$ & $18.94$ & $18.92$ \\
SF-TD-LC-$\omega$PBE08 & $19.05$ & $18.98$ & $19.74$ & $19.71$ \\[0.1cm]
SF-TD-M11 & $11.03$ & $10.25$ & $11.22$ & $11.12$ \\
SF-TD-LC-$\omega$PBE08 & $19.05$ & $18.98$ & $19.74$ & $19.71$ \\
SF-TD-M11 & $11.03$ & $10.25$ & $11.22$ & $11.12$ \\[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$ & \\
@ -382,8 +382,8 @@ CCSDTQ & $7.51$ & $[ 7.89]$\fnm[4]& $[\bf 8.93]$\fnm[5]& $[ 9.11]$\fnm[6]\\
%%% FIGURE II %%%
\begin{figure*}
\includegraphics[width=0.8\linewidth]{AB_AVTZ}
\caption{Automerization barrier (in \si{\kcalmol}) of CBD at various levels of theory using the aug-cc-pVTZ basis.
See {\SupInf} for the total energies.}
\caption{\titou{Automerization barrier (in \si{\kcalmol}) of CBD at various levels of theory using the aug-cc-pVTZ basis.
See {\SupInf} for the total energies.}}
\label{fig:AB}
\end{figure*}
%%% %%% %%% %%%
@ -532,7 +532,7 @@ CASPT2(4,4) &6-31+G(d)& $1.440$ & $3.162$ & $4.115$ \\
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]
& aug-cc-pVQZ & $1.384$ & $2.408$ & $4.116$ \\[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$ \\
@ -626,12 +626,12 @@ Again, the extended version, SF-ADC(2)-x, does not seem to be relevant in the pr
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}).
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 multi-configurational calculations, the most striking result is the poor description of the {\sBoneg} ionic state, especially with the (4e,4o) active space where CASSCF predicts this state higher in energy than the {\twoAg} state.
Of course, the PT2 correction is able to correct the root-flipping problem but cannot provide quantitative excitation energies due to the poor zeroth-order treatment.
Regarding the multi-configurational calculations, the most striking result is the poor description of the {\sBoneg} \titou{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.
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.
Using a larger active space resolves most of these issues: CASSCF predicts the correct state ordering (though the \titou{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.
Finally, for the CC models (Table \ref{tab:D2h}), the two states with a large $\%T_1$ value, {\tBoneg} and {\sBoneg}, are already extremely accurate at the CC3 level, and systematically improved by CCSDT and CC4.
This trend is in line with the observations made on the QUEST database. \cite{Veril_2021}
@ -808,7 +808,7 @@ Unfortunately, due to technical limitations, we could not compute $\%T_1$ values
However, it is clear from the inspection of the wave function that, with respect to the {\sBoneg} ground state, {\Atwog} and {\Btwog} are dominated by single excitations, while {\Aoneg} has a strong double excitation character.
As for the previous geometry we start by discussing the SF-TD-DFT results (Table \ref{tab:sf_D4h}), and in particular the singlet-triplet gap, \ie, the energy difference between {\sBoneg} and {\Atwog}.
For all functionals, this gap is small (basically below \SI{0.1}{\eV} while the TBE value is \SI{0.144}{\eV}) but it is worth mentioning that B3LYP and PBE0 incorrectly deliver a negative singlet-triplet gap (hence a triplet ground state).
For all functionals, this gap is small (basically below \SI{0.1}{\eV} while the TBE value is \SI{0.144}{\eV}) but it is worth mentioning that B3LYP and PBE0 incorrectly deliver a negative singlet-triplet gap (hence a triplet ground state at this geometry).
Increasing the fraction of exact exchange in hybrids or relying on RSHs (even with a small amount of short-range exact exchange) allows to recover a positive gap and a singlet ground state.
At the SF-TD-DFT level, the energy gap between the two singlet excited states, {\Aoneg} and {\Btwog}, is particularly small and grows moderately with the amount of exact exchange at short range.
The influence of the exact exchange on the singlet energies is quite significant with an energy difference of the order of \SI{1}{\eV} between the functional with the smallest amount of exact exchange (B3LYP) and the functional with the largest amount (M06-2X).
@ -827,11 +827,11 @@ Although the basis set effects are larger than at the SF-TD-DFT level, they rema
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.
This is, of course, magnified with the (4e,4o) active space for which the second-order perturbative treatment is unable to provide a satisfying description due to the restricted active space.
This is, of course, magnified with the (4e,4o) active space for which the second-order perturbative treatment is unable to provide a satisfying description due to the limited active space.
In particular SC-NEVPT2(4,4)/aug-cc-pVTZ and PC-NEVPT2(4,4)/aug-cc-pVTZ underestimate the singlet-triplet gap by \SI{0.072}{} and \SI{0.097}{\eV} and, more importantly, flip the ordering of {\Aoneg} and {\Btwog}.
Although {\Aoneg} is not badly described, the excitation energy of the ionic state {\Btwog} is off by 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 multi-configurational 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.
Although {\Aoneg} is not badly described, the excitation energy of the \titou{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 \titou{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 multi-configurational approaches remains questionable for the \titou{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,
@ -865,7 +865,7 @@ However, it was satisfying to see that the spin-flip version of ADC can lower th
\item Concerning the multi-configurational 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.
The \titou{ionic} states remain a struggle for both CASPT2 and NEVPT2, even with the (12e,12o) active space.
\item In the context of CC methods, although the inclusion of triple excitations (via CC3 or CCSDT) yields very satisfactory results in most cases, the inclusion of quadruples excitation (via CC4 or CCSDTQ) is mandatory to reach high accuracy (especially in the case of doubly-excited states).
Finally, we point out that, considering the error bar related to the CIPSI extrapolation procedure, CCSDTQ and CIPSI yield equivalent excitation energies, hence confirming the outstanding accuracy of CCSDTQ in the context of molecular excited states.