saving work in theory section

This commit is contained in:
Pierre-Francois Loos 2021-07-19 13:57:42 +02:00
parent 07e6f4e4e8
commit 57d06f35f5
2 changed files with 218 additions and 114 deletions

View File

@ -1,13 +1,86 @@
%% This BibTeX bibliography file was created using BibDesk.
%% http://bibdesk.sourceforge.net/
%% Created for Pierre-Francois Loos at 2021-07-02 10:54:28 +0200
%% Created for Pierre-Francois Loos at 2021-07-19 13:49:57 +0200
%% Saved with string encoding Unicode (UTF-8)
@book{Nocedal_1999,
address = {New York, NY, USA},
author = {Nocedal, Jorge and Wright, Stephen J.},
date-added = {2021-07-19 13:49:42 +0200},
date-modified = {2021-07-19 13:49:49 +0200},
doi = {10.1007/b98874},
issn = {1431-8598},
publisher = {Springer, New York, NY},
title = {{Numerical Optimization}},
year = {1999},
Bdsk-Url-1 = {https://doi.org/10.1007/b98874}}
@article{Henderson_2014a,
author = {Henderson, Thomas M. and Bulik, Ireneusz W. and Stein, Tamar and Scuseria, Gustavo E.},
date-added = {2021-07-19 13:37:56 +0200},
date-modified = {2021-07-19 13:37:56 +0200},
doi = {10.1063/1.4904384},
file = {/home/antoinem/Zotero/storage/D9XMVKVQ/Henderson et al. - 2014 - Seniority-based coupled cluster theory.pdf},
journal = {J. Chem. Phys.},
pages = {244104},
publisher = {{American Institute of Physics}},
title = {Seniority-Based Coupled Cluster Theory},
volume = {141},
year = {2014},
Bdsk-Url-1 = {https://doi.org/10.1063/1.4904384}}
@book{Jensen_2017,
address = {{Chichester, UK ; Hoboken, NJ}},
author = {Jensen, Frank},
date-added = {2021-07-19 09:46:10 +0200},
date-modified = {2021-07-19 09:46:10 +0200},
edition = {Third edition},
isbn = {978-1-118-82599-0},
lccn = {QD455.3.E4 J46 2017},
publisher = {{John Wiley \& Sons}},
title = {Introduction to Computational Chemistry},
year = {2017}}
@book{Helgaker_2013,
author = {T. Helgaker and P. J{\o}rgensen and J. Olsen},
date-added = {2021-07-19 09:45:42 +0200},
date-modified = {2021-07-19 09:45:57 +0200},
owner = {joshua},
publisher = {John Wiley \& Sons, Inc.},
timestamp = {2014.11.24},
title = {Molecular Electronic-Structure Theory},
year = {2013}}
@book{Szabo_1996,
address = {{Mineola, N.Y}},
author = {Szabo, Attila and Ostlund, Neil S.},
date-added = {2021-07-19 09:45:07 +0200},
date-modified = {2021-07-19 09:45:07 +0200},
isbn = {978-0-486-69186-2},
lccn = {QD462 .S95 1996},
publisher = {{Dover Publications}},
title = {Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory},
year = {1996}}
@article{Born_1927,
author = {Born, M. and Oppenheimer, R.},
date-added = {2021-07-19 09:42:46 +0200},
date-modified = {2021-07-19 09:43:17 +0200},
doi = {https://doi.org/10.1002/andp.19273892002},
journal = {Ann. Phys.},
number = {20},
pages = {457-484},
title = {Zur Quantentheorie der Molekeln},
volume = {389},
year = {1927},
Bdsk-Url-1 = {https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.19273892002},
Bdsk-Url-2 = {https://doi.org/10.1002/andp.19273892002}}
@misc{Chilkuri_2021,
archiveprefix = {arXiv},
author = {Vijay Gopal Chilkuri and Thomas Applencourt and Kevin Gasperich and Pierre-Fran{\c c}ois Loos and Anthony Scemama},

View File

@ -1,5 +1,5 @@
\documentclass[aps,prb,reprint,noshowkeys,superscriptaddress]{revtex4-1}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,wrapfig,txfonts}
\usepackage{graphicx,dcolumn,bm,xcolor,microtype,multirow,amscd,amsmath,amssymb,amsfonts,physics,wrapfig,txfonts,siunitx}
\usepackage[version=4]{mhchem}
\newcommand{\ie}{\textit{i.e.}}
@ -13,29 +13,51 @@
\newcommand{\trashAS}[1]{\textcolor{green}{\sout{#1}}}
\newcommand{\AS}[1]{\toto{(\underline{\bf AS}: #1)}}
\newcommand{\Ec}{E_\text{c}}
\newcommand{\mEh}{$mE_h$}
\newcommand{\Eh}{$E_h$}
% useful shortcuts
\newcommand{\mc}{\multicolumn}
\newcommand{\fnm}{\footnotemark}
\newcommand{\fnt}{\footnotetext}
\newcommand{\tabc}[1]{\multicolumn{1}{c}{#1}}
\newcommand{\QP}{\textsc{quantum package}}
% key quantities
\newcommand{\Nel}{n}
\newcommand{\Norb}{N}
\newcommand{\Ndet}{N_\text{det}}
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bc}{\boldsymbol{c}}
\newcommand{\bR}{\boldsymbol{R}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bH}{\boldsymbol{H}}
% operators
\newcommand{\hT}{\Hat{T}}
\newcommand{\hH}{\Hat{H}}
\newcommand{\hh}{\Hat{h}}
\newcommand{\hX}{\Hat{X}}
\newcommand{\hR}{\Hat{R}}
\newcommand{\hk}{\Hat{\kappa}}
\newcommand{\cre}[1]{\Hat{a}_{#1}^\dagger}
\newcommand{\ani}[1]{\Hat{a}_{#1}}
% vectors and matrices
\newcommand{\bO}{\boldsymbol{0}}
\newcommand{\bI}{\boldsymbol{1}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\bH}{\boldsymbol{H}}
\newcommand{\bg}{\boldsymbol{g}}
\newcommand{\bk}{\boldsymbol{\kappa}}
\newcommand{\bc}{\boldsymbol{c}}
\newcommand{\br}{\boldsymbol{r}}
% curly letters
\newcommand{\cA}{\mathcal{A}}
\newcommand{\cI}{\mathcal{I}}
\newcommand{\cP}{\mathcal{P}}
% energies
\newcommand{\Ec}{E_\text{c}}
\newcommand{\Evar}{E_\text{var}}
\newcommand{\EPT}{E_\text{PT2}}
\newcommand{\ECIPSI}{E_\text{CIPSI}}
% wave functions and orbitals
\newcommand{\PsiO}{\Psi_0}
\newcommand{\Psivar}{\Psi_\text{var}}
\newcommand{\MO}[1]{\phi_{#1}}
\usepackage[
colorlinks=true,
@ -67,7 +89,6 @@
% Abstract
\begin{abstract}
%We report (frozen-core) full configuration interaction (FCI) energies in finite Hilbert spaces for various five- and six-membered rings.
In the continuity of our recent work on the benzene molecule [\href{https://doi.org/10.1063/5.0027617}{J. Chem. Phys. \textbf{153}, 176101 (2020)}], itself motivated by the blind challenge of Eriksen \textit{et al.} [\href{https://doi.org/10.1021/acs.jpclett.0c02621}{J. Phys. Chem. Lett. \textbf{11}, 8922 (2020)}] on the same system, we report reference frozen-core correlation energies for twelve five- and six-membered ring molecules (cyclopentadiene, furan, imidazole, pyrrole, thiophene, benzene, pyrazine, pyridazine, pyridine, pyrimidine, tetrazine, and triazine) in the standard correlation-consistent double-$\zeta$ Dunning basis set (cc-pVDZ).
This corresponds to Hilbert spaces with sizes ranging from $10^{28}$ (for thiophene) to $10^{36}$ (for benzene).
Our estimates are based on energetically optimized-orbital selected configuration interaction (SCI) calculations performed with the \textit{Configuration Interaction using a Perturbative Selection made Iteratively} (CIPSI) algorithm.
@ -82,18 +103,18 @@ The performance of the ground-state gold standard CCSD(T) is also investigated.
%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%
Electronic structure theory relies heavily on approximations.
Electronic structure theory relies heavily on approximations. \cite{Szabo_1996,Helgaker_2013,Jensen_2017}
Loosely speaking, to make any theory useful, three main approximations must be enforced.
The first fundamental approximation, known as the Born-Oppenheimer approximation, usually consists in assuming that the motion of nuclei and electrons are decoupled.
The first fundamental approximation, known as the Born-Oppenheimer approximation, usually consists in assuming that the motion of nuclei and electrons are decoupled. \cite{Born_1927}
The nuclei coordinates can then be treated as parameters in the electronic Hamiltonian.
The second central approximation which makes calculations feasible by a computer is the basis set approximation where one introduces a set of pre-defined basis functions to represent the many-electron wave function of the system.
In most molecular calculations, a set of one-electron, atom-centered gaussian basis functions are introduced to expand the so-called one-electron molecular orbitals which are then used to build the many-electron Slater determinants.
The third and most relevant approximation in the present context is the ansatz (or form) of the electronic wave function $\Psi$.
For example, in configuration interaction (CI) methods, the wave function is expanded as a linear combination of Slater determinants, while in (single-reference) coupled-cluster (CC) theory, \cite{Cizek_1966,Paldus_1972,Crawford_2000,Piecuch_2002,Bartlett_2007,Shavitt_2009} a reference Slater determinant $\Psi_0$ [usually taken as the Hartree-Fock (HF) wave function] is multiplied by a wave operator defined as the exponentiated excitation operator $\Hat{T} = \sum_{k=1}^\Nel \Hat{T}_k$ (where $\Nel$ is the number of electrons).
For example, in configuration interaction (CI) methods, the wave function is expanded as a linear combination of Slater determinants, while in (single-reference) coupled-cluster (CC) theory, \cite{Cizek_1966,Paldus_1972,Crawford_2000,Piecuch_2002,Bartlett_2007,Shavitt_2009} a reference Slater determinant $\PsiO$ [usually taken as the Hartree-Fock (HF) wave function] is multiplied by a wave operator defined as the exponentiated excitation operator $\hT = \sum_{k=1}^\Nel \hT_k$ (where $\Nel$ is the number of electrons).
The truncation of $\Hat{T}$ allows to define a hierarchy of non-variational and size-extensive methods with improved accuracy:
The truncation of $\hT$ allows to define a hierarchy of non-variational and size-extensive methods with improved accuracy:
CC with singles and doubles (CCSD), \cite{Cizek_1966,Purvis_1982} CC with singles, doubles, and triples (CCSDT), \cite{Noga_1987a,Scuseria_1988} CC with singles, doubles, triples, and quadruples (CCSDTQ), \cite{Oliphant_1991,Kucharski_1992} with corresponding computational scalings of $\order*{\Norb^{6}}$, $\order*{\Norb^{8}}$, and $\order*{\Norb^{10}}$, respectively (where $\Norb$ denotes the number of orbitals).
Parallel to the ``complete'' CC series presented above, an alternative series of approximate iterative CC models have been developed by the Aarhus group in the context of CC response theory \cite{Christiansen_1998} where one skips the most expensive terms and avoids the storage of the higher-excitation amplitudes: CC2, \cite{Christiansen_1995a} CC3, \cite{Christiansen_1995b,Koch_1997} and CC4 \cite{Kallay_2005,Matthews_2021}
Parallel to the ``complete'' CC series presented above, an alternative series of approximate iterative CC models have been developed by the Aarhus group in the context of CC response theory \cite{Christiansen_1998} where one skips the most expensive terms and avoids the storage of the higher-excitation amplitudes: CC2, \cite{Christiansen_1995a} CC3, \cite{Christiansen_1995b,Koch_1997} and CC4. \cite{Kallay_2005,Matthews_2021,Loos_2021}
These iterative methods scale as $\order*{\Norb^{5}}$, $\order*{\Norb^{7}}$, and $\order*{\Norb^{9}}$, respectively, and can be seen as cheaper approximations of CCSD, CCSDT, and CCSDTQ.
Coupled-cluster methods have been particularly successful at computing accurately various properties for small- and medium-sized molecules.
\cite{Kallay_2003,Kallay_2004a,Gauss_2006,Kallay_2006,Gauss_2009}
@ -102,9 +123,9 @@ A similar systematic truncation strategy can be applied to CI methods leading to
Except for full CI (FCI) where all determinants from the Hilbert space (\ie, with excitation degree up to $\Nel$) are considered, truncated CI methods are variational but lack size-consistency.
The non-variationality of truncated CC methods being less of an issue than the size-inconsistency of the truncated CI methods, the formers have naturally overshadowed the latters in the electronic structure landscape.
However, a different strategy has recently made a come back in the context of CI methods. \cite{Bender_1969,Whitten_1969,Huron_1973}
Indeed, selected CI (SCI) methods, \cite{Booth_2009,Giner_2013,Evangelista_2014,Giner_2015,Holmes_2016,Tubman_2016,Liu_2016,Ohtsuka_2017,Zimmerman_2017,Coe_2018,Garniron_2018} where one iteratively selects the energetically relevant determinants from the FCI space (usually) based on a perturbative criterion, has been recently shown to be highly successful in order to produce reference energies for ground and excited states in small- and medium-size molecules \cite{Holmes_2017,Li_2018,Li_2020,Loos_2018a,Chien_2018,Loos_2019,Loos_2020b,Loos_2020c,Loos_2020e,Garniron_2019,Eriksen_2020,Yao_2020,Veril_2021,Loos_2021} thanks to efficient deterministic, stochastic or hybrid algorithms well suited for massive parallelization.
Indeed, selected CI (SCI) methods, \cite{Booth_2009,Giner_2013,Evangelista_2014,Giner_2015,Holmes_2016,Tubman_2016,Liu_2016,Ohtsuka_2017,Zimmerman_2017,Coe_2018,Garniron_2018} where one iteratively selects the important determinants from the FCI space (usually) based on a perturbative criterion, has been recently shown to be highly successful in order to produce reference energies for ground and excited states in small- and medium-sized molecules \cite{Holmes_2017,Li_2018,Li_2020,Loos_2018a,Chien_2018,Loos_2019,Loos_2020b,Loos_2020c,Loos_2020e,Garniron_2019,Eriksen_2020,Yao_2020,Veril_2021,Loos_2021} thanks to efficient deterministic, stochastic or hybrid algorithms well suited for massive parallelization.
We refer the interested reader to Refs.~\onlinecite{Loos_2020a,Eriksen_2021} for recent reviews.
SCI methods are based on a well-known fact: amongst the very large number of determinants belonging to the FCI space, only a relative small fraction of them significantly contributes to the energy.
SCI methods are based on a well-known fact: amongst the very large number of determinants contained in the FCI space, only a relative small fraction of them significantly contributes to the energy.
Accordingly, the SCI+PT2 family of methods performs a sparse exploration of the FCI space by selecting iteratively only the most energetically relevant determinants of the variational space and supplementing it with a second-order perturbative correction (PT2). \cite{Huron_1973,Garniron_2017,Sharma_2017,Garniron_2018,Garniron_2019}
Although the formal scaling of such algorithms remain exponential, the prefactor is greatly reduced which explains their current attractiveness in the electronic structure community and much wider applicability than their standard FCI parent.
Note that, very recently, several groups \cite{Aroeira_2021,Lee_2021,Magoulas_2021} have coupled CC and SCI methods via the externally-corrected CC methodology, \cite{Paldus_2017} showing promising performances for weakly and strongly correlated systems.
@ -141,139 +162,149 @@ The performance of the ground-state gold standard CCSD(T) is also investigated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The geometries of the twelve systems considered in the present study have been all obtained at the CC3/aug-cc-pVTZ level of theory and have been extracted from a previous study. \cite{Loos_2020a}
Note that, for the sake of consistency, the geometry of benzene considered here is different from one of Ref.~\onlinecite{Loos_2020e} which has been computed at a lower level of theory [MP2/6-31G(d)]. \cite{Schreiber_2008}
The MP2, MP3, MP4, CC2, CC3, CC4, CCSD, CCSDT, and CCSDTQ calculations have been performed with CFOUR, \cite{Matthews_2020} while the CCSD(T) and MP5 calculations have been computed with Gaussian 09. \cite{g09}
The CIPSI calculations have been performed with {\QP}. \cite{Garniron_2019}
In the current implementation, the selection step and the PT2 correction are computed simultaneously via a hybrid semistochastic algorithm. \cite{Garniron_2017,Garniron_2019} (which explains the statistical error associated with the PT2 correction in the following).
The MP2, MP3, MP4, CC2, CC3, CC4, CCSD, CCSDT, and CCSDTQ calculations have been performed with CFOUR, \cite{Matthews_2020} while the CCSD(T) and MP5 calculations have been computed with GAUSSIAN 09. \cite{g09}
The CIPSI calculations have been performed with QUANTUM PACKAGE. \cite{Garniron_2019}
In the current implementation, the selection step and the PT2 correction are computed simultaneously via a hybrid semistochastic algorithm \cite{Garniron_2017,Garniron_2019} (which explains the statistical error associated with the PT2 correction in the following).
Here, we employ the renormalized version of the PT2 correction which has been recently implemented and tested for a more efficient extrapolation to the FCI limit thanks to a partial resummation of the higher-order of perturbation. \cite{Garniron_2019}
We refer the interested reader to Ref.~\onlinecite{Garniron_2019} where one can find all the details regarding the implementation of the PT2 correction and the CIPSI algorithm.
For all these calculations, we consider Dunning's correlation-consistent double-$\zeta$ basis (cc-pVDZ).
Although the FCI energy has the enjoyable property of being independent of the set of one-electron orbitals used to construct the many-electron Slater determinants, as a truncated CI method, the convergence properties of CIPSI strongly dependent on this orbital choice.
In the present study, we investigate the convergence behavior of the CIPSI energy for two sets of orbitals in particular: natural orbitals (NOs) and optimized orbitals (OOs).
In the present study, we investigate, in particular, the convergence behavior of the CIPSI energy for two sets of orbitals: natural orbitals (NOs) and optimized orbitals (OOs).
Following our usual procedure, \cite{Scemama_2018,Scemama_2018b,Scemama_2019,Loos_2018a,Loos_2019,Loos_2020a,Loos_2020b,Loos_2020c,Loos_2020e} we perform first a preliminary SCI calculation using HF orbitals in order to generate a SCI wave function with at least $10^7$ determinants.
Natural orbitals are computed based on this wave function and they are used to perform a new CIPSI run.
Successive orbital optimizations are then performed, which consist in minimizing the variational CIPSI energy at each iteration up to approximately $2 \times 10^5$ determinants.
When convergence is achieved in terms of orbital optimization, as our ``production'' run, we perform a new CIPSI calculation from scratch using this set of optimized orbitals.
In some cases, we also explore the use of localized orbitals (LOs) which are produced with the help of the Boys-Foster localization procedure \cite{Boys_1960} that we apply to the natural orbitals in several orbital windows in order to preserve a strict $\sigma$-$\pi$ separation in the planar systems considered here.
Using optimized orbitals has the undeniable advantage to produce, for a given variational energy, more compact CI expansions.
In some cases, we also explore the use of localized orbitals (LOs) which are produced with the Boys-Foster localization procedure \cite{Boys_1960} that we apply to the natural orbitals in several orbital windows in order to preserve a strict $\sigma$-$\pi$ separation in the planar systems considered here. \cite{Loos_2020e}
Because they take advantage of the local character of electron correlation, localized orbitals have been shown to provide faster convergence towards the FCI limit compared to natural orbitals. \cite{Angeli_2003,Angeli_2009,BenAmor_2011,Suaud_2017,Chien_2018,Eriksen_2020,Loos_2020e}
As we shall see below, employing optimized orbitals has the advantage to produce an even smoother and faster convergence of the SCI energy toward the FCI limit.
Note that, unlike excited-state calculations where it is important to enforce that the wave functions are eigenfunctions of the $\Hat{S}^2$ spin operator, \cite{Chilkuri_2021} the present wave functions do not fulfil this property as we aim for the lowest possible energy of a singlet state.
Note that, unlike excited-state calculations where it is important to enforce that the wave functions are eigenfunctions of the $\Hat{S}^2$ spin operator, \cite{Chilkuri_2021} the present wave functions do not fulfill this property as we aim for the lowest possible energy of a closed-shell singlet state.
We have found that $\expval*{\Hat{S}^2}$ is, nonetheless, very close to zero for each system.
%%%%%%%%%%%%%%%%%%%%%%%%%
\section{CIPSI with optimized orbitals}
%%%%%%%%%%%%%%%%%%%%%%%%%
Here, we provide key details about the CIPSI method. \cite{Huron_1973,Garniron_2019}
Here, we provide key details about the CIPSI method \cite{Huron_1973,Garniron_2019} as well as the orbital optimization procedure which has been shown to be highly effective in the context of SHCI by Umrigar and coworkers. \cite{Eriksen_2020,Yao_2020,Yao_2021}
Note that we focus on the ground state but the present discussion can be easily extended to excited states. \cite{Scemama_2019,Veril_2021}
At the $k$th iteration, the total CIPSI energy $E_\text{CIPSI}^{(k)}$ is defined as the sum of the variational energy
At the $k$th iteration, the total CIPSI energy $\ECIPSI^{(k)}$ is defined as the sum of the variational energy
\begin{equation}
E_\text{var}^{(k)} = \frac{\mel*{\Psi_\text{var}^{(k)}}{\Hat{H}}{\Psi_\text{var}^{(k)}}}{\braket*{\Psi_\text{var}^{(k)}}{\Psi_\text{var}^{(k)}}}
\Evar^{(k)} = \frac{\mel*{\Psivar^{(k)}}{\hH}{\Psivar^{(k)}}}{\braket*{\Psivar^{(k)}}{\Psivar^{(k)}}}
\end{equation}
and a second-order perturbative correction
\begin{equation}
E_\text{PT2}^{(k)}
= \sum_{\alpha \in \mathcal{A}_k} e_{\alpha}
= \sum_{\alpha \in \mathcal{A}_k} \frac{\mel*{\Psi_\text{var}^{(k)}}{\Hat{H}}{\alpha}}{E_\text{var}^{(k)} - \mel*{\alpha}{\Hat{H}}{\alpha}}
\EPT^{(k)}
= \sum_{\alpha \in \cA_k} e_{\alpha}
= \sum_{\alpha \in \cA_k} \frac{\mel*{\Psivar^{(k)}}{\hH}{\alpha}}{\Evar^{(k)} - \mel*{\alpha}{\hH}{\alpha}}
\end{equation}
where $\hH$ is the (non-relativistic) electronic Hamiltonian,
\begin{equation}
\label{eq:Psivar}
\Psi_\text{var}^{(k)} = \sum_{I \in \mathcal{I}_k} c_I^{(k)} \ket*{I}
\Psivar^{(k)} = \sum_{I \in \cI_k} c_I^{(k)} \ket*{I}
\end{equation}
is the variational wave function, $\mathcal{I}_k$ is the set of internal determinants $\ket*{I}$ and $\mathcal{A}_k$ is the set of external determinants $\ket*{\alpha}$ which do not belong to the variational space but are linked to it via a nonzero matrix element, \ie, $\mel*{\Psi_\text{var}^{(k)}}{\Hat{H}}{\alpha} \neq 0$.
The sets $\mathcal{I}_k$ and $\mathcal{A}_k$ define, at the $k$th iteration, the internal and external spaces, respectively.
In practice, $E_\text{var}^{(k)}$ is computed by diagonalizing the $\Ndet^{(k)} \times \Ndet^{(k)}$ CI matrix $\bH$ with elements $H_{IJ} = \mel{I}{\hH}{J}$ via Davidson's algorithm \cite{Davidson_1975} and the magnitude of $E_\text{PT2}$ provides a qualitative idea of the ``distance'' to the FCI limit. \cite{Garniron_2018}
We then linearly extrapolate, using large variational space, the CIPSI energy to $E_\text{PT2} = 0$ (which effectively corresponds to the FCI limit).
is the variational wave function, $\cI_k$ is the set of internal determinants $\ket*{I}$ and $\cA_k$ is the set of external determinants $\ket*{\alpha}$ which do not belong to the variational space but are linked to it via a nonzero matrix element, \ie, $\mel*{\Psivar^{(k)}}{\hH}{\alpha} \neq 0$.
The sets $\cI_k$ and $\cA_k$ define, at the $k$th iteration, the internal and external spaces, respectively.
Hereafter, we will label these iterations over the number of determinants $\Ndet^{(k)}$ as \textit{macroiterations}.
In practice, $\Evar^{(k)}$ is computed by diagonalizing the $\Ndet^{(k)} \times \Ndet^{(k)}$ CI matrix $\bH$ with elements $H_{IJ} = \mel{I}{\hH}{J}$ via Davidson's algorithm. \cite{Davidson_1975}
The magnitude of $\EPT^{(k)}$ provides a qualitative idea of the ``distance'' to the FCI limit. \cite{Garniron_2018}
We then linearly extrapolate, using large variational space, the CIPSI energy to $\EPT = 0$ (which effectively corresponds to the FCI limit).
Further details concerning the extrapolation procedure are provided below (see Sec.~\ref{sec:res}).
Orbital optimization techniques at the SCI level are theoretically straightforward, but practically challenging. \cite{Yao_2020,Yao_2021}
Here, we detail our orbital optimization procedure within the CIPSI algorithm and we assume that the variational wave function is normalized, \ie, $\braket*{\Psi_\text{var}}{\Psi_\text{var}} = 1$.
Then, the variational energy can be written as
\begin{equation}
E_\text{var}(\bc,\bX) = \mel{\Psi_\text{var}}{e^{\hX} \hH e^{-\hX}}{\Psi_\text{var}},
\end{equation}
where $\bc$ gathers the CI coefficients, $\bX$ the orbital rotation parameters and $\hX$ is a one-electron anti-hermitian operator, which creates a rotation matrix when exponentiated, \ie, $\bR = e^{\bX}$.
From a more general point of view, the variational energy $E_\text{var}$ depends on both the coefficient $\{ c_I \}_{1 \le I \le \Ndet^{(k)}}$ [see Eq.~\eqref{eq:Psivar}] but also on the orbital rotation parameter $\{X_{pq}\}_{1 \le p,q \le \Norb}$.
Orbital optimization techniques at the SCI level are theoretically straightforward, but practically challenging.
Here, we detail our orbital optimization procedure within the CIPSI algorithm and we assume that the variational wave function is normalized, \ie, $\braket*{\Psivar}{\Psivar} = 1$.
For a given set of orbitals, The diagonalization of the CI matrix ensures that
From a more general point of view, $\Evar$ depends on both the CI coefficients $\{ c_I \}_{1 \le I \le \Ndet^{(k)}}$ [see Eq.~\eqref{eq:Psivar}] but also on the orbital rotation parameters $\{\kappa_{pq}\}_{1 \le p,q \le \Norb}$.
Then, one can conveniently rewrite the variational energy as
\begin{equation}
\pdv{E_\text{var}(\bc,\bX)}{c_I} = 0,
\label{eq:Evar_c_k}
\Evar(\bc,\bk) = \mel{\Psivar}{e^{\hk} \hH e^{-\hk}}{\Psivar},
\end{equation}
but, a priori, we have
where $\bc$ gathers the CI coefficients, $\bk$ the orbital rotation parameters, and
\begin{equation}
\pdv{E_\text{var}(\bc,\bX)}{X_{pq}} \neq 0,
\hk = \sum_{p < q} \sum_{\sigma} \kappa_{pq} \qty(\cre{p\sigma} \ani{q\sigma} - \cre{q\sigma} \ani{p\sigma})
\end{equation}
Although one could use a second order method to minimize the corresponding energy, one has to realize that the size of the CI space is much larger than orbital space.
It is therefore more appropriate to perform a minimization of the variational energy with respect to the orbital rotation parameters and then compute the new CI coefficients by re-diagonalizing the CI matrix.
is a real-valued one-electron anti-hermitian operator, which creates a unitary transformation of the orbital coefficients when exponentiated, $\ani{p\sigma}$ ($\cre{p\sigma}$) being the second quantization annihilation (creation) operator which annihilates (creates) a spin-$\sigma$ electron in the (real-valued) spatial orbital $\MO{p}(\br)$.
To understand the concept of orbital rotation, we look at this operator $\bX$ in more details,
Applying the Newton-Raphson method by Taylor-expanding the variational energy to second order around $\bk = \bO$, \ie,
\begin{equation}
\hX = \sum_{p > q} \sum_{\sigma} X_{pq} (\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma} - \hat{a}_{q \sigma}^{\dagger} \hat{a}_{p \sigma}),
\label{eq:energy_expansion}
\Evar(\bc,\bk) \approx \Evar(\bc,\bO) + \bg \cdot \bk + \frac{1}{2} \bk^{\dag} \cdot \bH \cdot \bk,
\end{equation}
To do so, we need the first and the second derivatives of the energy with respect to the orbital rotations. Also, we would need the coupling of these last ones with the changes in the CI coefficients, but we will see later that this last point is problematic.
Alternatively, the variational energy $E$ can be expressed as a function of the one- and two-electron density matrices \cite{Helgaker2000Aug}
we have
\begin{equation}
E = \sum_{pq} h_p^q \gamma_p^q + \frac{1}{2} \sum_{pqrs} v_{pq}^{rs} \Gamma_{pq}^{rs},
\bk = - \bH^{-1} \cdot \bg,
\end{equation}
where $\boldsymbol{\gamma}$ is the one-electron density matrix with elements
where $\bg$ and $\bH$ are the orbital gradient and Hessian, respectively, both evaluated at $\bk = \bO$.
Their elements are explicitly given by the following expressions: \cite{Henderson_2014a}
\begin{equation}
\gamma_p^q = \sum_{\sigma} \mel{\Psi_\text{var}}{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{}}{\Psi_\text{var}},
\end{equation}
$\boldsymbol{\Gamma}$ is the two-electron density matrix with elements
\begin{equation}
\Gamma_{pq}^{rs} = \sum_{\sigma \sigma'} \mel{\Psi_\text{var}}{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{r \sigma'}^{\dagger} \hat{a}_{s \sigma'}^{} \hat{a}_{q \sigma}^{}}{\Psi_\text{var}},
\end{equation}
$\boldsymbol{h}$ is the core Hamiltonian matrix containing the one-electron integrals with elements
\begin{equation}
h_p^q = \int \phi_p(\br) \hh(\br) \phi_q(\br) d\br,
\end{equation}
and $\boldsymbol{v}$ is the matrix containing the two-electron repulsion integrals with elements
\begin{equation}
v_{pq}^{rs} = \iint \phi_p(\br_1) \phi_q(\br_2) \frac{1}{\abs*{\br_1 - \br_2}} \phi_r(\br_1) \phi_s(\br_2) d\br_1 d\br_2.
\end{equation}
%with $\phi_k$ the $k$-th molecular orbital defined as a linear combination of the atomic orbitals, $\br$ the spatial coordinates of the electron and $\hat{h}$: the one-electron part of the Hamiltonian.
Then, with the differentiation of this variational energy $E$, we obtain the gradient $\bm{g}$ and the Hessian $\bm{H}$. \cite{Henderson2014Dec} The gradient of the energy with respect to the orbital rotation, $g_{pq}$, around $\bm{x} = 0$,
\begin{align*}
\begin{split}
g_{pq}
&= \left. \pdv{E(\bm{c,X})}{X_{pq}}\right|_{\bm{X}=0}
&= \left. \pdv{\Evar(\bc,\bk)}{\kappa_{pq}}\right|_{\bk=\bO}
\\
&= \sum_{\sigma} \mel{\Psi}{\comm*{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{} - \hat{a}_{q \sigma}^{\dagger} \hat{a}_{p \sigma}^{}}{\hH}}{\Psi}
&= \sum_{\sigma} \mel{\Psivar}{\comm*{\cre{p\sigma} \ani{q\sigma} - \cre{q\sigma} \ani{p\sigma}}{\hH}}{\Psivar}
\\
&= \mathcal{P}_{pq} \left[ \sum_r \left( h_p^r \ \gamma_r^q - h_r^q \ \gamma_p^r \right) + \sum_{rst} \left( v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs} \right) \right]
\end{align*}
And the Hessian of the energy with respect to the orbital rotation, $H_{pq,rs}$, around $\bX$ = 0,
\begin{align*}
H_{pq,rs} =
& \left. \pdv{E(\bm{c,X})}{X_{pq}}{X_{rs}}\right|_{\bm{X}=0} \\
=& \ \frac{1}{2} \mathcal{P}_{pq} \mathcal{P}_{rs} \sum_{\sigma, \sigma'} \bra*{\Psi} \comm*{\hat{a}_{r \sigma'}^{\dagger} \hat{a}_{s \sigma'}^{}}{\comm*{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{}}{\hH}} \ket*{\Psi} \\
&+ \frac{1}{2} \mathcal{P}_{pq} \mathcal{P}_{rs} \sum_{\sigma, \sigma'} \bra*{\Psi} \comm*{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{}}{\comm*{\hat{a}_{r \sigma'}^{\dagger} \hat{a}_{s \sigma'}^{}}{\hH}} \ket*{\Psi} \\
=& \ \mathcal{P}_{pq} \mathcal{P}_{rs} \frac{1}{2} \sum_u [\delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) + \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]-(h_p^s \gamma_r^q + h_r^q \gamma_p^s) \\
&+ \frac{1}{2} \sum_{tuv} [\delta_{qr}(v_{pt}^{uv}
\Gamma_{uv}^{st} +v_{uv}^{st} \Gamma_{pt}^{uv}) + \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})] \\
&+ \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs} \Gamma_{ps}^{uv}) - \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts})],
\end{align*}
with $\delta_{ij}$ is the Kronecker delta, $\mathcal{P}_{pq}$ the permutation operator $\mathcal{P}_{pq} = 1
- (p \leftrightarrow q)$ and where $(p \leftrightarrow q)$ applied to an equation returns the same equation with the indices $p$ and $q$ swapped.
&= \cP_{pq} \qty[ \sum_r \left( h_p^r \ \gamma_r^q - h_r^q \ \gamma_p^r \right) + \sum_{rst} \qty( v_{pt}^{rs} \Gamma_{rs}^{qt} - v_{rs}^{qt} \Gamma_{pt}^{rs} ) ]
\end{split}
\end{equation}
and
\begin{widetext}
\begin{equation}
\begin{split}
H_{pq,rs}
& = \left. \pdv{\Evar(\bc,\bk)}{\kappa_{pq}}{\kappa_{rs}}\right|_{\bk=\bO}
\\
& = \cP_{pq} \cP_{rs} \qty{
\frac{1}{2} \sum_{\sigma\sigma'} \mel*{\Psivar}{\comm*{\cre{r \sigma'} \ani{s \sigma'}}{\comm*{\cre{p \sigma} \ani{q \sigma}}{\hH}}}{\Psivar}
+ \frac{1}{2} \sum_{\sigma\sigma'} \mel*{\Psivar}{\comm*{\cre{p \sigma} \ani{q \sigma}}{\comm*{\cre{r \sigma'} \ani{s \sigma'}}{\hH}}}{\Psivar}
}
\\
& = \cP_{pq} \cP_{rs} \Bigg\{
\frac{1}{2} \sum_u \qty[ \delta_{qr}(h_p^u \gamma_u^s + h_u^s \gamma_p^u) + \delta_{ps}(h_r^u \gamma_u^q + h_u^q \gamma_u^r)]
- (h_p^s \gamma_r^q + h_r^q \gamma_p^s)
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} + \frac{1}{2} \sum_{tuv} \qty[ \delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} + v_{uv}^{st} \Gamma_{pt}^{uv})
+ \delta_{ps}(v_{uv}^{qt} \Gamma_{rt}^{uv} + v_{rt}^{uv}\Gamma_{uv}^{qt})]
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} + \sum_{uv} (v_{pr}^{uv} \Gamma_{uv}^{qs} + v_{uv}^{qs} \Gamma_{ps}^{uv})
- \sum_{tu} (v_{pu}^{st} \Gamma_{rt}^{qu}+v_{pu}^{tr} \Gamma_{tr}^{qu}+v_{rt}^{qu}\Gamma_{pu}^{st} + v_{tr}^{qu}\Gamma_{pu}^{ts})]
\Bigg\}
\end{split}
\end{equation}
\end{widetext}
where $\delta_{pq}$ is the Kronecker delta, $\cP_{pq} = 1 - (p \leftrightarrow q)$ is a permutation operator,
\begin{subequations}
\begin{gather}
\gamma_p^q = \sum_{\sigma} \mel{\Psivar}{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{}}{\Psivar},
\\
\Gamma_{pq}^{rs} = \sum_{\sigma \sigma'} \mel{\Psivar}{\cre{p\sigma} \cre{r\sigma'} \ani{s\sigma'} \ani{q\sigma}}{\Psivar},
\end{gather}
\end{subequations}
are the elements of the one- and two-electron density matrices, and
\begin{subequations}
\begin{gather}
h_p^q = \int \MO{p}(\br) \, \hh(\br) \, \MO{q}(\br) d\br,
\\
v_{pq}^{rs} = \iint \MO{p}(\br_1) \MO{q}(\br_2) \frac{1}{\abs*{\br_1 - \br_2}} \MO{r}(\br_1) \MO{s}(\br_2) d\br_1 d\br_2.
\end{gather}
\end{subequations}
are the one- and two-electron integrals, respectively.
Trust region method and Newton's method are very similar, they use a quadratical approximation of a real function, i.e.,
Taylor expansion truncated at the second order. In the Newton method, the step is given by the minimizer of the
quadratical approximation, contrary to the trust region method which gives the minimizer of the quadratical approximation
in the trust region. The trust region defines a region where the quadratical approximation is a adequate representation
of the real function and it evolves during the optimization process in order to preserve the adequacy. The constraint
for the step size is the following, $\norm{\bm{X}_{k+1}} \leq \Delta_k$ with $\Delta_k$ the trust radius. By putting the
constraint with a Lagrange multiplier $\lambda$ and derivating the Lagrangian, the solution is
$\bm{X}_{k+1} = - (\bm{H_k} + \lambda \bm{I})^{-1} \cdot \bm{g}_k$.
The addition of a constant $\lambda \geq 0$ on the diagonal of the hessian removes the negative eigenvalues and
reduce the size of the step since the calculation uses its inverse. By choosing the right $\lambda$ the step size is contraint
into a hypersphere of radius $\Delta_k$. In addition, the evolution of $\Delta_k$ during the optimization and the use of
a condition to cancel a step ensure the convergence of the algorithm.
More details could be found in Numerical Optimization, Nocedal \& Wright (1999)
Because the size of the CI space is much larger than the orbital space, for each macroiteration, we perform multiple \textit{microiterations} which consist in minimizing the variational energy \eqref{eq:Evar_c_k} with respect to the $\Norb(\Norb-1)/2$ independent orbital rotation parameters.
Micoriterations are stopped when a stationary point is found, \ie, $\norm{\bg}_\infty < \tau$, where $\tau$ is a user-defined threshold which has been set to $10^{-4}$ a.u.~in the present study, and a new CIPSI selection step is performed.
%\begin{equation}
% \Evar = \sum_{pq} h_p^q \gamma_p^q + \frac{1}{2} \sum_{pqrs} v_{pq}^{rs} \Gamma_{pq}^{rs},
%\end{equation}
\titou{To enhance convergence, we here employ a variant of the Newton-Raphson method known as ``trust region''. \cite{Nocedal_1999} which defines a region where the quadratic approximation is a adequate representation of the real function and it evolves during the optimization process in order to preserve the adequacy.
The constraint for the step size is the following, $\norm{\bk^{(\ell+1)}} \leq \Delta^{(\ell)}$ with $\Delta^{(\ell)}$ the trust radius at the $\ell$th microiteration.
By putting the constraint with a Lagrange multiplier $\lambda$ and derivating the Lagrangian, the solution is $\bk^{(\ell+1)} = - (\bH^{(\ell)} + \lambda \bI)^{-1} \cdot \bg^{(\ell)}$.
The addition of a constant $\lambda \geq 0$ on the diagonal of the hessian removes the negative eigenvalues and reduce the size of the step since the calculation uses its inverse.
By choosing the right $\lambda$ the step size is constraint into a hypersphere of radius $\Delta^{(\ell)}$.
In addition, the evolution of $\Delta^{(\ell)}$ during the optimization and the use of a condition to cancel a step ensure the convergence of the algorithm.
More details could be found in Ref.~\onlinecite{Nocedal_1999}.}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
@ -281,7 +312,7 @@ More details could be found in Numerical Optimization, Nocedal \& Wright (1999)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{table*}
\caption{Total energy $E$ (in \Eh) and correlation energy $\Delta E$ (in \mEh) for the frozen-core ground state of five-membered rings in the cc-pVDZ basis set.
\caption{Total energy $E$ (in \SI{}{\hartree}) and correlation energy $\Delta E$ (in \SI{}{\milli\hartree}) for the frozen-core ground state of five-membered rings in the cc-pVDZ basis set.
\label{tab:Tab5-VDZ}}
\begin{ruledtabular}
\begin{tabular}{lcccccccccc}
@ -314,7 +345,7 @@ More details could be found in Numerical Optimization, Nocedal \& Wright (1999)
\begin{squeezetable}
\begin{table*}
\caption{Total energy $E$ (in \Eh) and correlation energy $\Delta E$ (in \mEh) for the frozen-core ground state of six-membered rings in the cc-pVDZ basis set.
\caption{Total energy $E$ (in \SI{}{\hartree}) and correlation energy $\Delta E$ (in \SI{}{\milli\hartree}) for the frozen-core ground state of six-membered rings in the cc-pVDZ basis set.
\label{tab:Tab6-VDZ}}
\begin{ruledtabular}
\begin{tabular}{lcccccccccccccc}
@ -361,7 +392,7 @@ More details could be found in Numerical Optimization, Nocedal \& Wright (1999)
\includegraphics[width=0.24\textwidth]{Pyrimidine_EvsNdet}
\includegraphics[width=0.24\textwidth]{Tetrazine_EvsNdet}
\includegraphics[width=0.24\textwidth]{Triazine_EvsNdet}
\caption{$\Delta E_\text{var}$ (solid) and $\Delta E_\text{var} + E_\text{PT2}$ (dashed) as functions of the number of determinants $\Ndet$ in the variational space for the twelve cyclic molecules represented in Fig.~\ref{fig:mol}.
\caption{$\Delta \Evar$ (solid) and $\Delta \Evar + \EPT$ (dashed) as functions of the number of determinants $\Ndet$ in the variational space for the twelve cyclic molecules represented in Fig.~\ref{fig:mol}.
Two sets of orbitals are considered: natural orbitals (NOs, in red) and optimized orbitals (OOs, in blue).
The CCSDTQ correlation energy is represented as a thick black line.
\label{fig:vsNdet}}
@ -382,7 +413,7 @@ More details could be found in Numerical Optimization, Nocedal \& Wright (1999)
\includegraphics[width=0.24\textwidth]{Pyrimidine_EvsPT2}
\includegraphics[width=0.24\textwidth]{Tetrazine_EvsPT2}
\includegraphics[width=0.24\textwidth]{Triazine_EvsPT2}
\caption{$\Delta E_\text{var}$ as a function of $E_\text{PT2}$ for the twelve cyclic molecules represented in Fig.~\ref{fig:mol}.
\caption{$\Delta \Evar$ as a function of $\EPT$ for the twelve cyclic molecules represented in Fig.~\ref{fig:mol}.
Two sets of orbitals are considered: natural orbitals (NOs, in red) and optimized orbitals (OOs, in blue).
The four-point linear fit using the four largest variational wave functions for each set is depicted as a dashed black line.
The CCSDTQ correlation energy is also represented as a thick black line.