%aimed at recovering both static and dynamic correlation,
Here we propose a novel partitioning of the Hilbert space, hierarchy configuration interaction (hCI),
where the degree of excitation (with respect to a given reference) and the seniority number (number of unpaired electrons) are combined in a single hierarchy parameter.
The key appealing feature of hCI is that it includes all classes of determinants that share the same scaling with the number of electrons and basis functions.
In this way, it accounts for low-seniority high-excitation determinants lacking in excitation-based CI, while keeping the same computational scaling with system size.
By surveying the dissociation of multiple molecular systems, we examined how fast hCI and their excitation-based and seniority-based parents converge as
we step up towards the exact full CI limit.
We found that the overall performance of hCI usually exceeds or at least parallels that of excitation-based CI.
For small systems and basis sets, doubly-occupied CI (the first level of seniority-based CI) often remains the best option.
However, for larger systems or basis sets, and for higher accuracy, seniority-based CI becomes impractical.
However, some of its interesting features, particularly the small non-parallelity errors, are partially recovered with hCI, at only a polynomical cost.
We have futher explored the role of optimizing the orbitals at several levels of CI.
For higher orders of hCI and excitation-based CI,
the additional computational burden and other known issues related to orbital optimization usually do not compensate the marginal improvements often observed,
when compared with results obtained with canonical Hartree-Fock orbitals.
The exception is orbital-optimized CI with single excitations, a minimally correlated model displaying the qualitatively correct description of single bond breaking,
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
Importantly, the number of determinants $N_{det}$ (which is the key parameter governing the computational cost) scales polynomially with the number of electrons $N$ as $N^{2d}$.
By truncating at the seniority zero ($s =0$) sector, one obtains the doubly occupied CI (DOCI) method \cite{Bytautas_2011,Allen_1962,Smith_1965,Veillard_1967},
which 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}.
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.
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.
\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.}
We know that low degree excitations and low seniority sectors, when looked at individually, often have the most important contribution to the FCI expansion.
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$.
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.
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.
Bytautas et al.\cite{Bytautas_2015} explored a different hybrid scheme combining determinants from a complete active space and with a maximum seniority number.
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.
We evaluated the convergence of four observables: the nonparalellity error (NPE), the distance error, the equilibrium geometries, and the vibrational frequencies.
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},
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.
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).
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.
%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 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 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.
\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).
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).
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.
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},
Results for \ce{HF} with larger basis sets (see Fig.Sx in the \SI) show very similar convergence behaviours, though with less oscillations for the hCI methods.
\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).
\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).
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},
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.
Orbital optimization usually reduces the NPE for seniority-based CI (in this case we only considered oo-DOCI) as well.
The gain is specially noticeable for \ce{H4} and \ce{H8} (where the orbitals become symmetry-broken \cite{}),
and much less so for \ce{HF}, ethylene, and \ce{N2} (where the orbitals remain symmetry-preserving).
This is in line with what has been observed before for \ce{N2}\cite{Bytautas_2011}.
For \ce{F2}, we found that orbital optimization actually increases the NPE (though by a small amount),
due to the larger energy lowering at the Franck-Condon region than at the dissociating region.
These results suggest that, when bond breaking involves one site, orbital optimization at the DOCI level does not have such an important role,
at least in the sense of decreasing the NPE.
Optimizing the orbitals at the CI level also tends to benefit the convergence of vibrational frequencies and equilibrium geometries (shown in Fig.Sx of the \SI).
The impact is often somewhat larger for hCI than for excitation-based CI, by a small margin.
The large oscillations observed in the hCI convergence with HF orbitals (for \ce{HF} and \ce{F2}) are significantly suppressed upon orbital optimization.
We come back to the surprisingly good performance of oo-CIS, which is interesting due to its low computational cost.
The PECs are compared with those of HF and FCI in Fig.Sx of the \SI.
At this level, the orbital rotations provide an optimized reference (different from the HF solution), from which only single excitations are performed.
Since the reference is not the HF one, Briluoin's theorem no longer holds, and single excitations actually connect with the reference.
Thus, with only single excitations (and a referece that is optimized in the presence of these excitations), one obtains a minimally correlated model.
Surprisingly, oo-CIS recovers a non-negligible fraction (15\%-40\%) of the correlation energy around the equilibrium geometries.
For all systems, significantly more correlation energy (25\%-65\% of the total) is recovered at dissociation.
In fact, the larger account of correlation at dissociation implies in the relatively small NPEs encountered at the oo-CIS level.
We also found that the NPE drops more significantly (with respect to the HF one) for the single bond breaking cases (\ce{HF} and \ce{F2}),
followed by the double (ethylene) and triple (\ce{N2}) bond breaking, then \ce{H4}, and finally \ce{H8}.
The above findings can be understood by looking at the character of the oo-CIS orbitals.
At dissociation, the closed-shell reference is actually ionic, with orbitals assuming localized atomic-like characters.
The reference has a decreasing weight in the CI expansion as the bond is stretched, becoming virtually zero at dissociaion.
However, that is the reference one needs to achieve the correct open-shell character of the fragments when the single excitations of oo-CIS are accounted for.
Indeed, the most important single excitations promote the electron from the negative to the positive fragment, thus leading to two singly open-shell radicals.
This is enough to obtain a qualitatively correct description of single bond breaking, hence the relatively low NPEs observed for \ce{HF} and \ce{F2}.
In constrast, the oo-CIS method can only explicitly account for one unpaired electron at each fragment, such that multiple bond breaking become insufficiently described.
Nevertheless, double (ethylene) and even triple (\ce{N2}) bond breaking still appear to be reasonably well-described at the oo-CIS level.
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 and derivied quantities (non-parallelity errors, distance errors, vibrational frequencies, and equilibrium geometries)
for six systems, ranging from single to multiple bond breaking.
Our key finding is that the overall performance of hCI either surpasses or equals that of excitation-based CI,
in the sense that convergence with respect to the number of determinants is usually faster.
The superiority of hCI methods is more noticeable for the nonparallelity errors, but also observed to a lesser extent for the vibrational frequencies and equilibrium geometries.
The comparison to seniority-based CI is less trivial.
DOCI (the first level of seniority-based CI) often provides even lower NPEs for a similar number of determinants, but it falls short in describing the other properties investigated here.
If higher accuracy is desired, than the convergence is faster with hCI (and also excitation-based CI) than seniority-based CI, at least for HF orbitals.
Finally, the exponential scaling of seniority-based CI in practice precludes this approach for larger systems and larger basis sets,
while the favourable polynomial scaling and encouraging performance of hCI as an alternative.
We found surprisingly good results for the first level of hCI (hCI1) and the orbital optimized version of CIS (oo-CIS), two methods with very favourable computational scaling.
In particular, oo-CIS correctly describes single bond breaking.
We hope to report on generalizations to excited states in the future.
%For the challenging cases of \ce{H4} and \ce{H8}, hCI and excitation-based CI perform similarly.
%We have also performed orbital optimization at several CI levels,
An important conclusion is that orbital optimization at the CI level is not necessarily a recommended strategy,
given the overall modest improvement in convergence when compared to results with canonical HF orbitals.
One should bear in mind that 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),
One interesting possibility to explore is to first optimize the orbital at a lower level of CI, and then to employ this set of orbitals at a higher level of CI.
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).
%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}.
