small modifications in III

Manuscript/Ec.bib

@@ -1,13 +1,79 @@
 %% This BibTeX bibliography file was created using BibDesk.
 %% http://bibdesk.sourceforge.net/
 
-%% Created for Pierre-Francois Loos at 2021-07-21 08:17:37 +0200 
+%% Created for Pierre-Francois Loos at 2021-07-21 13:08:32 +0200 
 
 
 %% Saved with string encoding Unicode (UTF-8) 
 
 
 
+@article{Kreplin_2020,
+	author = {Kreplin,David A. and Knowles,Peter J. and Werner,Hans-Joachim},
+	date-added = {2021-07-21 13:06:31 +0200},
+	date-modified = {2021-07-21 13:08:29 +0200},
+	doi = {10.1063/1.5142241},
+	journal = {J. Chem. Phys.},
+	number = {7},
+	pages = {074102},
+	title = {MCSCF optimization revisited. II. Combined first- and second-order orbital optimization for large molecules},
+	volume = {152},
+	year = {2020},
+	Bdsk-Url-1 = {https://doi.org/10.1063/1.5142241}}
+
+@article{Kreplin_2019,
+	author = {Kreplin,David A. and Knowles,Peter J. and Werner,Hans-Joachim},
+	date-added = {2021-07-21 13:06:07 +0200},
+	date-modified = {2021-07-21 13:07:46 +0200},
+	doi = {10.1063/1.5094644},
+	journal = {J. Chem. Phys.},
+	number = {19},
+	pages = {194106},
+	title = {Second-order MCSCF optimization revisited. I. Improved algorithms for fast and robust second-order CASSCF convergence},
+	volume = {150},
+	year = {2019},
+	Bdsk-Url-1 = {https://doi.org/10.1063/1.5094644}}
+
+@article{Sun_2017,
+	abstract = {We present a new second order complete active space self-consistent field implementation to converge wavefunctions for both large active spaces and large atomic orbital (AO) bases. Our algorithm decouples the active space wavefunction solver from the orbital optimization in the microiterations, and thus may be easily combined with various modern active space solvers. We also introduce efficient approximate orbital gradient and Hessian updates, and step size determination. We demonstrate its capabilities by calculating the low-lying states of the Fe(II)-porphine complex with modest resources using a density matrix renormalization group solver in a CAS(22,27) active space and a 3000 AO basis.},
+	author = {Qiming Sun and Jun Yang and Garnet Kin-Lic Chan},
+	date-added = {2021-07-21 13:05:05 +0200},
+	date-modified = {2021-07-21 13:05:21 +0200},
+	doi = {https://doi.org/10.1016/j.cplett.2017.03.004},
+	journal = {Chem. Phys. Lett.},
+	pages = {291-299},
+	title = {A general second order complete active space self-consistent-field solver for large-scale systems},
+	volume = {683},
+	year = {2017},
+	Bdsk-Url-1 = {https://www.sciencedirect.com/science/article/pii/S0009261417302166},
+	Bdsk-Url-2 = {https://doi.org/10.1016/j.cplett.2017.03.004}}
+
+@article{Werner_1985,
+	author = {Werner,Hans‐Joachim and Knowles,Peter J.},
+	date-added = {2021-07-21 13:04:44 +0200},
+	date-modified = {2021-07-21 13:07:18 +0200},
+	doi = {10.1063/1.448627},
+	journal = {J. Chem. Phys.},
+	number = {11},
+	pages = {5053-5063},
+	title = {A second order multiconfiguration SCF procedure with optimum convergence},
+	volume = {82},
+	year = {1985},
+	Bdsk-Url-1 = {https://doi.org/10.1063/1.448627}}
+
+@article{Werner_1980,
+	author = {Werner,Hans‐Joachim and Meyer,Wilfried},
+	date-added = {2021-07-21 13:04:17 +0200},
+	date-modified = {2021-07-21 13:06:47 +0200},
+	doi = {10.1063/1.440384},
+	journal = {J. Chem. Phys.},
+	number = {5},
+	pages = {2342-2356},
+	title = {A quadratically convergent multiconfiguration--self‐consistent field method with simultaneous optimization of orbitals and CI coefficients},
+	volume = {73},
+	year = {1980},
+	Bdsk-Url-1 = {https://doi.org/10.1063/1.440384}}
+
 @article{Magoulas_2021,
 	author = {Magoulas, Ilias and Gururangan, Karthik and Piecuch, Piotr and Deustua, J. Emiliano and Shen, Jun},
 	date-added = {2021-07-21 08:17:09 +0200},

Manuscript/Ec.tex

@@ -210,7 +210,7 @@ and a second-order perturbative correction
 \begin{equation}
 	\EPT^{(k)}
 	= \sum_{\alpha \in \cA_k} e_{\alpha}^{(k)}
-	= \sum_{\alpha \in \cA_k} \frac{\abs*{\mel*{\Psivar^{(k)}}{\hH}{\alpha}}^2}{\Evar^{(k)} - \mel*{\alpha}{\hH}{\alpha}}
+	= \sum_{\alpha \in \cA_k} \frac{\abs*{\mel*{\Psivar^{(k)}}{\hH}{\alpha}}^2}{\Evar^{(k)} - \mel*{\alpha}{\hH}{\alpha}},
 \end{equation}
 where $\hH$ is the (non-relativistic) electronic Hamiltonian,
 \begin{equation}
@@ -228,10 +228,12 @@ We then linearly extrapolate, using large variational wave functions, the CIPSI
 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.
+Most of the technology presented here has been borrowed from complete-active-space self-consistent-field (CASSCF) methods \cite{Werner_1980,Werner_1985,Sun_2017,Kreplin_2019,Kreplin_2020} but one of the strength of SCI methods is that one does not need to select an active space and to classify orbitals as active, inactive, and virtual orbitals.
 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$.
 
 As stated in Sec.~\ref{sec:intro}, $\Evar$ depends on both the CI coefficients $\{ c_I \}_{1 \le I \le \Ndet}$ [see Eq.~\eqref{eq:Psivar}] but also on the orbital rotation parameters $\{\kappa_{pq}\}_{1 \le p,q \le \Norb}$.
-Here, motivated by cost saving arguments, we have chosen to optimize separately the CI and orbital coefficients by alternatively diagonalizing the CI matrix after each selection step and then rotating the orbitals until the variational energy for a given number of determinants is minimal.
+Motivated by cost saving arguments, we have chosen to optimize separately the CI and orbital coefficients by alternatively diagonalizing the CI matrix after each selection step and then rotating the orbitals until the variational energy for a given number of determinants is minimal.
+(For a detailed comparison of coupled, uncoupled, and partially-coupled optimizations within SCI methods, we refer the interested reader to the recent work of Yao and Umrigar. \cite{Yao_2021})
 To do so, we conveniently rewrite the variational energy as
 \begin{equation}
 \label{eq:Evar_c_k}
@@ -241,7 +243,7 @@ where $\bc$ gathers the CI coefficients, $\bk$ the orbital rotation parameters,
 \begin{equation}
 	\hk = \sum_{p < q} \sum_{\sigma} \kappa_{pq} \qty(\cre{p\sigma} \ani{q\sigma} - \cre{q\sigma} \ani{p\sigma})
 \end{equation}
-is a real-valued one-electron antisymmetric operator, which creates an orthogonal 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)$. \cite{Helgaker_2013}
+is a real-valued one-electron antisymmetric operator, which creates an orthogonal 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)$. \cite{Helgaker_2013}
 
 Applying the Newton-Raphson method by Taylor-expanding the variational energy to second order around $\bk = \bO$, \ie,
 \begin{equation}
@@ -250,6 +252,7 @@ Applying the Newton-Raphson method by Taylor-expanding the variational energy to
 \end{equation}
 one can iteratively minimize the variational energy with respect to the parameters $\kappa_{pq}$ by setting
 \begin{equation}
+	\label{eq:kappa_newton}
 	\bk = - \bH^{-1} \cdot \bg,
 \end{equation}
 where $\bg$ and $\bH$ are the orbital gradient and Hessian, respectively, both evaluated at $\bk = \bO$.
@@ -294,8 +297,10 @@ and
 where $\delta_{pq}$ is the Kronecker delta, $\cP_{pq} = 1 - (p \leftrightarrow q)$ is a permutation operator,
 \begin{subequations}
 \begin{gather}
+	\label{eq:one_dm}
 	\gamma_p^q = \sum_{\sigma} \mel{\Psivar}{\hat{a}_{p \sigma}^{\dagger} \hat{a}_{q \sigma}^{}}{\Psivar},
 	\\
+	\label{eq:two_dm}
 	\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}
@@ -311,20 +316,23 @@ are the elements of the one- and two-electron density matrices, and
 \end{subequations}
 are the one- and two-electron integrals, respectively.
 
-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 iteratively minimizing the variational energy \eqref{eq:Evar_c_k} with respect to the $\Norb(\Norb-1)/2$ independent orbital rotation parameters.
+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 iteratively minimizing the variational energy \eqref{eq:Evar_c_k} with respect to the $\Norb(\Norb-1)/2$ independent orbital rotation parameters for a fixed set of determinants.
+After each microiteration (\ie, orbital rotation), the one- and two-electron integrals [see Eqs.~\eqref{eq:one} and \eqref{eq:two}] have to be updated. 
+Moreover, the CI matrix must be re-diagonalized and new one- and two-electron density matrices [see Eqs.~\eqref{eq:one_dm} and \eqref{eq:two_dm}] are computed.
 Microiterations 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^{-3}$ a.u.~in the present study, and a new CIPSI selection step is performed.
 Note that a tight convergence is not critical here as a new set of microiterations is performed at each macroiteration and a new production CIPSI run is performed from scratch using the final set of orbitals.
-It is also worth pointing out that, after each orbital rotation, the one- and two-electron integrals defined in Eqs.~\eqref{eq:one} and \eqref{eq:two} have to be updated for the next iteration.
+
+This procedure might sound computationally expensive but one has to realize that it is usually performed only for relatively compact variational space 
 %\begin{equation}
 %	\Evar =  \sum_{pq} h_p^q \gamma_p^q + \frac{1}{2} \sum_{pqrs} v_{pq}^{rs} \Gamma_{pq}^{rs},
 %\end{equation}
 
 To enhance the convergence of the microiteration process, we employ a variant of the Newton-Raphson method known as ``trust region''. \cite{Nocedal_1999}
-This popular variant defines a region where the quadratic approximation \eqref{eq:EvarTaylor} is an adequate representation of the objective energy function \eqref{eq:Evar_c_k} and it evolves during the optimization process in order to preserve the adequacy via a constraint on the step size preventing it from overstepping, \ie, $\norm{\bk^{(\ell+1)}} \leq \Delta^{(\ell)}$, where $\Delta^{(\ell)}$ is the trust radius at the $\ell$th microiteration.
-By introducing a Lagrange multiplier $\lambda$ to control the trust-region size, one obtains $\bk^{(\ell+1)} = - (\bH^{(\ell)} + \lambda \bI)^{-1} \cdot \bg^{(\ell)}$.
-The addition of the level shift $\lambda \geq 0$ removes the negative eigenvalues and ensure the positive definiteness of the Hessian matrix by reducing the step size.
-By choosing the right value of $\lambda$, the step size is constrained into a hypersphere of radius $\Delta^{(\ell)}$ and is able to evolve from the Newton direction at $\lambda = 0$ to the steepest descent direction as $\lambda$ grows.
-The evolution of the trust radius during the optimization and the use of a condition to cancel the step where the energy rises ensure the convergence of the algorithm.
+This popular variant defines a region where the quadratic approximation \eqref{eq:EvarTaylor} is an adequate representation of the objective energy function \eqref{eq:Evar_c_k} and it evolves during the optimization process in order to preserve the adequacy via a constraint on the step size preventing it from overstepping, \ie, $\norm{\bk} \leq \Delta$, where $\Delta$ is the trust radius.
+By introducing a Lagrange multiplier $\lambda$ to control the trust-region size, one replaces Eq.~\eqref{eq:kappa_newton} by $\bk = - (\bH + \lambda \bI)^{-1} \cdot \bg$.
+The addition of the level shift $\lambda \geq 0$ removes the negative eigenvalues and ensures the positive definiteness of the Hessian matrix by reducing the step size.
+By choosing the right value of $\lambda$, the step size is constrained into a hypersphere of radius $\Delta$ and is able to evolve from the Newton direction at $\lambda = 0$ to the steepest descent direction as $\lambda$ grows.
+The evolution of the trust radius during the optimization and the use of a condition to cancel the step when the energy rises ensure the convergence of the algorithm.
 More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -332,6 +340,7 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 \label{sec:res}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+%%% FIG 2 %%%
 \begin{figure*}
 	\includegraphics[width=0.24\textwidth]{Cyclopentadiene_EvsNdet}
 	\includegraphics[width=0.24\textwidth]{Furan_EvsNdet}
@@ -352,7 +361,9 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 	The CCSDTQ correlation energy is represented as a thick black line.
 	\label{fig:vsNdet}}
 \end{figure*}
+%%% %%% %%%
 
+%%% FIG 3 %%%
 \begin{figure*}
 	\includegraphics[width=0.24\textwidth]{Cyclopentadiene_EvsPT2}
 	\includegraphics[width=0.24\textwidth]{Furan_EvsPT2}
@@ -375,7 +386,9 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 	The CCSDTQ correlation energy is also represented as a thick black line.
 	\label{fig:vsEPT2}}
 \end{figure*}
+%%% %%% %%%
 
+%%% TABLE I %%%
 \begin{squeezetable}
 \begin{table*}
 	\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.
@@ -410,7 +423,9 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 	\end{ruledtabular}
 \end{table*}
 \end{squeezetable}
+%%% %%% %%%
 
+%%% TABLE II %%%
 \begin{squeezetable}
 \begin{table*}
 	\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.
@@ -445,8 +460,10 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 	\end{ruledtabular}
 \end{table*}
 \end{squeezetable}
+%%% %%% %%%
 
 
+%%% TABLE III %%%
 \begin{squeezetable}
 \begin{table}
 	\caption{
@@ -547,13 +564,14 @@ More details can be found in Ref.~\onlinecite{Nocedal_1999}.
 	\end{ruledtabular}
 \end{table}
 \end{squeezetable}
+%%% %%% %%%
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \subsection{CIPSI estimates}
 \label{sec:cipsi_res}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-We first study the convergence of the CIPSI correlation energy $\Delta \Evar = \Evar - \EHF$ as a function of the number of determinants.
+We first study the convergence of the CIPSI correlation energy $\Delta \Evar = \Evar - \EHF$ (where $\EHF$ is the HF energy) as a function of the number of determinants.
 Our motivation here is to generate FCI-quality reference correlation energies for the twelve cyclic molecules represented in Fig.~\ref{fig:mol} in order to benchmark, in a second time, the performance and convergence properties of various mainstream MP and CC methods (see Sec.~\ref{sec:mpcc_res}).
 For the natural and optimized orbital sets, we report, in Fig.~\ref{fig:vsNdet}, the evolution of the variational correlation energy $\Delta \Evar$ and its perturbatively corrected value $\Delta \Evar + \EPT$ with respect to the number of determinants for each cyclic molecule.
 As compared to natural orbitals (solid red lines), one can see that, for a given number of determinants, the use of optimized orbitals greatly lowers $\Delta \Evar$ (solid blue lines).
@@ -561,8 +579,9 @@ Adding the perturbative correction $\EPT$ yields similar curves for both sets of
 This indicates that, for a given number of determinants, $\EPT$ (which provides a qualitative idea to the distance to the FCI limit) is much smaller for optimized orbitals than for natural orbitals.
 This is further evidenced in Fig.~\ref{fig:vsEPT2} where we show the behavior of $\Delta \Evar$ as a function of $\EPT$ for both sets of orbitals.
 The four-point weighted linear fit using the four largest variational wave functions are also represented (dashed black lines), while the CCSDTQ correlation energy (solid black line) is also reported for comparison purposes.
-From Fig.~\ref{fig:vsEPT2}, it is clear that the behavior of $\Delta \Evar$ is much more linear and produces smaller $\EPT$ values, hence facilitating the extrapolation procedure to the FCI limit.
+From Fig.~\ref{fig:vsEPT2}, it is clear that the behavior of $\Delta \Evar$ is much more linear and produces smaller $\EPT$ values, hence facilitating the extrapolation procedure to the FCI limit (see below).
 
+%%% FIG 4 %%%
 \begin{figure}
 	\includegraphics[width=\linewidth]{Benzene_EvsNdetLO}
 	\caption{$\Delta \Evar$ (solid) and $\Delta \Evar + \EPT$ (dashed) as functions of the number of determinants $\Ndet$ in the variational space for the benzene molecule.
@@ -570,13 +589,14 @@ From Fig.~\ref{fig:vsEPT2}, it is clear that the behavior of $\Delta \Evar$ is m
 	The CCSDTQ correlation energy is represented as a thick black line.
 	\label{fig:BenzenevsNdet}}
 \end{figure}
+%%% %%% %%%
 
-Figure~\ref{fig:BenzenevsNdet} compares the convergence of $\Delta \Evar$ for the natural, localized, and optimized sets of orbitals.
+Figure~\ref{fig:BenzenevsNdet} compares the convergence of $\Delta \Evar$ for the natural, localized, and optimized sets of orbitals in the case of benzene.
 As mentioned in Sec.~\ref{sec:compdet}, although both sets break the spatial symmetry to take advantage of the local nature of electron correlation, optimized orbitals further improve on the use of localized orbitals.
 More quantitatively, optimized orbitals produce the same variational energy as localized orbitals with, roughly, a ten-fold reduction in the number of determinants.
 A similar improvement is observed going from natural to localized orbitals.
-\titou{T2: comment on PT2 for localized orbitals?}
-Accordingly to these observations, all our FCI correlation energy estimates have been produced with the set of optmized orbitals.
+\titou{Comment on PT2 for localized orbitals.}
+Accordingly to these observations, all our FCI correlation energy estimates have been produced with the set of optimized orbitals.
 
 To do so, we have then extrapolated to the orbital-optimized CIPSI calculations to $\EPT = 0$ via a weighted linear extrapolation using the four largest variational wave functions.
 These estimates are reported in Tables \ref{tab:Tab5-VDZ} and \ref{tab:Tab6-VDZ} for the five- and six-membered rings, respectively.
@@ -586,6 +606,7 @@ These estimates are reported in Tables \ref{tab:Tab5-VDZ} and \ref{tab:Tab6-VDZ}
 \label{sec:mpcc_res}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+%%% FIG 5 %%%
 \begin{figure*}
 	\includegraphics[width=0.32\textwidth]{Cyclopentadiene_MPCC}
 	\includegraphics[width=0.32\textwidth]{Furan_MPCC}
@@ -607,6 +628,7 @@ These estimates are reported in Tables \ref{tab:Tab5-VDZ} and \ref{tab:Tab6-VDZ}
 	The CIPSI estimate of the correlation energy is represented as a black line.
 	\label{fig:MPCC}}
 \end{figure*}
+%%% %%% %%%
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%