Merge branch 'master' of git.irsamc.ups-tlse.fr:loos/Ec

This commit is contained in:
Anthony Scemama 2021-07-05 10:57:32 +02:00
commit bafd18ce0f
2 changed files with 91116 additions and 47183 deletions

133322
Ec.nb

File diff suppressed because it is too large Load Diff

View File

@ -27,8 +27,13 @@
\newcommand{\br}{\boldsymbol{r}}
\newcommand{\bc}{\boldsymbol{c}}
\newcommand{\bC}{\boldsymbol{C}}
\newcommand{\bR}{\boldsymbol{R}}
\newcommand{\bX}{\boldsymbol{X}}
\newcommand{\hH}{\Hat{H}}
\newcommand{\hh}{\Hat{h}}
\newcommand{\hX}{\Hat{X}}
\newcommand{\hR}{\Hat{R}}
\usepackage[
colorlinks=true,
@ -65,7 +70,7 @@ In the continuity of our recent work on the benzene molecule [\href{https://doi.
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.
In particular, we study the convergence properties of i) the M{\o}ller-Plesset perturbation series up to fifth-order (MP2, MP3, MP4, and MP5), ii) the iterative approximate single-reference coupled-cluster series CC2, CC3, and CC4, and ii) the single-reference coupled-cluster series CCSD, CCSDT, and CCSDTQ.
In particular, we study the convergence properties of i) the M{\o}ller-Plesset perturbation series up to fifth-order (MP2, MP3, MP4, and MP5), ii) the iterative approximate single-reference coupled-cluster series CC2, CC3, and CC4, and iii) the single-reference coupled-cluster series CCSD, CCSDT, and CCSDTQ.
The performance of the ground-state gold standard CCSD(T) is also investigated.
\end{abstract}
@ -79,7 +84,7 @@ Electronic structure theory relies heavily on approximations.
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 nuclei coordinates can then be treated as parameters in the electronic Hamiltonian.
The second central approximation which makes calculations feasable 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.
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}^n \Hat{T}_k$ (where $n$ is the number of electrons).
@ -140,10 +145,10 @@ The CIPSI calculations have been performed with {\QP}. \cite{Garniron_2019}
\titou{The particularity of the current implementation is that the selection step and the PT2 correction are computed \textit{simultaneously} via a hybrid semistochastic algorithm \cite{Garniron_2017,Garniron_2019} (which explains the statistical error associated with the PT2 correction in the following).
Moreover, a renormalized version of the PT2 correction (dubbed rPT2 below) 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 rPT2 correction and the CIPSI algorithm.}
For all these calculations, we consider Dunning's correlation-consistent double-$\zeta$ basis (cc-pVDZ) which consists of Hilbert space sizes ranging from $10^{20}$ (for thiophene) to $10^{36}$ (for benzene).
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 four distinct sets: HF orbitals, natural orbitals (NOs), localized orbitals (LOs), and optimized orbitals (OOs).
In the present study, we investigate the convergence behavior of the CIPSI energy for four distinct sets: natural orbitals (NOs), localized orbitals (LOs), 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.
Then, natural orbitals (NOs) are computed based on this wave function, and subsequently localized orbitals.
The Boys-Foster localization procedure \cite{Boys_1960} that we apply to the natural orbitals is performed in several orbital windows: \titou{i) core, ii) valence $\sigma$, iii) valence $\pi$, iv) valence $\pi^*$, v) valence $\sigma^*$, vi) the higher-lying $\sigma$ orbitals, and vii) the higher-lying $\pi$ orbitals.}
@ -152,6 +157,8 @@ Because they take advantage of the local character of electron correlation, loca
Using these localized orbitals as starting point, we also perform successive orbital optimizations, 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.
As we shall see below, employing optimized orbitals has the advantage to produce a 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.
We have found that $\expval*{\Hat{S}^2}$ is, nonetheless, very close to zero for each system.
%%%%%%%%%%%%%%%%%%%%%%%%%
@ -175,25 +182,85 @@ where $\Psi_\text{var}^{(k)} = \sum_{I \in \mathcal{I}_k} c_I^{(k)} \ket*{I}$ is
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 CI matrix in the reference space and the magnitude of $E_\text{PT2}$ provides a qualitative idea of the ``distance'' to the FCI limit.
We then linearly extrapolate, using large variational space, the CIPSI energy to $E_\text{PT2} = 0$ (which effectively corresponds to 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.
We have found that $\expval*{\Hat{S}^2}$ is, nonetheless, very close to zero for each system.
From a more general point of view, the variational energy $E_\text{var}^{(k)}$ depends on both the coefficient $\{ c_I \}_{1 \le I \le \Ndet^{(k)}}$ but also on the orbital coefficients $\{C_{\mu p}\}_{1 \le \mu,p \le \Norb,1 \le \mu,p \le \Norb}$ such that the $p$th orbital is
\begin{equation}
\phi_p(\br) = \sum_{\mu} C_{\mu p} \chi_{\mu}(\br)
\end{equation}
where $\chi_{\mu}(\br)$ is a basis function.
Orbital optimization techniques at the SCI level are theoretically straightforward, but practically challenging.
Here, we detail our orbital optimization procedure with the CIPSI algorithm.
From a more general point of view, the variational energy $E_\text{var}^{(k)}$ depends on both the coefficient $\{ c_I \}_{1 \le I \le \Ndet^{(k)}}$ but also on the orbital rotation parameter $\{X_{pq}\}_{1 \le \mu p,q \le \Norb}$.
%such that the $p$th orbital is
%\begin{equation}
% \phi_p(\br) = \sum_{\mu} C_{\mu p} \chi_{\mu}(\br)
%\end{equation}
%where $\chi_{\mu}(\br)$ is a basis function.
The diagonalization of the CI matrix ensures that
\begin{equation}
\pdv{E_\text{var}(\bc,\bC)}{c_I} = 0,
\pdv{E_\text{var}(\bc,\bX)}{c_I} = 0,
\end{equation}
but, a priori, we have
\begin{equation}
\pdv{E_\text{var}(\bc,\bC)}{C_{\mu p}} \neq 0,
\pdv{E_\text{var}(\bc,\bX)}{X_{pq}} \neq 0,
\end{equation}
Here, we assume that the variational wave function is normalized, \ie, $\braket*{\Psi_\text{var}}{\Psi_\text{var}} = 1$.
Then, the previous equation can be rewritten,
\begin{equation}
E(\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}$.
The energy $E$ depends on the determinants and their number, thus the orbital optimization will not be the same for two different sets of determinants.
To understand the concept of orbital rotation, we look at this operator $\bX$ in more details,
\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}),
\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}
\begin{equation}
E = \sum_{pq} h_p^q \gamma_p^q + \frac{1}{2} \sum_{pqrs} v_{pq}^{rs} \Gamma_{pq}^{rs},
\end{equation}
where $\boldsymbol{\gamma}$ is the one-electron density matrix with elements
\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*}
g_{pq}
&= \left. \pdv{E(\bm{c,X})}{X_{pq}}\right|_{\bm{X}=0}
\\
&= \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}
\\
&= \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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results and discussion}
@ -226,7 +293,7 @@ but, a priori, we have
CCSD(T) & $-193.5439$ & $-735.6$ & $-229.4073$ & $-764.0$ & $-225.6099$ & $-774.5$ & $-209.5836$ & $-754.9$ & $-552.0458$ & $-724.8$
\\
\hline
CIPSI & & & & & & & & & & \\
CIPSI & & $-739.3$ & & $-768.1$ & & $-778.3$ & & $-758.4$ & & $-729.1$\\
\end{tabular}
\end{ruledtabular}
\end{table*}
@ -259,7 +326,7 @@ but, a priori, we have
\hline
CCSD(T) & $-231.5798$ & $-857.5$ & $-263.6024$ & $-899.4$ & $-263.5740$ & $-904.1$ & $-247.5929$ & $-877.7$ & $-263.6099$ & $-896.2$ & $-295.5680$ & $-952.2$ & $-279.6305$ & $-913.1$ \\
\hline
CIPSI & & & & & & & & & & \\
CIPSI & & $-863.0$ & & & & $-908.8$ & & $-883.4$ & & $-900.4$ & & $-957.3$ & & $-918.5$\\
\end{tabular}
\end{ruledtabular}
\end{table*}