saving work

This commit is contained in:
Pierre-Francois Loos 2021-07-19 19:54:44 +02:00
parent 57d06f35f5
commit e419811229

View File

@ -89,7 +89,7 @@
% Abstract
\begin{abstract}
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).
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.
The performance and convergence properties of several series of methods are investigated.
@ -102,6 +102,7 @@ The performance of the ground-state gold standard CCSD(T) is also investigated.
%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{sec:intro}
%%%%%%%%%%%%%%%%%%%%%%%%%
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.
@ -196,7 +197,7 @@ At the $k$th iteration, the total CIPSI energy $\ECIPSI^{(k)}$ is defined as the
and a second-order perturbative correction
\begin{equation}
\EPT^{(k)}
= \sum_{\alpha \in \cA_k} e_{\alpha}
= \sum_{\alpha \in \cA_k} e_{\alpha}^{(k)}
= \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,
@ -204,18 +205,20 @@ where $\hH$ is the (non-relativistic) electronic Hamiltonian,
\label{eq:Psivar}
\Psivar^{(k)} = \sum_{I \in \cI_k} c_I^{(k)} \ket*{I}
\end{equation}
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$.
is the variational wave function, $\cI_k$ is the set of internal determinants $\ket*{I}$ and $\cA_k$ is the set of external determinants (or perturbers) $\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}
The perturbers corresponding to the largest $\abs*{e_{\alpha}^{(k)}}$ values are then added to the variational space at iteration $k+1$.
In our implementation, the size of the variational space is roughly doubled at each iteration.
Hereafter, we label these iterations over $k$ which consist in enlarging the variational space as \textit{macroiterations}.
In practice, $\Evar^{(k)}$ is computed by diagonalizing the $\Ndet^{(k)} \times \Ndet^{(k)}$ CI matrix with elements $\mel{I}{\hH}{J}$ via Davidson's algorithm. \cite{Davidson_1975}
The magnitude of $\EPT^{(k)}$ provides, at iteration $k$, 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.
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$.
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}$.
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}$.
Then, one can conveniently rewrite the variational energy as
\begin{equation}
\label{eq:Evar_c_k}
@ -229,10 +232,10 @@ is a real-valued one-electron anti-hermitian operator, which creates a unitary t
Applying the Newton-Raphson method by Taylor-expanding the variational energy to second order around $\bk = \bO$, \ie,
\begin{equation}
\label{eq:energy_expansion}
\label{eq:EvarTaylor}
\Evar(\bc,\bk) \approx \Evar(\bc,\bO) + \bg \cdot \bk + \frac{1}{2} \bk^{\dag} \cdot \bH \cdot \bk,
\end{equation}
we have
one can iteratively minimize the variational energy with respect to the parameters $\kappa_{pq}$ by setting
\begin{equation}
\bk = - \bH^{-1} \cdot \bg,
\end{equation}
@ -249,30 +252,32 @@ Their elements are explicitly given by the following expressions: \cite{Henderso
\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{
& = \cP_{pq} \cP_{rs} \Bigg\{
\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}
}
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} + \frac{1}{2} \sum_{\sigma\sigma'} \mel*{\Psivar}{\comm*{\cre{p \sigma} \ani{q \sigma}}{\comm*{\cre{r \sigma'} \ani{s \sigma'}}{\hH}}}{\Psivar}
\Bigg\}
\\
& = \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\{} - (h_p^s \gamma_r^q + h_r^q \gamma_p^s)
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} + \frac{1}{2} \sum_{tuv} \delta_{qr}(v_{pt}^{uv} \Gamma_{uv}^{st} + v_{uv}^{st} \Gamma_{pt}^{uv})
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} + \frac{1}{2} \sum_{tuv} \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})]
\\
& \phantom{\cP_{pq} \cP_{rs} \Bigg\{} - \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}
@ -291,20 +296,20 @@ 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 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.
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.
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^{-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.
%\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)}$.
To enhance convergence, we here 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: $\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 differentiating 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}.}
More details can be found in Ref.~\onlinecite{Nocedal_1999}.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}