seniority/Manuscript/seniority.tex

\documentclass[aip,jcp,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,wrapfig,txfonts}
\usepackage[version=4]{mhchem}
%\usepackage{natbib}
%\bibliographystyle{achemso}

%fk
\usepackage{caption}
\usepackage{subcaption}
\usepackage{comment}
\newcommand{\fk}[1]{\textcolor{blue}{#1}}

\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\alert}[1]{\textcolor{black}{#1}}
\usepackage[normalem]{ulem}
\newcommand{\titou}[1]{\textcolor{red}{#1}}
\newcommand{\trashPFL}[1]{\textcolor{red}{\sout{#1}}}
\newcommand{\PFL}[1]{\titou{(\underline{\bf PFL}: #1)}}
\newcommand{\toto}[1]{\textcolor{green}{#1}}
\newcommand{\trashAS}[1]{\textcolor{green}{\sout{#1}}}
\newcommand{\AS}[1]{\toto{(\underline{\bf AS}: #1)}}
\newcommand{\ant}[1]{\textcolor{orange}{#1}}
\newcommand{\SI}{\textcolor{blue}{Supporting Information}}

\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\QP}{\textsc{quantum package}}

\newcommand{\hI}{\Hat{1}}
\newcommand{\hT}{\Hat{T}}
\newcommand{\hH}{\Hat{H}}
\newcommand{\bH}{\Bar{H}}
\newcommand{\ERI}[2]{v_{#1}^{#2}}

\newcommand{\EHF}{E_\text{HF}}
\newcommand{\EDOCI}{E_\text{DOCI}}
\newcommand{\EFCI}{E_\text{FCI}}
\newcommand{\ECC}{E_\text{CC}}
\newcommand{\EVCC}{E_\text{VCC}}
\newcommand{\EpCCD}{E_\text{pCCD}}

\newcommand{\si}{\sigma}

\usepackage[
	colorlinks=true,
    citecolor=blue,
    breaklinks=true
	]{hyperref}
\urlstyle{same}

\begin{document}

\newcommand{\LCPQ}{Laboratoire de Chimie et Physique Quantiques (UMR 5626), Universit\'e de Toulouse, CNRS, UPS, France}

\title{Configuration interaction with seniority number and excitation degree}

\author{F\'abris Kossoski}
\email{fkossoski@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}
\author{Yann Damour}
\affiliation{\LCPQ}
\author{Pierre-Fran\c{c}ois Loos}
\email{loos@irsamc.ups-tlse.fr}
\affiliation{\LCPQ}

% Abstract
\begin{abstract}
Here comes the abstract.
%\bigskip
%\begin{center}
%        \boxed{\includegraphics[width=0.4\linewidth]{TOC}}
%\end{center}
%\bigskip
\end{abstract}

% Title
\maketitle


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In electronic structure theory, configuration interaction (CI) methods allow for a systematic way to obtain approximate or exact solutions of the electronic Hamiltonian,
by expanding the wave function as a linear combination of Slater determinants (or configuration state functions).
At the full CI (FCI) level, the complete Hilbert space is spanned in the wave function expansion, leading to the exact solution for a given one-particle basis set.
Except for very small systems, the FCI limit is unnatainable, and in practice the expansion of the CI wave function must be trunctated.
The question is then how to construct an effective and computationally tractable hierarchy of truncated CI methods
that best recover the correlation energy, understood as the energy difference between the FCI and the mean-field restricted Hartree-Fock (HF) solutions.
%that lead as fast as possible to the FCI limit.

The most well-known and popular class of CI methods is excitation-based,
where one accounts for all determinants generated by exciting up to $e$ electrons from a given close-shell reference, which is usually the restricted HF solution, but does not have to.
In this way, the excitation degree $e$ parameter defines the sequence
CI with single excitations (CIS), CI with single and double excitations (CISD), CI with single, double, and triple excitations (CISDT), and so on.
Excitation-based CI manages to quickly recover weak (dynamic) correlation effects, but struggles at strong (static) correlation regimes.
Importantly, the number of determinants $N_{det}$ (which control the computational cost) scale polynomially with the number of electrons $N$ as $N^{2d}$.
%This means that the contribution of higher excitations become progressively smaller.

%In turn, seniority-based CI is specially targeted to describe static correlation.
%\fk{Still have to work in this paragraph.}
Alternatively, CI methods based on the seniority number \cite{Ring_1980} have been proposed \cite{Bytautas_2011}.
In short, the seniority number $s$ is the number of unpaired electrons in a given determinant.
The seniority zero ($s = 0$) sector has been shown to be the most important for static correlation, while higher sectors tend to contribute progressively less ~\cite{Bytautas_2011,Bytautas_2015,Alcoba_2014b,Alcoba_2014}.
% scaling
However, already at the sCI0 level the number of determinants scale exponentially with $N$, since excitations of all excitation degrees $e$ are included.
Therefore, despite the encouraging successes of seniority-based CI methods, their unfourable computational scaling restricts applications to very small systems.
Besides CI, other methods that exploit the concpet of seniority number have been pursued. \cite{Henderson_2014,Chen_2015,Bytautas_2018}
% Seniority Number in Valence Bond Theory
%https://doi.org/10.1021/acs.jctc.5b00416
% Seniority based energy renormalization group (Ω-ERG) approach in quantum chemistry: Initial formulation and application to potential energy surfaces
%https://doi.org/10.1016/j.comptc.2018.08.011

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hierarchy configuration interaction}
\label{sec:hCI}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

At this point, we notice the current dicothomy.
When targeting static correlation, seniority-based CI methods tend to have a better performance than excitation-based CI, despite the higher computational cost.
The latter class of methods, in contrast, are well-suited for recovering dynamic correlation, and only at polynomial cost with system size.
% tackling
Ideally, we aim for a method that captures most of both static and dynamic correlation, with as few determinants as possible.
With this goal in mind, we propose a new partitioning of the Hilbert space, named hierarchy configuration interaction (hCI).
It combines both the excitation degree $e$ and the seniority number $s$ into one single hierarchy parameter $h$,
\begin{equation}
	\label{eq:h}
	h = \frac{e+s/2}{2},
\end{equation}
which assumes half-integer values.
% open-shell
Fig.~\ref{fig:allCI} shows how the Hilbert space is populated in excitation-based CI, seniority-based CI, and our hybrid hCI methods.

%%% FIG 1 %%%
\begin{figure}[h!]
%\centering
\begin{subfigure}[b]{0.48\linewidth}
\includegraphics[width=\linewidth]{table_exc_full}
        \caption{Excitation-based CI.}
        \label{fig:exc3}
\end{subfigure}
%    \hfill
\begin{subfigure}[b]{0.48\linewidth}
\includegraphics[width=\linewidth]{table_sen_full}
        \caption{Seniority-based CI.}
        \label{fig:sen}
\end{subfigure}
%    \hfill
\begin{subfigure}[b]{0.48\linewidth}
\includegraphics[width=\linewidth]{table_hCI}
        \caption{Hierarchy-based CI.}
        \label{fig:hCI}
\end{subfigure}
	\caption{Partionining of the full Hilbert space into blocks of specific excitation degree $e$ (with respect to a closed-shell determinant) and seniority number $s$.
	Each of three classes of CI methods truncate this $e$-$s$ map differently, and each color tone represents the added determinants at a given CI level.}
        \label{fig:allCI}
\end{figure}
%%% %%% %%%

We have three key justifications for this new CI hierarchy.
The first one is physical.
We know that low degree excitations and low seniority sectors, when looked at individually, often have the most important contribution to the FCI expansion.
%carry the most important weights.
By combining $e$ and $s$ as is eq.~\ref{eq:h}, we ensure that both directions in the excitation-seniority map (see Fig.~\ref{fig:allCI}) will be contemplated.
Rather than filling the map top-bottom (as in excitation-based CI) or left-right (as in seniority-based CI), the hCI methods fills it diagonally.
In this sense, we hope to recover dynamic correlation by moving right in the map (increasing the excitation degree while keeping a low seniority number),
at the same time as static correlation, by moving down (increasing the seniority number while keeping a low excitation degree).
%dynamic correlation is recovered with traditional CI.

The second justification is computational.
%Fig.~\ref{fig:scaling} also illustrates how the number of determinants within each block scales with the number of occupied orbitals $O$ and the number of virtual orbitals $V$.
In the hCI class of methods, each next level of theory accomodates additional determinants from different excitation-seniority sectors (each block of Fig.~\ref{fig:allCI}).
The key realization behind hCI is that the number of additional determinants presents the same scaling with respect to $N$, for all excitation-seniority sectors entering at a given hierarchy $h$.
%to $O$ and $V$, for all excitation-seniority sectors of a given hierarchy $h$.
%This computational realization represents the second justification for the introduction of the hCI method.
This further justifies the parameter $h$ as being the simple average between $e$ and $s/2$.

Each level of excitation-based CI has a hCI counterpart with the same scaling of $N_{det}$ with respect to $N$.
%However, hCI counts with additional half-integer levels of theory, with no parallel in excitation-based CI.
For example, in both hCI2 and CISD we have $N_{det} \sim N^4$, whereas in hCI3 and CISDT, $N_{det} \sim N^6$, and so on.
%the number of determinants of hCI2 and CISD scale as $O^2V^2$, those of hCI3 and CISDT scale as $O^3V^3$, and so on.
From this computational perspective, hCI can be seen as a more natural choice than the traditional excitation-based CI,
because if one can afford for, say, the $N_{det} \sim N^6$ cost of a CISDT calculation, than one can probably afford a hCI3 calculation, which has the same computational scaling.
Of course, in practice an integer-$h$ hCI method will have more determinants than its excitation-based counterpart (though the same scaling),
and thus one should first ensure whether including the lower-triangular blocks (going from CISDT to hCI3 in our example)
is a better strategy than adding the next column (going from CISDT to CISDTQ).
Therefore, here we decided to discuss the results in terms of the number of determinants, rather than the computational scaling,
which could make the comparison somewhat biased toward hCI.
%
The lowest level in hCI (hCI1) parallels CIS of excitation-based CI. 
However, the single excitations do not connect with the reference, at least for HF orbitals, and therefore CIS provides the same energy as HF.
In contrast, the paired doubles excitations of hCI1 do connect with the reference (as well as the singles, indireclty via the doubles).
Therefore, while the HF-based lowest level of excitation-based CI (CIS) does not improve with respect to the mean-field HF wave function,
the hCI1 counterpart already represents a minimally correlated model, with the very favourable $N_{det} \sim N^2$ scaling.
%number of determinants scaling only as $OV$.
%
In addition, hCI allows for half-integer values of $h$, with no parallel in excitation-based CI.
This gives extra flexibility in terms of choice of method.
%when evaluating the computational cost and desired accuracy of a calculation.
For a particular application with excitation-based CI, CISD might be too inaccurate, for example, while the price for the improved accuracy of CISDT might be too high.
hCI2.5 could represent an alternative, being more accurate than hCI2 and less expensive than hCI3.

Finally, the third justification for our hCI method is empirical and closely related to the computational motivation.
There are many possible ways to populate the Hilbert space starting from the a given reference determinant, 
and one can in principle formulate any systematic recipe that includes progressively more determinants.
Besides a physical or computational perspective, the question of what makes for a good recipe can be framed empirically.
Does our hCI class of methods perform better than excitation-based or seniority-based CI,
in the sense of recovering most of the correlation energy with the least computational effort?

A hybrid approach based on both excitation degree and seniority number has been proposed. \cite{Alcoba_2014,Raemdonck_2015,Alcoba_2018}
In these works, the authors established separate maximum values for the excitation and the seniority,
and either the union or the intersection between the two sets of determinants have been considered.
For the union case, the number of determinants grows exponentially with $N$,
while in the intersection approach the Hilbert space is filled rectangle-wise in our excitation-seniority map.
In the latter case, the scaling of $N_{det}$ would be dominated by the rightmost bottom block.
Bytautas et al.\cite{Bytautas_2015} explored a different hybrid scheme combining determinants from a complete active space and with a maximum seniority number.
In comparison to previous approaches, our hybrid hCI scheme has two key advantages.
First, it is defined by a single parameter that unifies excitation degree and seniority number (eq.\ref{eq:h}).
And second, each next level includes all classes of determinants sharing the same scaling with system size, as discussed before, thus keeping the method at a polynomial scaling.

Our main goal here is to assess the performance of hCI against excitation-based and seniority-based CI.
To do so, we evaluated how fast different observables converge to the FCI limit as a function of the number of determinants.
We have calculated the potential energy curves (PECs) for a total of 6 systems, 
\ce{HF}, \ce{F2}, \ce{N2}, 
%\ce{Be2}, \ce{H2O}, 
ethylene, \ce{H4}, and \ce{H8},
which display a variable number of bond breaking.
For the latter two, we considered linearly arranged and equally spaced hydrogen atoms, and computed PECs along the symmetric dissociation coordinates.
%For \ce{H2O}, we considered the symmetric stretching of the O$-$H bonds, 
For ethylene, we considered the C$=$C double bond stretching, while freezing the remaining internal coordinates.
%in both cases freezing the remaining internal coordinates.
Due to the (multiple) bond breaking, these are challenging systems for electronic structure methods,
and are often considered when assessing novel methodologies.
%
From the PECs, we evaluated the convergence of four observables: the nonparalellity error (NPE), the distance error, the equilibrium geometries, and the vibrational frequencies.
The NPE is defined as the maximum minus the minimum differences between the PECs obtained at given CI level and the exact FCI result.
We define the distance error as the the maximum plus and the minimum differences between a given PEC and the FCI result.
%This is an important metric because it captures the resemblance between the shape of the two PECs, 
%which in turn determine the relevant physical observables, as equilibrium geometries, vibrational frequencies, and dissociation energies.
Thus, while the NPE probes the similarity regarding the shape of the PECs, the distance error provides a measure of how their magnitudes compare.
From the PECs, we have also extrated the equilibrium geometries and vibrational frequencies (details can be found in the SI).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational details}
\label{sec:compdet}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The hCI method was implemented in {\QP} via a straightforward adaptation of the \textit{configuration interaction using a perturbative selection made iteratively} (CIPSI) algorithm \cite{Huron_1973,Giner_2013,Giner_2015},
by allowing only for determinants at a given hierarchy $h$.
The excitation-based CI, seniority-based CI, and FCI calculations presented here were also performed with the CIPSI algorithm implemented in {\QP}. \cite{Huron_1973,Giner_2013,Giner_2015,Garniron_2019}
In practice, we consider the CI energy to be converged when the second-order perturbation correction lies below $10^{-5}$ Hartree,
which requires considerably fewer determinants than the formal number of determinants (understood as all those that belong to a given CI level, regardless of their weight or symmetry).
Nevertheless, we decided to present the results as functions of the formal number of determinants, 
which are not related to the particular algorithmic choices of the CIPSI calculations.
%
All CI calculations were performed for the cc-pVDZ basis set and with frozen core orbitals.
For \ce{HF} we have also considered the cc-pVTZ and cc-pVQZ basis sets.

The CI calculations were performed with both canonical Hartree-Fock (HF) orbitals and optimized orbitals (oo).
In the latter case, the energy is obtained variationally in both the CI space and in the orbital parameter space, given rise to what may be called an orbital-optimized CI (oo-CI) method.
%We have also performed orbital optimized CI (oo-CI) calculations, where the energy is obtained variationally both in the CI space and in the orbital parameter space.
We employed the algorithm described elsewhere \cite{Damour_2021} and also implemented in {\QP} for optimizing the orbitals within a CI wave function.
In order to avoid converging to a saddle point solution, we employed a similar strategy as recently described in Ref. \cite{Hollett_2022}.
Namely, whenever the eigenvalue of the orbital rotation Hessian is negative and the corresponding gradient component $g_i$ lies below a given threshold $g_0$,
then this gradient component is replaced by $g_0 |g_i|/g_i$.
While we can never ensure that the obtained solutions are global minima in the orbital parameter space, we verified that in all cases surveyed here, the stationary solutions are real minima (rather than maxima or stationary points).

It is worth mentioning that obtaining smooth PECs for the orbital optimized calculations proved to be far from trivial.
First, the orbital optimization started from the HF orbitals of each geometry.
This usually lead to discontinuous PECs, meaning that distinct solutions of the orbital optimization have been found with our algorithm.
Then, at some geometry or geometries that seem to present the lowest lying solution,
the optimized orbitals were employed as the guess orbitals for the neighbouring geometries, and so on, until a new PEC is attained.
%orthonormalized
This protocol is repeated until the PEC built from the lowest lying orbital optimized solutions becomes continuous.
While we cannot guarantee that the presented solutions represent the global minima, we believe that in most cases the above protocol provides at least close enough solutions.
%Multiple solutions for the orbital optimization are usually found, meaning several local minimal in the orbital parameter landscape.
%meaning that the set of orbitals are stationary with respect to the energy.
We recall that saddle point solutions were purposedly avoided in our orbital optimization algorithm. If that was not the case, then even more stationary solutions would have been found.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
\label{sec:res}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\subsection{Correlation energies}

%\ce{F2}
%\ce{N2}
%\ce{HF}
%ethylene
%Linear \ce{H4} and \ce{H8}
%\ce{Be2}

%\subsection{Potential energy curves}

%\subsection{Nonparallelity errors and dissociation energies}
\subsection{Nonparallelity errors}

In Fig.~\ref{fig:plot_stat} we present the NPEs for the six systems studied, and for the three classes of CI methods,
as functions of the number of determinants.
%the potential energy curves (PECs) and the corresponding differences with respect to the FCI result, as well as the nonparalellity error (NPE) and the distance error.
The corresponding PECs and the energy differences with respect to the FCI results can be found in the SI.
The main result contained in Fig.~\ref{fig:plot_stat} concerns the overall faster convergence of the hCI methods when compared to excitation-based and seniority-based CI methods.
This is observed both for single bond breaking (\ce{HF} and \ce{F2}) as well as the more challenging double (ethylene), triple (\ce{N2}), and quadruple (\ce{H4}) bond breaking.
For (\ce{H8}), hCI and excitation-based CI perform similarly.
The convergence with respect to the number of determinants is slower in the latter, more challenging cases, irrespective of the class of CI methods, as would be expected.
But more importantly, the superiority of the hCI methods appears to be highlighted in the multiple bond break systems (compare ethylene and \ce{N2} with \ce{HF} and \ce{F2} in Fig.~\ref{fig:plot_stat}).

For \ce{HF} we also evaluated the convergence is affected by increasing the basis sets, going from cc-pVDZ to cc-pVTZ and cc-pVQZ basis sets (see Fig.Sx).
While a larger number of determinants is required to achieve the same level of convergence, as expected,
the convergence profiles remain very similar for all basis sets.
We thus believe that the main findings discussed here for the other systems would be equally basis set independent.

%%% FIG 2 %%%
\begin{figure}[h!]
\includegraphics[width=\linewidth]{plot_stat}
	\caption{Nonparallelity errors as function of the number of determinants, for the three classes of CI methods: seniority-based CI (blue), excitation-based CI (red), and our proposed hybrid hCI (green).
	}
        \label{fig:plot_stat}
\end{figure}
%%% %%% %%%

%We start by discussing the dissociation of \ce{F2}, which involves a single bond breaking.
%Now moving to a more challenging problem, the dissociation of \ce{N2}, where three bonds are broken.
%For different CI approaches, Fig.~\ref{fig:N2_pes} shows PECs and their differences with respect to FCI, as well as the NPE and distance errors.
%The associated differences with respect to the FCI result can be seen in the Supporting Information.

%Unless stated otherwise, from here on the performance of each method is probed by their NPE. Later we discuss other metrics.
For all systems (specially ethylene and \ce{N2}), hCI2 is better than CISD, two methods where the number of determinants scales as $N^4$.
hCI2.5 is better than CISDT (except for \ce{H8}), despite its lower computational cost, whereas hCI3 is much better than CISDT, and comparable in accuracy with CISDTQ (again for all systems).
Inspection of the PECs (see SI) reveal that the lower NPE in the hCI results stem mostly from the contribution of the dissociation region.
This result demonstrates the importance of higher-order excitations with low seniority number in this strong correlation regime, 
which are accounted for in hCI but not in excitation-based CI (for a given scaling with the number of determinants).
%The situation at the Franck-Condon region will be discussed later.
 
Meanwhile, the first level of seniority-based CI (sCI0, which is the same as doubly-occupied CI\cite{}) tends to offer a rather low NPE when compared to the other CI methods with a similar number of determinants (hCI2.5 and CISDT).
However, convergence is clearly slower for the next levels (sCI2 and sCI4), whereas excitation-based CI and specially hCI methods converge faster.
Furthermore, seniority-based CI becomes less atractive for larger basis set in view of its exponential scaling.
This can be seen in Fig.Sx of the SI, which shows that seniority-based CI has a much more pronounced increase in the number of determinants for larger basis sets.

It is worth mentioning the surprisingly good performance of the hCI1 and hCI1.5 methods. 
For \ce{HF}, \ce{F2}, and ethylene, they presented lower NPEs than the much more expensive CISDT method, being slightly higher in the case of \ce{N2}.
For the same systems we also see the NPEs increase from hCI1.5 to hCI2, and decreasing to lower values only at the hCI3 level. 
(Even than, it is important to remember that the hCI2 results remain overall superior to their excitation-based counterparts.)
Both findings are not observed for \ce{H4} and \ce{H8}.
It seems that both the relative worsening of the hCI2 method and success of hCI1 and hCI1.5 methods 
become less apparent as progressively more bonds are being broken (compare for instance \ce{F2}, \ce{N2}, and \ce{H8} in Fig.~\ref{fig:plot_stat}).
This reflects the fact that higher-order excitations are needed to properly describe multiple bond breaking.
%Whereas in excitation-based CI, the NPE always decrease as one moves to higher orders.

In Fig.Sx of the SI we present the distance error, which is also found to decrease faster with the hCI methods.
Most of observations discussed for the NPE also hold for the distance error, with two main differences.
The convergence is always monotonic for the latter observable (which is expected from the definition of the observable),
and the performance of seniority-based CI is much poorer (due to the slow recovery of dynamic correlation).


\subsection{Equilibrium geometries and vibrational frequencies}

In Fig.~\ref{fig:xe} and \ref{fig:freq}, we present the convergence of the equilibrium geometries and vibrational frequencies, respectively,
as functions of the number of determinants, for the three classes of CI methods.
%, vibrational frequencies, and dissociation energies,
%For \ce{HF}, \ce{F2}, \ce{N2}, and ethylene,
%For both observables, the overall performance of hCI either exceeds or is comparable to that of excitation-based CI, 
For the equilibrium geometries, hCI performs slightly better overall than excitation-based CI.
%, being better for \ce{F2}, ethylene, and \ce{N2}, and 
A more significant advantage of hCI can be seen for the vibrational frequencies.
For both observables, hCI and excitation-based CI largely outperform seniority-based CI.

Similarly to what we observed for the NPEs, the convergence of hCI was also found to be non-monotonic in some cases.
This oscillatory behavior is particularly evident for \ce{F2}, also noticeable for \ce{HF}, becoming less aparent for ethylene, virtually absent for \ce{N2},
and showing up again for \ce{H4} and \ce{H8}.
Interstingly, equilibrium geometries and vibrational frequencies of \ce{HF} and \ce{F2} (single bond breaking), 
are rather accurate when evaluated at the hCI1.5 level, bearing in mind its relatively modest computational cost.

%%% FIG 3 %%%
\begin{figure}[h!]
	\includegraphics[width=\linewidth]{xe}
	\caption{Equilibrium geometries as function of the number of determinants, for the three classes of CI methods: seniority-based CI (blue), excitation-based CI (red), and our proposed hybrid hCI (green).
        }
        \label{fig:xe}
\end{figure}
%%% %%% %%%

%%% FIG 4 %%%
\begin{figure}[h!]
\includegraphics[width=\linewidth]{freq}
	\caption{Vibrational frequencies (or force constants) as function of the number of determinants, for the three classes of CI methods: seniority-based CI (blue), excitation-based CI (red), and our proposed hybrid hCI (green).
	}
        \label{fig:freq}
\end{figure}
%%% %%% %%%


\subsection{Orbital optimized configuration interaction}

Up to this point, all results and discussions have been based on CI calculations with HF orbitals.
Now we discuss the role of further optimizing the orbitals at each given CI calculation.
Due to the significantly higher computational cost and numerical difficulties for optimizing the orbitals at higher levels of CI,
such calculations were typically limited up to oo-CISD (for excitation-based), oo-DOCI (for seniority-based), and oo-hCI2 (for hCI).
The PECs and analogous results to those of Figs. \ref{fig:plot_stat}, \ref{fig:xe}, and ~\ref{fig:freq} are shown in the SI.

At a given CI level, orbital optimization will lead to lower energies than with HF orbitals. 
However, even though the energy is lowered (thus improved) at each geometry, such improvement may vary largely along the PEC, which may or may not decrease the NPE.
More often than not, the NPEs do decrease upon orbital optimization, though not always.
%For example, oo-hCI1 presents smaller NPEs than their non-optimized counterparts for \ce{N2}, \ce{H4}, \ce{H8}, but not for \ce{HF}, \ce{F2}, and ethylene.
For example, compared with their non-optimized counterparts, oo-hCI1 and oo-hCI1.5 provide somewhat larger NPEs for \ce{HF} and \ce{F2},
% oo-hCI2
comparable for ethylene, and smaller for \ce{N2}, \ce{H4}, and \ce{H8}.
Following the same trend, oo-CISD presents smaller NPEs than HF-CISD for the multiple bond breaking systems, but very similar ones for the single bond breaking cases.
oo-CIS has significantly smaller NPEs than HF-CIS, being comparable to oo-hCI1 for all systems except for \ce{H4} and \ce{H8}. We will come back to oo-CIS latter.
Based on the present oo-CI results, hCI still has the upper hand when compared with excitation-based CI, though by a much smaller margin.
%It does, however, lead to a more monotonic convergence in the case of hCI, which is not necessarily an advantage.
%Also, hCI is slightly better than excitation-based CI for HF orbitals, whereas both are equally good with orbital optimization.
Orbital optimization also reduces the NPE for seniority-based CI (in this case we only considered oo-DOCI).
The gain is specially noticeable for \ce{H4} and \ce{H8}, and much less so for \ce{HF}, ethylene, and \ce{N2}.
For \ce{F2}, we found the interesting situation where orbital optimization actually increases the NPE (though by a small amount).
\fk{in progress...}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion and perspectives}
\label{sec:ccl}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\fk{in progress...}

Here we have proposed a new scheme for truncating the Hilbert space in configuration interaction calculations, named hCI.
By merging the excitation degree and the seniority number into a single hierarchy parameter, 
the hCI method ensures that all classes of determinants sharing the same scaling with the number of electrons are included in each level of the hierarchy.
We evaluated the performance of hCI against the traditional excitation-based CI and seniority-based CI,
by comparing PECs for dissociation of 6 systems, going from single to multiple bond breaking.

Our key finding is that the overall performance of hCI either surpasses or equals those of excitation-based CI and seniority-based CI.
The superiority of hCI methods is more noticeable for the nonparallelity errors, but and also be seen for the equilibrium geometries and vibrational frequencies.
For the challenging cases of \ce{H4} and \ce{H8}, hCI and excitation-based CI perform similarly.
We also found surprisingly good performances for the first level of hCI (hCI1) and the orbital optimized version of CIS (oo-CIS),
given their very favourable computational scaling.

One important conclusion is that orbital optimization is not necessarily a recommended strategy, depending on the properties one is interested in.
One should bear in mind that the optimizing the orbitals is always accompanied with well-known challenges (several solutions, convergence issues)
and may imply in a significant computational burden (associated with the calculations of the orbital gradient, Hessian, and the many iterations that are often required).
In this sense, stepping up in the CI hierarchy might be a more straightforward and possibly cheaper alternative than optimizing the orbitals.
One interesting possibility to explore is to first employ a low order CI method to optimize the orbitals, and then to employ this set of orbitals at a higher level of CI.

The hCI pathways presented here offers several interesting possibilities to explore.
One is to investigate the performance of hCI or some adaptation of it for excited states.
Another is to develop coupled cluster methods based on an analogous hybrid excitation-seniority truncation of the excitation operator.
One could also test the performance of hCI wave functions for Quantum Monte Carlo simulations.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{acknowledgements}
This work was performed using HPC resources from CALMIP (Toulouse) under allocation 2021-18005.
This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No.~863481).
\end{acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Supporting information available}
\label{sec:SI}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PECs, energy differences with respect to FCI results, NPE, closeness errors, equilibrium geometries, vibrational frequencies,
according to the three classes of CI methods (excitation-based CI, seniority-based CI, and hierarchy-CI),
computed for \ce{HF}, \ce{F2}, ethylene, \ce{N2}, \ce{H4}, and \ce{H8}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section*{Data availability statement}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The data that support the findings of this study are openly available in Zenodo at \href{http://doi.org/XX.XXXX/zenodo.XXXXXXX}{http://doi.org/XX.XXXX/zenodo.XXXXXXX}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{seniority}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\end{document}