#  ${\Delta\text{NO}}$ and the complexities of electron correlation in simple
hydrogen clusters

Ismael A. Elayan  Department of Chemistry, University of Winnipeg, Winnipeg,
Manitoba, R3B 2E9, Canada  Rishabh Gupta  Department of Chemistry, University
of Winnipeg, Winnipeg, Manitoba, R3B 2E9, Canada  Joshua W. Hollett
<EMAIL_ADDRESS>Department of Chemistry, University of Winnipeg,
Winnipeg, Manitoba, R3B 2E9, Canada  Department of Chemistry, University of
Manitoba, Winnipeg, Manitoba, R3T 2N2, Canada

###### Abstract

The $\Delta$ natural orbital ($\Delta\text{NO}$) two-electron density matrix
(2-RDM) and energy expression are derived from a multideterminantal wave
function. The approximate $\Delta\text{NO}$ 2-RDM is combined with an on-top
density functional and a double-counting correction to capture electron
correlation. A trust-region Newton’s method optimization algorithm for the
simultaneous optimization of $\Delta\text{NO}$ orbitals and occupancies is
introduced and compared to the previous iterative diagonalization algorithm.
The combination of $\Delta\text{NO}$ and two different on-top density
functionals, Colle-Salvetti (CS) and Opposite-spin exponential cusp and Fermi-
hole correction (OF), is assessed on small hydrogen clusters and compared to
density functional, single-reference coupled cluster, and multireference
perturbation theory (MRMP2) methods. The $\Delta\text{NO-CS}$ and
$\Delta\text{NO-OF}$ methods outperform the single-reference methods, and are
comparable to MRMP2. However, there is a distinct qualitative error in the
$\Delta\text{NO}$ potential energy surface for H4 compared to the exact. This
discrepancy is explained through analysis of the $\Delta\text{NO}$ orbitals,
occupancies and the two-electron density.

††preprint: Hn v1

##  I Introduction

In an effort to create efficient electron structure methods, it has become
commonplace to combine contrasting models of electron correlation (e.g. wave
function, density functional, density matrix functional) in order to exploit
their most effective attributes.Cizek (1969); Jeziorski and Monkhorst (1981);
Andersson  _et al._ (1990); Mok, Neumann, and Handy (1996); Handy and Cohen
(2001); Cremer (2001); Pulay (2011); Becke (2013); Crittenden (2013);
Tsuchimochi and Voorhis (2014); Wallace and Crittenden (2014); Hollett,
Hosseini, and Menzies (2016); Ramos-Cordoba, Salvador, and Matito (2016);
Benavides-Riveros, Lathiotakis, and Marques (2017); Evangelista (2018); Via-
Nadal  _et al._ (2019); Fromager, Toulouse, and Jørgen Aa. Jensen (2007);
Weimer, Sala, and Görling (2008); Toulouse  _et al._ (2009); Janesko,
Henderson, and Scuseria (2009); Chai and Head-Gordon (2009); Stoyanova  _et
al._ (2013); Yamanaka  _et al._ (2006); Rohr, Toulouse, and Pernal (2010);
Piris (2017, 2018); Pérez-Jiménez, Pérez-Jordá, and Sancho-Garcia (2007);
Manni  _et al._ (2014); Carlson, Truhlar, and Gagliardi (2015); Gagliardi  _et
al._ (2017); Malcolm and McDouall (1998); Köhn  _et al._ (2013); Hubert,
Hedegård, and H. J (2016); Hollett and Loos (2019); Hirao (1992a); Grimme and
Waletzke (1999); Gräfenstein and Cremer (2000); Angeli, Cimiraglia, and
Malrieu (2001, 2002); Mostafanejad, Liebenthal, and DePrince III (2020);
Furche (2008) Current methods that employ such a strategy, vary in
computational cost and accuracy, and in most cases their remains room for
improvement in both departments. The most common strategy is to decompose
electron correlation into static (strong) and dynamic (weak) correlation, and
approach each differently. The development of these methods involves the
design of models that precisely capture each component, and also blend the
approaches together in a seamless fashion. The former aspect can be thoroughly
evaluated by applying the method to relatively simple systems, such as small
hydrogen clusters.

Although they are not chemically complex, the electronic structure of H atom
clusters, particularly linear chains and regular lattices of H atoms, can
serve as a systematic testing ground for new treatments of electron
correlation.Paldus  _et al._ (1993); Kowalski and Jankowski (1998a, b);
Schipper, Gritsenko, and Baerends (1999); Jankowski and Kowalski (1999); Van
Voorhis and Head-Gordon (2000); Hachmann, Cardoen, and Chan (2006); Al-Saidi,
Zhang, and Krakauer (2007); Rassolov and Xu (2007); Sinitskiy, Greenman, and
Mazziotti (2010); Stella  _et al._ (2011); Bytautas  _et al._ (2011); Ramos-
Cordoba  _et al._ (2015); Bulik, Henderson, and Scuseria (2015); Robinson and
Knowles (2012a, b, c, d); Sand and Mazziotti (2013); Limacher  _et al._
(2013); Kats and Manby (2013); Henderson  _et al._ (2014); Motta  _et al._
(2017, 2020); Giner  _et al._ (2020); Marie, Kossoski, and Loos (2021)
Creating a model that accurately describes the surprisingly complex features
of H cluster potential energy surfaces can prove challenging. There have been
several studies of such clusters with an array of methods, including density
functional approximations, natural orbital functionals, and truncated post-HF
methods. These systems have been shown to be effective for evaluating the
ability of methods to describe both ground and excited state properties. While
it is important for these methods to accurately describe the potential energy
surface of these “simple” systems, it is also important to accurately describe
the underlying two-electron density. A method that correctly models the two-
electron density, as well as the energy and one-electron density, is more
likely to be universally successful.

The $\Delta\text{NO}$ method combines a cumulant functional approach for the
treatment of static correlation with an on-top density functional (or post-HF)
to treat dynamic correlation.Hollett, Hosseini, and Menzies (2016); Hollett
and Loos (2019) Amongst other approaches, it most resembles
multiconfigurational-pair density funcitonal theory (MC-PDFT).Manni  _et al._
(2014); Carlson, Truhlar, and Gagliardi (2015); Gagliardi  _et al._ (2017);
Sand, Truhlar, and Gagliardi (2017); Mostafanejad, Liebenthal, and DePrince
III (2020) In which, a multireference one-electron density matrix (1-RDM), and
subsequent one-electron density, is used to determine the kinetic, nuclear-
potential, and classical Coulomb energy. The on-top pair density, obtained
from the multireference two-electron density matrix (2-RDM), is used to create
a modified spin-densityManni  _et al._ (2014) which is inserted into the usual
exchange-correlation density functionals (e.g. PBEPerdew, Burke, and Ernzerhof
(1996)) to determine the exchange and correlation energy. The 2-RDM can be
obtained from a traditional multireference self-consistent-field Manni  _et
al._ (2014); Carlson, Truhlar, and Gagliardi (2015); Gagliardi  _et al._
(2017); Sand, Truhlar, and Gagliardi (2017), a variational 2-RDM-driven
complete active-space self-consistent fieldMostafanejad, Liebenthal, and
DePrince III (2020), or another economical alternative such as a pair-coupled
cluster doubles calculation.Garza  _et al._ (2015) The $\Delta\text{NO}$
method uses an on-top density functional for only the dynamic correlation
energy, where the functional depends directly on the on-top density. The
remaining components of the energy are determined using the approximate
$\Delta\text{NO}$ 2-RDM.

This work begins by introducing a $\Delta\text{NO}$ wave function (Subsection
II.1) from which the $\Delta\text{NO}$ two-electron density matrix and energy
expression can be derived. The $\Delta\text{NO}$ 2-RDM is then supplemented
with an ad hoc term to correct for missing static correlation between electron
pairs (Subsection II.2). This is followed by a description of the
determination of the dynamic correlation energy from the $\Delta\text{NO}$
2-RDM (Subsection II.3), and then an approach to reduce correlation double-
counting (Subsection II.4). Then a new optimization algorithm for the
$\Delta\text{NO}$ 2-RDM is described (Subsection II.5), which involves the
simultaneous optimization of the orbitals and the electron transfer variables.
The performance of the optimization method is then evaluated relative to the
previous $\Delta\text{NO}$ 2-RDM optimization algorithm. The performance of
$\Delta\text{NO}$, in describing potential energy surfaces of small H
clusters, is then compared to established single- and multireference methods
(Section IV). The two-electron density is then used to analyze the static
correlation provided by the $\Delta\text{NO}$ model in the rectangular H4
cluster (Subsection IV.7). Finally, the results of the $\Delta\text{NO}$
evaluation and analysis are summarized (Section V).

##  II Theory

###  II.1 The $\Delta\text{NO}$ wave function

The multideterminantal $\Delta\text{NO}$ wave function,

$\Psi_{\Delta\text{NO}}=\sum_{I}c_{I}\Psi_{I},$ (1)

is defined with respect to a reference configuration in which all open-shell
(singly occupied) orbitals are considered frozen (i.e. singly occupied in all
determinants). Only doubly occupied orbitals are active and all excitations
are $2n$-tuples (doubles, quadruples, …) where no electrons are unpaired. In
other words, with respect to the active orbitals, $\Psi_{\Delta\text{NO}}$ is
seniority-zero.Bytautas  _et al._ (2011); Stein, Henderson, and Scuseria
(2014); Henderson  _et al._ (2014); Poelmans  _et al._ (2015); Chen, Zhou, and
Wu (2015); Limacher (2015); Van Meer, Gritsenko, and Baerends (2018); Brzek
_et al._ (2019); Vu and DePrince III (2020) Also, each active pair of
electrons has its own active space, consisting of the doubly occupied orbital
in the reference wave function and some number of vacant (virtual) orbitals.
Each determinant of the expansion is then defined by $2n$-tuple excitations of
the reference configuration,

$\Psi_{I}=\Psi_{i\bar{i}j\bar{j}k\bar{k}\dots}^{a\bar{a}b\bar{b}c\bar{c}\dots}$
(2)

where $i,j$ and $k$ refer to reference occupied orbitals, and $a,b$ and $c$
refer to all active orbitals. The absence of a bar, $i$, indicates an
$\alpha$-spin orbital while the presence of a bar, $\bar{i}$, indicates a
$\beta$-spin orbital. Considering this, it is useful to define the
determinantal expansion coefficients, $c_{I}$, as a product of coefficients
corresponding to specific double excitations of the active electron pairs,

$c_{I}=d_{i}^{a}d_{j}^{b}d_{k}^{c}\dots$ (3)

In order to arrive at the $\Delta\text{NO}$ formalism, the double excitation
coefficients, $d_{i}^{a}$, are defined in terms of orbital occupancies,
$n_{i}$, and electron transfer variables, $\Delta_{ia}$,

$d_{i}^{a}=\delta_{ia}\sqrt{n_{i}}-\sqrt{\Delta_{ia}},$ (4)

where $\delta_{ia}$ is the Kroenecker delta and $\Delta_{ii}=0$. Orbital
occupancies, occupied and vacant (with respect to the reference
configuration), are defined as

$\displaystyle n_{i}$ $\displaystyle=1-\sum_{a}\Delta_{ia},$ (5)
$\displaystyle n_{a}$ $\displaystyle=\sum_{i}\Delta_{ia}.$ (6)

The $N$-representable components of the $\Delta\text{NO}$ energy expression,
$E^{\Delta\text{NO}^{\prime}}$, can then be obtained by expanding the
expectation value of the $N$-electron hamiltonian in terms of $n_{i}$,
$\Delta_{ia}$, and one and two-electron integrals over orbitals. Initially,
the energy may be divided into diagonal and off-diagonal contributions,

$\displaystyle E^{\Delta\text{NO}^{\prime}}$
$\displaystyle=\langle\Psi_{\Delta\text{NO}}|\hat{H}|\Psi_{\Delta\text{NO}}\rangle$
$\displaystyle=\sum_{I}c_{I}^{2}H_{II}+\sum_{I\neq J}c_{I}c_{J}H_{IJ},$ (7)

where

$H_{IJ}=\langle\Psi_{I}|\hat{H}|\Psi_{J}\rangle,$ (8)

and the electronic hamiltonian (in atomic units) is given by

$\hat{H}=\sum_{i=1}^{N}-\tfrac{1}{2}\nabla^{2}_{i}-\sum_{i=1}^{N}\sum_{A=1}^{M}\frac{Z_{A}}{r_{iA}}+\sum_{i=2}^{N}\sum_{j<i}\frac{1}{r_{ij}},$
(9)

where $M$ is the number of nuclei, $Z_{A}$ is the nuclear charge, $r_{iA}$ is
the electron-nucleus distance, and $r_{ij}$ is the electron-electron distance.
The diagonal contributions can be expressed as,

$\sum_{I}c_{I}^{2}H_{II}=E^{\Delta\text{NO}}_{\text{0-1RDM}}+E^{\Delta\text{NO}}_{\text{pair}},$
(10)

where $E^{\Delta\text{NO}}_{\text{0-1RDM}}$ is the energy associated with the
zeroth-order term in the expansion of the two-electron density matrix (2-RDM)
in terms of the one-electron density matrix (1-RDM),Hollett and Loos (2019)
which is the basis of cumulant functional theory.Kutzelnigg and Mukherjee
(1999) The 0-1RDM energy is expressed as,

$\displaystyle E^{\Delta\text{NO}}_{\text{0-1RDM}}$
$\displaystyle=\sum^{\text{cl}}_{p}2n_{p}h_{pp}+\sum^{\text{op}}_{p}n_{p}h_{pp}$
$\displaystyle+\sum^{\text{cl}}_{p}\sum^{\text{cl}}_{q}n_{p}n_{q}\left(2J_{pq}-K_{pq}\right)$
$\displaystyle+\sum^{\text{cl}}_{p}\sum^{\text{op}}_{q}n_{p}n_{q}\left(2J_{pq}-K_{pq}\right)$
$\displaystyle+\sum^{\text{op}}_{p}\sum^{\text{op}}_{q}\frac{n_{p}n_{q}}{2}\left(J_{pq}-K_{pq}\right),$
(11)

where cl and op denote the set of closed and open-shell orbitals,
respectively. The one-electron integrals over spatial orbitals, $\phi_{p}$,
are defined as,

$h_{pp}=\int\phi^{*}_{p}(\mathbf{r})\left[-\frac{1}{2}\nabla^{2}-\sum_{A}\frac{Z_{A}}{r_{A}}\right]\phi^{*}_{p}(\mathbf{r})d\mathbf{r}$
(12)

where $\mathbf{r}=(x,y,z)$ and the two-electron Coulomb, $J_{pq}$, and
exchange, $K_{pq}$, integrals are defined as,

$J_{pq}=\int\frac{\phi^{*}_{p}(\mathbf{r}_{1})\phi^{*}_{q}(\mathbf{r}_{2})\phi_{p}(\mathbf{r}_{1})\phi_{q}(\mathbf{r}_{2})}{r_{12}}d\mathbf{r}_{1}d\mathbf{r}_{2}$
(13)

and

$K_{pq}=\int\frac{\phi^{*}_{p}(\mathbf{r}_{1})\phi^{*}_{q}(\mathbf{r}_{2})\phi_{q}(\mathbf{r}_{1})\phi_{p}(\mathbf{r}_{2})}{r_{12}}d\mathbf{r}_{1}d\mathbf{r}_{2}.$
(14)

After expanding the expectation values for the diagonal contributions to
$E^{\Delta\text{NO}^{\prime}}$ [first term of Equation (II.1)], simplifying,
and removing the contribution of $E^{\Delta\text{NO}}_{\text{0-1RDM}}$, the
pair correction energy, $E^{\Delta\text{NO}}_{\text{pair}}$, remains. The pair
correction energy is the first term associated with the cumulant and is given
by,

$E^{\Delta\text{NO}}_{\text{pair}}=\sum_{pq}\eta_{pq}\left(2J_{pq}-K_{pq}\right),$
(15)

where

$\displaystyle\eta_{pq}$ $\displaystyle=\delta_{pq}n_{p}(1-n_{p})$
$\displaystyle+\left(1-\delta_{pq}\right)\Bigg{[}O_{p}V_{q}\Delta_{pq}\left(n_{q}-n_{p}-\Delta_{pq}\right)$
$\displaystyle\hphantom{\left(1-\delta_{pq}\right)\Bigg{[}}+V_{p}O_{q}\Delta_{qp}\left(n_{p}-n_{q}-\Delta_{qp}\right)$
$\displaystyle\hphantom{\left(1-\delta_{pq}\right)\Bigg{[}}-V_{p}V_{q}\sum_{r}\Delta_{rp}\Delta_{rq}\Bigg{]}.$
(16)

The $O_{p}$ and $V_{p}$ coefficients are elements of orbital-basis-sized
vectors consisting of 1s and 0s,

$O_{p}=\begin{cases}1\;,\;p\in\mathcal{O}\\\
0\;,\;\text{otherwise}\end{cases}$ (17)
$V_{p}=\begin{cases}1\;,\;p\in\mathcal{V}\\\
0\;,\;\text{otherwise}\end{cases}$ (18)

where $\mathcal{O}$ denotes the set of active orbitals that are occupied in
the reference wave function, and $\mathcal{V}$ denotes the set of active
vacant orbitals.

The off-diagonal terms of the expectation value expansion are responsible for
electron correlation. In the case of $\Delta\text{NO}$, excitations are
limited to low-lying near-degenerate orbitals which is aimed at capturing the
so-called “static” correlation energy,

$\sum_{IJ}c_{I}c_{J}H_{IJ}=E^{\Delta\text{NO}}_{\text{stat}}$ (19)

By employing Slater-Condon rulesSzabo and Ostlund (1996) for the evaluation of
matrix elements between Slater determinants, and considering that the active
spaces for the electron pairs are disjoint, the static correlation energy can
be simplified to give,

$E^{\Delta\text{NO}}_{\text{stat}}=\sum_{pq}\left(\zeta_{pq}-\xi_{pq}\right)L_{pq},$
(20)

where

$\zeta_{pq}=V_{p}V_{q}\sum_{r}\sqrt{\Delta_{rp}\Delta_{rq}},$ (21)

and

$\xi_{pq}=O_{p}V_{q}\sqrt{n_{p}\Delta_{pq}}+O_{q}V_{p}\sqrt{n_{q}\Delta_{qp}}.$
(22)

The time-inversion exchange integral is defined as

$L_{pq}=\int\frac{\phi^{*}_{p}(\mathbf{r}_{1})\phi^{*}_{p}(\mathbf{r}_{2})\phi_{q}(\mathbf{r}_{1})\phi_{q}(\mathbf{r}_{2})}{r_{12}}d\mathbf{r}_{1}d\mathbf{r}_{2}.$
(23)

The $\Delta\text{NO}$ wave function contains only $2n$-tuple excitations, in
which the electrons remained paired, and each electron pair is excited to its
own (restricted) active space. Therefore, $E^{\Delta\text{NO}^{\prime}}$ is an
upper bound to the more general doubly-occupied configuration interaction
(DOCI) energy.Bytautas  _et al._ (2011); Bytautas, Scuseria, and Ruedenberg
(2015)

###  II.2 The high-spin correction

As defined above, the $\Delta\text{NO}$ wave function and the subsequent 2-RDM
do not account for the correlation between active electron pairs (interpair
correlation). The high-spin correction is added to $\Delta\text{NO}$ to
account for this correlation exactly in the strong correlation limit (e.g.
multiple bond dissociation). The term is ad hoc, and therefore does not
guarantee $N$-representability below the strong correlation limit.
Nevertheless, it proved effective previouslyHollett and Loos (2019) and is
added to the $N$-representable $\Delta\text{NO}$ energy, derived in Subsection
II.1, to give the total $\Delta\text{NO}$ energy,

$E^{\Delta\text{NO}}=E^{\Delta\text{NO}^{\prime}}+E^{\Delta\text{NO}}_{\text{HSC}}.$
(24)

The high-spin correction energy is defined as

$E^{\Delta\text{NO}}_{\text{HSC}}=-\sum_{pq}\kappa_{pq}K_{pq},$ (25)

where

$\kappa_{pq}=(1-\delta_{pq})\Bigg{[}\sum_{r\neq s}\xi_{pr}\xi_{ps}\\\
+W_{p}\frac{n_{p}}{2\sqrt{2}}\sum_{r}\xi_{qr}+W_{q}\frac{n_{q}}{2\sqrt{2}}\sum_{r}\xi_{pr}\Bigg{]},$
(26)

where the last two terms have been added here for modelling open-shell
systems. The coefficient $W_{p}$ indicates if the orbital is singly-occupied,

$W_{p}=\begin{cases}1\;,\;p\in\mathcal{S}\\\
0\;,\;\text{otherwise}\end{cases}$ (27)

where $\mathcal{S}$ is the set of singly-occupied orbitals. By following the
same reasoning as the high-spin correction for the doubly-occupied active
orbitals (first term of Equation 26), a factor of $\frac{1}{2}$ for the last
two terms would be predicted using the strong-correlation limit.Hollett and
Loos (2019) However, that results in too severe of a correction and therefore,
for this first application of $\Delta\text{NO}$ to open-shell systems, it has
been reduced by a factor of $\sqrt{2}$.

###  II.3 Dynamic correlation

As mentioned above, the $\Delta\text{NO}$ wave function is designed to capture
static correlation, which means dynamic correlation must be accounted for by
some other means. In this study, the dynamic correlation energy is captured
using an on-top density functional (ODF). The dynamic correlation energy,
$E^{\Delta\text{NO}}_{\text{dyn}}$, is added to the total $\Delta\text{NO}$
energy to give the total energy including static and dynamic correlation,

$E^{\Delta\text{NO}\text{-dyn}}=E^{\Delta\text{NO}}+E^{\Delta\text{NO}}_{\text{dyn}}.$
(28)

Within the current formalism, it is necessary to remove extra correlation that
is erroneously included, or “double counted”. Therefore, the total dynamic
correlation energy is the sum of the dynamic correlation energy supplied by
the ODF and a double-counting correction,

$E^{\Delta\text{NO}}_{\text{dyn}}=E_{\text{ODF}}[\Gamma^{\Delta\text{NO}}(\mathbf{r},\mathbf{r})]+E_{\text{DC}}[\Gamma^{\Delta\text{NO}}(\mathbf{r},\mathbf{r})],$
(29)

where both are a function of the two-electron on-top density,
$\Gamma(\mathbf{r},\mathbf{r})$. The on-top density is the probability that
two-electrons may be found at the same coordinate, $\mathbf{r}$. It may be
calculated from the 2-RDM, which is given by

$\tilde{\Gamma}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime})=\frac{N(N-1)}{2}\int\Psi^{*}(\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime},\mathbf{x}_{3},\dots,\mathbf{x}_{N})\\\
\times\Psi(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{3},\dots,\mathbf{x}_{N})d\mathbf{x}_{3}\dots
d\mathbf{x}_{N},$ (30)

where $\Psi$ is the $N$-electron wave function and
$\mathbf{x}=(\mathbf{r},\omega)$ is a combination of spatial $\mathbf{r}$ and
spin $\omega$ coordinates. The spinless 2-RDM may be obtained through
integration over the spin-coordinates,

$\Gamma(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})=\int\left.\tilde{\Gamma}(\mathbf{x}_{1},\mathbf{x}_{2},\mathbf{x}_{1}^{\prime},\mathbf{x}_{2}^{\prime})\right|_{\begin{subarray}{c}\omega_{1}^{\prime}=\omega_{1}\\\
\omega_{2}^{\prime}=\omega_{2}\end{subarray}}d\omega_{1}d\omega_{2},$ (31)

and the aforementioned two-electron density is the diagonal,
$\Gamma(\mathbf{r}_{1},\mathbf{r}_{2})=\Gamma(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1},\mathbf{r}_{2})$.
The spinless 2-RDM can then be resolved into different spin-components
corresponding to the relative spin of electrons 1 and 2,

$\begin{split}\Gamma(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})&=\Gamma^{\alpha\alpha}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})+\Gamma^{\beta\beta}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})\\\
&+\Gamma^{\alpha\beta}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})+\Gamma^{\beta\alpha}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime}),\end{split}$
(32)

where $\alpha$ and $\beta$ denote spin-up and spin-down electrons,
respectively.

The 2-RDM, and its components, can also be expanded in an orbital basis,

$\Gamma(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime})=\sum_{pqrs}\Gamma_{pqrs}\phi^{*}_{p}(\mathbf{r}_{1}^{\prime})\phi^{*}_{q}(\mathbf{r}_{2}^{\prime})\phi_{r}(\mathbf{r}_{1})\phi_{s}(\mathbf{r}_{2}).$
(33)

where, in this case, the $\phi_{p}$ are the $\Delta\text{NO}$ orbitals.
Following a procedure similar to the derivation of
$E^{\Delta\text{NO}^{\prime}}$, but using second quantization, one can arrive
at an expression for the elements of the spin-resolved $\Delta\text{NO}$ 2-RDM
over the $\Delta\text{NO}$ basis,

$\Gamma^{\Delta\text{NO},\sigma\sigma}_{pqrs}=\frac{f_{pq}^{\sigma\sigma}n_{p}n_{q}+\eta_{pq}}{2}\delta_{pr}^{qs}$
(34)

and

$\Gamma^{\Delta\text{NO},\sigma\sigma^{\prime}}_{pqrs}=\frac{f_{pq}^{\sigma\sigma^{\prime}}n_{p}n_{q}+\eta_{pq}}{2}\delta_{pr}\delta_{qs}\\\
-\frac{\kappa_{pq}}{2}\delta_{ps}\delta_{qr}+\frac{\zeta_{pr}-\xi_{pr}}{2}\delta_{pq}\delta_{rs},$
(35)

where $\sigma=\alpha$ or $\beta$ (with $\sigma\neq\sigma^{\prime}$),
$\delta_{pr}^{qs}=\delta_{pr}\delta_{qs}-\delta_{ps}\delta_{qr}$ and the
$f^{\sigma\sigma^{\prime}}_{pq}$ ensures there are only interactions with
$\alpha$ electrons in singly-occupied orbitals,

$\displaystyle f^{\alpha\alpha}_{pq}$ $\displaystyle=1,$ (36a) $\displaystyle
f^{\alpha\beta}_{pq}$ $\displaystyle=1-W_{q},$ (36b) $\displaystyle
f^{\beta\alpha}_{pq}$ $\displaystyle=1-W_{p},$ (36c) $\displaystyle
f^{\beta\beta}_{pq}$ $\displaystyle=\left(1-W_{p}\right)\left(1-W_{q}\right).$
(36d)

Two ODFs are used in this study. The first is a combination of the Opposite-
spin exponential cusp (OSEC) and Fermi-hole correction (FHC) functionals,
referred to as the OF functional.Hollett and Pegoretti (2018) The second is
the Colle-Salvetti (CS) functional.Colle and Salvetti (1975, 1979) The OSEC
functional depends on the total $\alpha\beta$ on-top density, whereas the FHC
functional depends on the $\alpha\alpha$ and $\beta\beta$ on-top densities
(more specifically, their Laplacians),

$E_{\text{OF}}[\Gamma(\mathbf{r},\mathbf{r})]=E_{\text{OSEC}}[\Gamma^{\alpha\beta}(\mathbf{r},\mathbf{r})+\Gamma^{\beta\alpha}(\mathbf{r},\mathbf{r})]\\\
+E_{\text{FHC}}[\Gamma^{\alpha\alpha}(\mathbf{r},\mathbf{r})]+E_{\text{FHC}}[\Gamma^{\beta\beta}(\mathbf{r},\mathbf{r})]$
(37)

The FHC functional effectively widens the Fermi-hole, which exists in any
2-RDM that obeys antisymmetry, using an exponential-cusp ansatz for the
second-order cusp. The functional is given by

$E_{\text{FHC}}[\Gamma^{\sigma\sigma}]=\int\Bigg{\\{}\pi\left(4\delta\left[3\sqrt{\pi}\left(\delta+16\delta^{3}\right)-40\delta^{2}-1\right]e^{-\frac{1}{16\delta^{2}}}\right.\\\
\left.+\sqrt{\pi}\left[3\sqrt{\pi}\left(\delta+24\delta^{3}\right)-48\left(\delta^{2}+4\delta^{4}\right)-1\right]\left[1+\text{erf}\left(\frac{1}{4\delta}\right)\right]\right)L^{\sigma\sigma}(\mathbf{r})\Bigg{\\}}\\\
/\Bigg{\\{}3\delta^{4}\left(4\left(\delta+40\delta^{3}\right)e^{-\frac{1}{16\delta^{2}}}+\sqrt{\pi}\left[1+48\left(\delta^{2}+4\delta^{4}\right)\right]\left[1+\text{erf}\left(\frac{1}{4\delta}\right)\right]\right)\Bigg{\\}}d\mathbf{r},$
(38)

where $L^{\sigma\sigma}(\mathbf{r})$ is the Laplacian of the two-electron
density with respect to the interelectronic coordinate,
$\bm{u}=\mathbf{r}_{1}-\mathbf{r}_{2}$, at the electron-electron coalescence
point ($u=0$),

$L^{\sigma\sigma}(\mathbf{r})=\nabla_{\bm{u}}^{2}\left.\Gamma^{\sigma\sigma}(\mathbf{r}+\bm{u}/2,\mathbf{r}-\bm{u}/2)\right|_{u=0},$
(39)

where $\mathbf{r}=\frac{\mathbf{r}_{1}+\mathbf{r}_{2}}{2}$. The $\delta$
determines the correlation length,

$\delta\equiv\delta(\mathbf{r})=\kappa L^{\sigma\sigma}(\mathbf{r})^{1/8}$
(40)

where $\kappa=2.30$.111This value differs from the value determined in
reference 81, $\kappa=1.73$, because it was found to overestimate the
correlation energy between electrons of parallel-spin in these H clusters.

The OSEC functional approximates the input $\alpha\beta$ two-electron density
as constant in the neighbourhood of the electron-electron coalescence point,
and determines the energy associated with introducing an exponential cusp. The
OSEC energy is given by the functional

$E_{\text{OSEC}}[\Gamma^{\alpha\beta}(\mathbf{r},\mathbf{r})]=\int\frac{2\pi\left(2\lambda\left(\sqrt{\pi}\lambda-1\right)e^{-\frac{1}{4\lambda^{2}}}+\sqrt{\pi}\left(\sqrt{\pi}\lambda-2\lambda^{2}-1\right)\left[1+\text{erf}\left(\frac{1}{2\lambda}\right)\right]\right)\Gamma^{\alpha\beta}(\mathbf{r},\mathbf{r})}{\lambda^{2}\left(2\lambda
e^{-\frac{1}{4\lambda^{2}}}+\sqrt{\pi}\left(1+2\lambda^{2}\right)\left[1+\text{erf}\left(\frac{1}{2\lambda}\right)\right]\right)}d\mathbf{r},$
(41)

where the correlation length is determined by $\lambda$,

$\lambda\equiv\lambda(\mathbf{r})=q_{\text{OSEC}}\rho(\mathbf{r})^{1/3}$ (42)

with $q_{\text{OSEC}}=2.54$ and $\rho(\mathbf{r})$ is the one-electron
density.Hollett and Pegoretti (2018)

The Colle-Salvetti ODF depends on the total (sum of all spin components) on-
top density,

$E_{\text{CS}}\left[\Gamma\left(\bm{r},\bm{r}\right)\right]=-4a\int\frac{\Gamma\left(\bm{r},\bm{r}\right)}{\rho\left(\bm{r}\right)}\\\
\left(\frac{1+b\rho\left(\bm{r}\right)^{-8/3}e^{-c\rho\left(\bm{r}\right)^{-1/3}}L(\mathbf{r})}{1+d\rho\left(\bm{r}\right)^{-1/3}}\right)d\bm{r}$
(43)

where $L(\mathbf{r})$ is the Laplacian of the total on-top density at $u=0$,
$a=0.049$, $b=0.132$, $c=0.2533$ and $d=0.349$.Colle and Salvetti (1975)

###  II.4 Double-counting correction

The ODFs, OF and CS, are both derived from correcting the two-electron density
in the neighbourhood of the electron-electron coalescence point, or in the
short-range of interelectronic distance. Multideterminantal treatments of
electron correlation, such as $\Delta\text{NO}$, capture the full range of
electron correlation. Therefore, there is overlap, or double-counting, that
can occur when wave function methods are combined with density functional
methods in the present manner. The main motivation for specifically using ODFs
with $\Delta\text{NO}$, is that as electrons become strongly correlated this
information is revealed to the ODF via $\Gamma^{\Delta\text{NO}}$. However,
the ODF only depends on the OD,
$\Gamma^{\Delta\text{NO}}(\mathbf{r},\mathbf{r})$, which does not possess
information regarding $\Delta\text{NO}$ short-range correlation beyond the
electron-electron coalescence point (in the neighbourhood of $u=0$). Because
the ODF is not sensitive to the effects of correlation at small $u\neq 0$, a
correction must be applied.

The general scheme behind the derivation of the OF functionalHollett and
Pegoretti (2018) (similar to that of the CS functionalColle and Salvetti
(1975)), is the approximation of the exact two-electron density
$\Gamma\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)$ by the uncorrelated two-
electron density $\Gamma^{0}\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)$ with a
cusp, $c(u)$, inserted at $u=0$,

$\Gamma\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)=\Gamma^{0}\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)\left(1-g(u)\left[1-\Phi(\mathbf{r})c(u)\right]\right)$
(44)

where $\Phi(\mathbf{r})$ is a function that maintains the normalization of
$\Gamma\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)$. The correction is applied
in the short-range (i.e. small $u$) and the range is controlled by the
function $g(u)$,

$g(u)=e^{-\lambda(\mathbf{r})u^{2}},$ (45)

which depends on the correlation length $\lambda(\mathbf{r})$. The derivation
is completed by Taylor-expanding
$\Gamma^{0}\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)$ about $u=0$. The
expansion is zeroth-order for the OSEC functional,

$\Gamma^{0}\left(\mathbf{r}_{1},\mathbf{r}_{2}\right)\approx\Gamma^{0}\left(\mathbf{r},\mathbf{r}\right),$
(46)

and second-order for the FHC functional. The CS functional also uses a second-
order expansion, but multiple approximations are made to ensure the integrand
is well-behaved.

If the ODF is to exclusively account for the correlation energy in the short-
range, then the short-range correlation already present in $\Gamma$ must be
removed,

$E^{\Delta\text{NO}}_{\text{dyn}}=E^{\Delta\text{NO}}_{\text{ODF}}[\Gamma]\\\
-\frac{1}{2}\int\frac{\left[\Gamma(\mathbf{r}_{1},\mathbf{r}_{2})-\Gamma^{0}(\mathbf{r}_{1},\mathbf{r}_{2})\right]g(u)}{u}d\mathbf{r}_{1}d\mathbf{r}_{2}.$
(47)

To arrive at a practical expression for the double-counting correction,
$\Gamma(\mathbf{r}_{1},\mathbf{r}_{2})$ and
$\Gamma^{0}(\mathbf{r}_{1},\mathbf{r}_{2})$, are treated with a zeroth-order
Taylor expansion about $u=0$ (i.e. constant), identical to the OSEC
functional. Also, the double-counting correction is not required in the strong
correlation limit, when correlation causes electrons to completely avoid each
other and $\Gamma(\mathbf{r},\mathbf{r})=0$, and must be turned off.
Therefore, the correction is mutliplied by an electron-pair specific factor,
$\frac{\Gamma_{pp}}{\Gamma^{0}_{pp}}$, and the resulting double-counting
correction is given as

$E_{\text{DC}}^{\Delta\text{NO}}\left[\Gamma\left(\mathbf{r},\mathbf{r}\right)\right]=2\pi
c_{\text{DC}}\\\
\int\sum_{p}O_{p}\left(\frac{\Gamma_{pp}}{\Gamma^{0}_{pp}}\right)\frac{\Gamma_{pp}^{0}\left(\mathbf{r},\mathbf{r}\right)-\Gamma_{pp}\left(\mathbf{r},\mathbf{r}\right)}{\lambda(\mathbf{r})^{2}}d\mathbf{r}$
(48)

where $\Gamma^{0}$ is the uncorrelated two-electron density, which for
$\Delta\text{NO}$ is given by

$\Gamma^{\Delta\text{NO},0}=\Gamma^{\Delta\text{NO}}_{\text{0-1RDM}}+\Gamma^{\Delta\text{NO}}_{\text{pair}}.$
(49)

The two-electron densities are decomposed into electron-pair components, which
for $u=0$, are given by

$\Gamma_{pp}^{0}\left(\bm{r},\bm{r}\right)=n_{p}\left|\phi_{p}(\mathbf{r})\right|^{4}+\sum_{q}\Delta_{pq}\left|\phi_{q}(\mathbf{r})\right|^{4},$
(50)

and

$\Gamma_{pp}\left(\bm{r},\bm{r}\right)=\Gamma_{pp}^{0}\left(\bm{r},\bm{r}\right)+\sum_{qs}\Delta_{pq}\Delta_{ps}\left|\phi_{q}(\mathbf{r})\right|^{2}\left|\phi_{s}(\mathbf{r})\right|^{2}\\\
-2\sum_{q}\xi_{pq}\left|\phi_{p}(\mathbf{r})\right|^{2}\left|\phi_{q}(\mathbf{r})\right|^{2}.$
(51)

The above correction is added to the OF functional energy to give the total
dynamic correlation energy. In the case of the CS functional,

$\lambda(\mathbf{r})=q_{\text{CS}}\rho(\mathbf{r})^{1/3}$ (52)

where $q_{\text{CS}}=2.29$ (original value from CS derivation).

###  II.5 2-RDM optimization

The $\Delta\text{NO}$ 2-RDM depends upon the $\Delta\text{NO}$ orbitals,
$\\{\phi_{p}\\}$, and the electron transfer variables, $\\{\Delta_{ia}\\}$.
The optimal $\Delta\text{NO}$ 2-RDM is found by minimizing
$E^{\Delta\text{NO}}$ with respect to both $\\{\phi_{p}\\}$ and
$\\{\Delta_{ia}\\}$. The optimization can be performed for both
$\\{\phi_{p}\\}$ and $\\{\Delta_{ia}\\}$ simultaneously, using the trust-
region method.

In order to obtain the optimal $\Delta\text{NO}$ orbitals, $\\{\phi_{p}\\}$,
the initial set of orthonormalized orbitals, $\\{\tilde{\phi}_{r}\\}$, are
rotated amongst themselves via a unitary transformation,

$\phi_{p}=\sum_{r}U_{rp}\tilde{\phi}_{r}.$ (53)

To determine such a transformation, the unitary matrix is parameterized with
the exponential of a skew matrix,Helgaker, Jorgensen, and Olsen (2014)
$\mathbf{Y}$,

$\mathbf{U}=e^{\mathbf{Y}},$ (54)

where

$Y_{pq}=y_{pq},\quad Y_{qp}=-y_{pq}.$ (55)

Any energy, $E[\mathbf{y}]$, that is a functional of a set of orbitals, may be
expanded about the optimal orbitals, in a Taylor series, in terms of the skew
matrix parameters, $\mathbf{y}$,

$E[\mathbf{y}]\approx
E[\mathbf{0}]+\mathbf{y}\cdot\mathbf{g}^{\phi}+\frac{1}{2}\mathbf{y}^{\dagger}\mathbf{H}^{\phi\phi}\mathbf{y},$
(56)

and approximated to second-order.Bozkaya  _et al._ (2011); Bozkaya and
Sherrill (2013, 2014); Bozkaya (2016) The gradient, $\mathbf{g}^{\phi}$, and
hessian, $\mathbf{H}^{\phi}$, consist of the first and second derivatives of
the energy with respect to the skew-matrix parameters,

$g^{\phi}_{pq}=\left.\frac{\partial E^{\Delta\text{NO}}}{\partial
y_{pq}}\right|_{\mathbf{y}=0},\quad
H^{\phi\phi}_{pq,rs}=\left.\frac{\partial^{2}E^{\Delta\text{NO}}}{\partial
y_{pq}\partial y_{rs}}\right|_{\mathbf{y}=0}.$ (57)

General expressions for the above derivatives may be found elsewhere,Bozkaya
_et al._ (2011) but for $\Delta\text{NO}$ specifically see Appendix A. In the
case of $E^{\Delta\text{NO}}$, the energy is invariant ($g^{\phi}_{pq}=0$)
with respect to rotations between orbital pairs where both belong to the
inactive (frozen) closed-shell, open-shell, or inactive vacant orbital
subspaces. Consequently, optimization must be performed with respect to all
other possible rotations.

In the case of $\Delta\text{NO}$, the energy is also a function of the
electron transfer variables, $\\{\Delta_{ia}\\}$. To simultaneously optimize
the energy with respect to both $\\{\phi_{p}\\}$ and $\\{\Delta_{ia}\\}$, the
orbital gradient and hessian are combined with those for $\\{\Delta_{ia}\\}$.
However, in order to avoid violation of $N$-representability of the 1-RDM, the
$\\{\Delta_{ia}\\}$ are first parameterized such that they remain between 0
and 1/2,

$\Delta_{ia}=\frac{\cos^{2}\theta_{ia}}{2}.$ (58)

The total gradient and hessian are then given by

$\mathbf{g}=\begin{bmatrix}\mathbf{g}^{\phi}\\\
\mathbf{g}^{\theta}\end{bmatrix},\quad\mathbf{H}=\begin{bmatrix}\mathbf{H}^{\phi\phi}&\mathbf{H}^{\phi\theta}\\\
\mathbf{H}^{\phi\theta}&\mathbf{H}^{\theta\theta}\end{bmatrix},$ (59)

where

$g^{\theta}_{pq}=\frac{\partial
E^{\Delta\text{NO}}}{\partial\theta_{pq}},\quad
H^{\theta\theta}_{pq,rs}=\frac{\partial^{2}E^{\Delta\text{NO}}}{\partial\theta_{pq}\partial\theta_{rs}},$
(60)

and

$H^{\phi\theta}_{pq,rs}=\frac{\partial^{2}E^{\Delta\text{NO}}}{\partial
y_{pq}\partial\theta_{rs}}.$ (61)

Minimization can be performed using the trust-region Newton method, in which
an augmented hessian is prepared and a step is taken within a dynamic trust-
region radius. A standard trust-region approach is described in Nocedal and
Wright (Algorithm 4.1).Nocedal and Wright (2006) The procedure begins with the
calculation of an augmented hessian, $\mathbf{A}$,

$\mathbf{A}=\mathbf{H}+\epsilon\mathbf{I}$ (62)

where $\mathbf{I}$ is the identity matrix and $\epsilon$ is chosen such that
$\mathbf{A}$ is positive definite. Then a Newton step, $\mathbf{p}$, is
calculated by solving

$\mathbf{A}\mathbf{p}=-\mathbf{g}.$ (63)

If the step is within the trust-region radius,
$p=\sqrt{\mathbf{p}\cdot\mathbf{p}}<t$, then the step is taken, otherwise a
subproblem must be solved (Equation 4.5 in reference 90) in which a step
within the trust-region radius is found iteratively. Each step, whether pure
Newton-Raphson or a solution to the subproblem, is evaluated by comparing the
actually change in the function ($E$) to the change predicted by a second-
order Taylor expansion. If the reduction is less than a minimum threshold
(e.g. 25% of predicted) then no step is taken and $t$ is reduced. If the
reduction is more than the minimum, a step is taken, and if it is more than a
maximum threshold (e.g. 75% of predicted) then $t$ is increased.

In order to avoid optimizations to saddle points in $\\{\phi_{p}\\}$ and
$\\{\Delta_{ia}\\}$ space, a modification can be made before commencing the
trust-region algorithm. The non-augmented hessian is diagonalized to evaluate
the eigenvalues,

$\mathbf{Z}^{T}\mathbf{H}\mathbf{Z}=\mathbf{H}^{\prime}$ (64)

where $\mathbf{H}^{\prime}$ is a diagonal matrix of hessian eigenvalues and
$\mathbf{Z}$ is the matrix that transforms to normal coordinates. For any
normal coordinate (mode) with a negative eigenvalue, $H_{ii}^{\prime}<0$, and
a gradient component with a magnitude less than a threshold,
$\left|g_{i}^{\prime}\right|<g_{0}$, the gradient is amplified to promote
escape from the saddle point,

$\mathbf{g}_{\text{mod}}=\mathbf{Z}\mathbf{g}_{\text{mod}}^{\prime}$ (65)

where

$g_{\text{mod},i}^{\prime}=\begin{cases}\text{sign}(g_{i}^{\prime})g_{0},&\text{
If }H_{ii}^{\prime}<0\text{ and }\left|g_{i}^{\prime}\right|<g_{0}\\\
g_{i}^{\prime},&\text{ otherwise}\end{cases}$ (66)

The trust-region method then continues as usual.

When an optimization step is calculated, the unitary matrix for the orbital
rotation is prepared by first approximating $e^{\mathbf{Y}}$ to second order,

$\mathbf{U}\approx\mathbf{I}+\mathbf{Y}+\frac{1}{2}\mathbf{Y}\mathbf{Y}^{{\dagger}},$
(67)

and then ensuring it is unitary by applying Gram-Schmidt orthonormalization.

###  II.6 Conditions of $N$-representability

Due to the addition of the high-spin correction, the $\Delta\text{NO}$ 2-RDM
is not strictly $N$-representable. Therefore, it is useful to have some means
of evaluating the extent to which $N$-representability might be violated. The
necessary, but not sufficient, $N$-representability conditions that depend
only on the 2-RDM are referred to as the $PQG$ (or $DQG$) conditions.Garrod
and Percus (1964); Mazziotti (2006, 2012) They require that the following
matrices, in a spin-orbital basis and defined in terms of the second
quantization, be positive semi-definite (i.e. all non-negative eigenvalues)

$\displaystyle P_{ijkl}$
$\displaystyle=\langle\Psi|a_{j}^{\dagger}a_{i}^{\dagger}a_{k}a_{l}|\Psi\rangle,$
(68a) $\displaystyle Q_{ijkl}$
$\displaystyle=\langle\Psi|a_{j}a_{i}a_{k}^{\dagger}a_{l}^{\dagger}|\Psi\rangle,$
(68b) $\displaystyle G_{ijkl}$
$\displaystyle=\langle\Psi|a_{j}^{\dagger}a_{i}a_{k}^{\dagger}a_{l}|\Psi\rangle,$
(68c)

where $\Psi$ is the $N$-electron wave function, $a_{i}^{\dagger}$ and $a_{i}$
are creation and annihilation operators, respectively. Sufficient conditions
for $N$-representability require the calculation of higher-order matrices and
are not considered here.Mazziotti (2012); Li, Liebenthal, and DePrince (2021)
For a spin-resolved 2-RDM over spatial orbitals [(LABEL:eq:2RDMps) and (35)],
and a diagonal 1-RDM, the elements of $\bm{P}$, $\bm{Q}$, and $\bm{G}$ can be
expressed as

$\displaystyle P^{\sigma\sigma^{\prime}\sigma\sigma^{\prime}}_{pqrs}=$
$\displaystyle\;\Gamma^{\sigma\sigma^{\prime}}_{pqrs},$ (69a) $\displaystyle
Q^{\sigma\sigma^{\prime}\sigma\sigma^{\prime}}_{pqrs}=$
$\displaystyle\;\Gamma^{\sigma\sigma^{\prime}}_{pqrs}+\frac{\delta_{pr}\delta_{qs}}{2}\left(1-n_{p}-n_{q}\right)$
$\displaystyle-\delta_{\sigma\sigma^{\prime}}\frac{\delta_{ps}\delta_{qr}}{2}\left(1-n_{p}-n_{q}\right),$
(69b) $\displaystyle G^{\sigma\sigma^{\prime}\sigma\sigma^{\prime}}_{pqrs}=$
$\displaystyle\;\frac{\delta_{pr}\delta_{qs}}{2}n_{q}-\Gamma^{\sigma\sigma^{\prime}}_{rqps},$
(69c) $\displaystyle G^{\sigma\sigma\sigma^{\prime}\sigma^{\prime}}_{pqrs}=$
$\displaystyle\;\Gamma^{\sigma\sigma^{\prime}}_{qrps}\quad(\sigma\neq\sigma^{\prime}),$
(69d)

where $\sigma=\alpha$ or $\beta$, and $\sigma$ can be equivalent to
$\sigma^{\prime}$ (except for
$G^{\sigma\sigma\sigma^{\prime}\sigma^{\prime}}_{pqrs}$). The expression of
$\bm{P}$, $\bm{Q}$, and $\bm{G}$ in terms of spatial-orbital indices
(subscript) requires specifying the spin-function (superscript) used to form
the corresponding spin orbital. The matrices can be decomposed into spin-
blocks for construction and analysis, which is evident from their
structure,Mazziotti (2002)

$\displaystyle\bm{P}$
$\displaystyle=\begin{bmatrix}\bm{P}^{\alpha\alpha\alpha\alpha}&0&0&0\\\
0&\bm{P}^{\beta\beta\beta\beta}&0&0\\\
0&0&\bm{P}^{\alpha\beta\alpha\beta}&0\\\
0&0&0&\bm{P}^{\beta\alpha\beta\alpha}\end{bmatrix},$ (70a)
$\displaystyle\bm{Q}$
$\displaystyle=\begin{bmatrix}\bm{Q}^{\alpha\alpha\alpha\alpha}&0&0&0\\\
0&\bm{Q}^{\beta\beta\beta\beta}&0&0\\\
0&0&\bm{Q}^{\alpha\beta\alpha\beta}&0\\\
0&0&0&\bm{Q}^{\beta\alpha\beta\alpha}\end{bmatrix},$ (70b)
$\displaystyle\bm{G}$
$\displaystyle=\begin{bmatrix}\bm{G}^{\alpha\alpha\alpha\alpha}&\bm{G}^{\alpha\alpha\beta\beta}&0&0\\\
\bm{G}^{\beta\beta\alpha\alpha}&\bm{G}^{\beta\beta\beta\beta}&0&0\\\
0&0&\bm{G}^{\alpha\beta\alpha\beta}&0\\\
0&0&0&\bm{G}^{\beta\alpha\beta\alpha}\end{bmatrix}.$ (70c)

The sum of the negative eigenvalues of these matrices can be used to quantify
the extent to which they are not positive semi-definite and subsequently
violate the $PQG$-conditions.Rodríguez-Mayorga  _et al._ (2017) These sums are
denoted by

$X_{\text{neg}}=\sum_{\begin{subarray}{c}i=1\\\
x_{i}<0\end{subarray}}^{K^{2}}x_{i},$ (71)

where $K$ is the size of the orbital basis, $\bm{X}=\bm{P},\bm{Q}$ or $\bm{G}$
and $x_{i}$ is an eigenvalue of $\bm{X}$.

##  III Method

###  III.1 Reference energies and properties

The energies from single-reference methods: completely renormalized coupled-
cluster (CR-CCL)Piecuch  _et al._ (2002); Piecuch and Włoch (2005) and spin-
unrestricted density functional approximation (B3LYP),Becke (1993); Stephens
_et al._ (1994); Hertwig and Koch (1997) and multireference energies:
multireference self-consistent field (MRSCF)Chaban, Schmidt, and Gordon (1997)
and multireference Möller-Plesset second-order perturbation theory
(MRMP2),Hirao (1992b, c, d, 1993) were calculated using GAMESS.Barca  _et al._
(2020) The active space for the MRSCF and MRMP2 calculations was chosen to be
all occupied orbitals of the reference configuration, and an equivalent number
of the lowest lying vacant (virtual) orbitals. Full configuration interaction
(FCI) energies, and 2-RDMs, were determined using Quantum Package 2.Garniron
_et al._ (2019) All calculations were performed using the cc-pVTZ basis
set.Pritchard  _et al._ (2019); Dunning Jr. (1989)

###  III.2 $\Delta\text{NO}$ energies and properties

The $\Delta\text{NO}$, $\Delta\text{NO-OF}$, and $\Delta\text{NO-CS}$
calculations were performed with MUNgauss.Poirier, Hollett, and Warburton
(2015) The active space for each electron pair consisted of the orbital they
occupy in the reference configuration and one vacant (virtual) orbital. The
$\Delta\text{NO}$ 2-RDMs were optimized using the the trust-region method
outlined in Subsection II.5. The iterative subproblem was solved using the
dgqt routine from MINPACK-2,Moré and Sorensen (1983) translated to Fortran 90.
Restricted Hartree-Fock orbitals are used as the intial guess, and the intial
trust-region radius was $t=0.05$ with a maximum radius value of 2.00. The
gradient threshold, required for saddle-point escape, was $g_{0}=10^{-4}$. The
convergence criteria was $\bar{g}<10^{-6}$ where
$\bar{g}=\frac{|\mathbf{g}|}{n}$ is the normalized total gradient length,
where $n$ is the total number of orbital optimization parameters,
$\\{y_{pq}\\}$ and $\\{\Delta_{ia}\\}$.

The OF and CS energies are calculated post-$\Delta\text{NO}$ using an SG-1
gridGill, Johnson, and Pople (1993) and Becke atomic weights.Becke (1988)

The value of the double-counting correction parameter, $c_{\text{DC}}$, was
chosen to minimize the error in the $\Delta\text{NO-OF}$, and $\Delta\text{NO-
CS}$, H2 potential energy curves, relative to FCI. The values were 0.40 and
0.35 for $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$, respectively.

##  IV Results

###  IV.1 Cluster geometries

The geometries of the H clusters studied are presented in Figure 1.

Figure 1:  Geometries of H clusters

The potential energy surfaces were calculated with established single-
reference (FCI, CR-CCL and B3LYP) and multireference methods (MRSCF and MRMP2)
as well as with $\Delta\text{NO}$, $\Delta\text{NO-OF}$ and $\Delta\text{NO-
CS}$. The potential energy surface is defined as,

$U(\mathbf{X})=E(\mathbf{X})+V_{nn}(\mathbf{X})$ (72)

where $E(\mathbf{X})$ is the electronic energy provided by the electronic
structure method, $V_{nn}(\mathbf{X})$, is the nuclear-repulsion energy, and
both are a function of the nuclear coordinates, $\mathbf{X}$.

###  IV.2 2-RDM optimization

The most time consuming part of the $\Delta\text{NO}$ 2-RDM optimization is
the transformation of the two-electron integrals from the atomic orbital basis
to the $\Delta\text{NO}$ orbital basis, which scales as $N_{\text{act}}K^{4}$,
where $N_{\text{act}}$ is the number of active orbitals and $K$ is the number
of atomic basis functions. Recently, this cost has been reduced for similar
calculations in natural orbital functional methods by introducing the
resolution-of-the-identity (DoNOF-RI).Lew-Yee, Piris, and M. del Campo (2021)

Table 1:  Comparison of the number of integral transformations required to optimize the $\Delta\text{NO}$ 2-RDM |  | integral transformations
---|---|---
cluster |  $R$ (Å) | ID222ID = iterative diagonalization of pseudo-Fock matrix | TR333TR = trust-region method
H2 | 1.2 | 170 | 10
H3 | 1.2 | 924 | 13
H4 linear | 1.2 | 252 | 47
H4 square | 1.2 | 320 | 20

For the systems studied, the use of the trust-region method, on average,
reduced the number of optimization steps (and therefore the number of integral
transformations) by an order of magnitude compared to the previous iterative
diagonalization algorithmPiris and Ugalde (2009); Hollett, Hosseini, and
Menzies (2016); Hollett and Loos (2019) (Table 1). Furthermore, the
introduction of an analytical hessian, and a saddle-point escape procedure,
enabled the trust-region method to find global minima that the previous
iterative-diagonalization algorithm could not locate.

In some cases, in order to determine whether optimized 2-RDMs were indeed
global minima, multiple different initial guesses were used. This includes
using orbitals different from RHF, including MRSCF orbitals and FCI natural
orbitals, choosing different vacant orbitals for the active space, and also
using 2-RDMs (both the orbitals and the electron-transfer variables) from
other cluster geometries. The minima found by the 2-RDM optimization algorithm
were invariant with respect to the choice of initial orbitals or which vacant
orbitals to include in the active space. Local minima were only located when
an alternative 2-RDM guess was supplied from another cluster geometry (e.g.
2-RDM from H4 at $\theta=90^{\circ}$ for H4 at $\theta=85^{\circ}$).

###  IV.3 H2

The error in the H2 potential energy curves relative to the FCI result is
presented in Figure 2.

Figure 2:  Error in H2 potential energy curves compared to FCI ($\Delta
U=U-U_{\text{FCI}}$).

To illustrate the importance of the double-counting correction, the
$\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ result without including the
double-counting correction ($c_{\text{DC}}=0$) is shown. Without the double-
counting correction, both $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$
predict $\sim$ 30 kJ mol-1 too much correlation energy near the equilibirum
bond length ($R=0.743$ Å). However, as the amount a static correlation grows
with increasing $R$, the error approaches zero. The static correlation
included in the $\Delta\text{NO}$ 2-RDM causes the on-top density to vanish,
and therefore there is no superfluous correlation energy at dissociation. If
the double-counting correction is included, the error in $\Delta\text{NO-OF}$
and $\Delta\text{NO-CS}$ is dramatically reduced with a maximum error of
$\sim$ 3 kJ mol-1 (a positive deviation near equilibrium, and negative at
stretched bond lengths). The error in the B3LYP potential energy curve is also
included, where the most notable deviation is the energy at the dissociation
limit. The dissociation limit predicted by B3LYP is 5.6 kJ mol-1 higher than
FCI.

The equilibrium bond lengths and dissociation energies predicted by each
method for the various H clusters are presented in Tables 2 and 3. In the case
of H2, $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ predict equilibrium bond
lengths within 0.004 Å of the FCI value, 0.0743 Å, while the other methods are
slightly more accurate. Also for H2, $\Delta\text{NO}$ and MRSCF are identical
and therefore predict that same bond length (0.756 Å) and dissociation energy
(399 kJ mol-1).

Table 2:  Equilibrium bond lengths (in Å) of H clusters | cluster
---|---
method |  H2 |  H3 |  H4 linear |  H4 square
$\Delta\text{NO}$ | 0.756 | 0.942 | 0.908 | 0.910
$\Delta\text{NO-OF}$ | 0.739 | 0.924 | 0.882 | 0.863
$\Delta\text{NO-CS}$ | 0.739 | 0.918 | 0.879 | 0.857
B3LYP | 0.744 | 0.931 | 0.885 | 0.841
CR-CCL | 0.743 | 0.930 | 0.887 | 0.878
MRSCF | 0.756 | 0.957 | 0.906 | 0.912
MRMP2 | 0.745 | 0.934 | 0.889 | 0.868
FCI | 0.743 | 0.930 | 0.887 | 0.864

Other than $\Delta\text{NO}$ and MRSCF, all methods predict fairly accurate
bond dissociation energies for H2. The high accuracy of $\Delta\text{NO-OF}$
(0 kJ mol-1 error) and $\Delta\text{NO-CS}$ (1 kJ mol-1 error) is due to the
calibration of the double-counting correction coefficient, $c_{\text{DC}}$.
The strategy is to use the bond of H2 as a prototype for the static and
dynamic correlation energy balance between paired electrons.

Table 3:  Dissociation energies (in kJ mol-1) of H clusters | cluster
---|---
method |  H2 |  H3 |  H4 linear |  H4 square
$\Delta\text{NO}$ | 399 | 317 | 588 | 168
$\Delta\text{NO-OF}$ | 454 | 411 | 722 | 288
$\Delta\text{NO-CS}$ | 455 | 413 | 718 | 288
B3LYP | 461 | 444 | 765 | 349
CR-CCL | 454 | 412 | 726 | 292
MRSCF | 399 | 327 | 610 | 202
MRMP2 | 444 | 400 | 713 | 278
FCI | 454 | 412 | 733 | 288

###  IV.4 H3

The potential energy curves (PECs) for H3, in a linear arrangement, are
presented in Figure 3.

Figure 3:  Potential energy of a linear H3 cluster with neighbouring H atoms
separated by a distance $R$.

The $\Delta\text{NO}$ and MRSCF PECs are above the PECs of the other methods
designed to also capture dynamic correlation. Unlike H2, the $\Delta\text{NO}$
and MRSCF PECs differ and it is apparent that, even with the high-spin
correction applied to the unpaired electron (Equation 25), $\Delta\text{NO}$
is missing correlation at moderately stretched bond lengths. Nevertheless, the
correct dissociation limit is reached by $\Delta\text{NO}$, $\Delta\text{NO-
OF}$ and $\Delta\text{NO-CS}$.

The error in the PECs, for the methods that include both static and dynamic
correlation ($\Delta\text{NO-OF}$, $\Delta\text{NO-CS}$, MRMP2 and B3LYP),
compared to the FCI PEC is shown in Figure 4.

Figure 4:  Error in H3 potential energy curves compared to FCI ($\Delta
U=U-U_{\text{FCI}}$).

At small $R$, $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ overestimate the
correlation energy by a maximum of 4 and 11 kJ mol-1, respectively. At
moderately stretched bond lengths, both $\Delta\text{NO-OF}$ and
$\Delta\text{NO-CS}$ underestimate the correlation energy (23 and 29 kJ mol-1,
respectively), which is due to the missing static correlation in
$\Delta\text{NO}$ between the electron pair and the unpaired electron.

The largest deviation of the MRMP2 PEC from the FCI PEC is at small $R$ (16 kJ
mol-1), and the error decreases monotonically as $R$ increases. In the case of
B3LYP, the correlation energy is significantly overestimated around
equilibrium $R$ (24 kJ mol-1) and like H2, the energy at the dissociation
limit is incorrect. These errors are reflected in the predicted equilibrium
$R$ values and dissociation energies (Tables 2 and 3).

Similar to H2, both $\Delta\text{NO}$ and MRSCF overestimate the equilibrium
separation. Whereas, $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ slightly
underestimate $R_{e}$, by 0.006 Å and 0.012 Å, respectively. The values of
$R_{e}$ from both MRMP2 (0.934 Å) and B3LYP (0.931 Å) differ only slightly
from the FCI value. Also, like H2, the CR-CCL PEC is equivalent to FCI.

For those methods designed to capture all the correlation energy, it is only
the dissociation energy predicted by B3LYP (444 kJ mol-1) that differs
substantially from the FCI value (412 kJ mol-1).

###  IV.5 Linear H4

The PECs for the linear H4 cluster, with interatomic spacing $R$, are
presented in Figure 5.

Figure 5:  Potential energy of a linear H4 cluster with neighbouring H atoms
separated by a distance $R$.

Similar to H3, the $\Delta\text{NO}$ and MRSCF PECs appear above the others
(as expected), and the B3LYP PEC lies below the others at equilibrium and
above at dissociation. Interestingly, while the CR-CCL PEC is essentially
exact close to $R_{e}$, and eventually approaches the correct dissociation
limit, there is a significant negative deviation from the FCI PEC at large
$R$.

The error in the linear H4 PECs, compared to FCI, are explicitly shown in
Figure 6.

Figure 6:  Error in linear H4 cluster potential energy curves compared to FCI
($\Delta U=U-U_{\text{FCI}}$).

The behaviour of the error in the $\Delta\text{NO-OF}$ and $\Delta\text{NO-
CS}$ PECs for linear H4 is the same as that seen for H2, although
substantially larger. Both $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$
underestimate the correlation energy at small $R$, and overestimate at large
$R$, while eventually reaching the correct dissociation limit. Like H3, the
error in the MRMP2 PEC is largest at small $R$ and decreases as $R$ is
increased. In the case of CR-CCL, the PEC is essentially exact for small $R$
and then begins overestimating the correlation energy at $\sim$ 1.4 Å. The
error in the B3LYP PEC is similar that seen for both H2 and H3.

With the exception of $\Delta\text{NO}$ and MRSCF, all methods provide good
estimates of $R_{e}$ compared to FCI, with $\Delta\text{NO-OF}$ and
$\Delta\text{NO-CS}$ showing the most deviation (0.005 and 0.008 Å,
respectively). Both $\Delta\text{NO}$ and MRSCF predict larger $R_{e}$ values,
and are within close agreement with each other (differing by 0.002 Å).

The dissociation energies predicted by each method are more varied. The
dissociation energies predicted by $\Delta\text{NO-OF}$ and $\Delta\text{NO-
CS}$, 722 kJ mol-1 and 718 kJ mol-1, respectively, are in decent agreement
with the FCI value, 733 kJ mol-1. Only CR-CCL provides a better estimate, 726
kJ mol-1.

###  IV.6 Square H4

The PECs for the dissociation of H4 in a square arrangement
($\theta=90^{\circ}$, see Figure 1) are shown in Figure 7.

Figure 7:  Potential energy of a square H4 cluster where each at is a distance
$R$ from the centre (see Figure 1).

The square arrangement of H atoms results in much more varied performance
amongst the methods assessed compared to the linear arrangement. The B3LYP PEC
is far below the others near $R_{e}$, but above again at the dissociation
limit. There is substantial separation between the $\Delta\text{NO}$ and MRSCF
PECs near $R_{e}$. Furthermore, this system is a known issue for single-
reference coupled-cluster methods, and while CR-CCL performs better here than
CCSD(T), the error is significant. The explicit errors are shown in Figure 8.

Figure 8:  Error in square H4 cluster potential energy curves compared to FCI
($\Delta U=U-U_{\text{FCI}}$).

The error in the $\Delta\text{NO-CS}$ PEC is relatively small, and provides,
arguably, the best agreement with the FCI PEC. The usual underestimation of
the correlation energy at small $R$ by MRMP2 is relatively small for this
system, and therefore MRMP2 also provides a good estimate of the FCI PES. The
negative deviation of $\Delta\text{NO-OF}$ at moderate $R$ values reaches a
maximum of $\sim$ 14 kJ mol-1 before going to the correct dissociation limit.
Whereas, the CR-CCL PES overestimates the correlation energy by more than 30
kJ mol-1 before approaching the dissociation limit.

Despite the errors in the PECs, accurate predictions of the FCI $R_{e}$, 0.864
Å, are provided by $\Delta\text{NO-OF}$, $\Delta\text{NO-CS}$, MRMP2 and CR-
CCL, with values of 0.863, 0.857, 0.868 and 0.878, respectively. As usual, the
dissociation energies vary more widely. Both $\Delta\text{NO-OF}$ and
$\Delta\text{NO-CS}$ provide very accurate estimates of the dissociation
energy, and MRMP2 and CR-CCL differ by only 10 and 4 kJ mol-1, respectively.
While B3LYP overestimates the dissociation energy by 61 kJ mol-1.

###  IV.7 Rectangle to square H4

At fixed $R$, if $\theta$ of the, $D_{4h}$ symmetry, square arrangement of H4
(Figure 1) is changed from 90∘, the result is a rectangular, $D_{2h}$
symmetry, arrangement. The PECs for $\theta=70^{\circ}$ to
$\theta=110^{\circ}$ ($D_{2h}$ $\to$ $D_{4h}$ $\to$ $D_{2h}$) for fixed $R$
values of 0.8, 1.2 and 1.7 Å are presented in Figure 9.

Figure 9:  Potential energy of H4 cluster as a function of $\theta$ and $R$
(see Figure 1).

Like the PECs for the other H clusters, the $\Delta\text{NO}$ and MRSCF PECs
lie above those of the methods designed to capture all of the correlation
energy. As $R$ increases, and static correlation increases while dynamic
correlation decreases, the $\Delta\text{NO}$ and MRSCF PECs approach the
others (note the scale of the $y$-axis in Figure 9).

In the case of CR-CCL, coupled-cluster methods are known to have a concave
(downward pointing) derivative discontinuity,Robinson and Knowles (2012c);
Ramos-Cordoba  _et al._ (2015); Marie, Kossoski, and Loos (2021) or cusp, at
$\theta=90$, which becomes more prominent as $R$ increases. This is in
contrast to the convex (upward pointing) cusp seen for RHF and spin-restricted
density functional methods. The cusp in the CR-CCL PEC is barely noticable at
$R=0.8$ Å, but is evident at the larger $R$ values. The MRMP2 PEC closely
resembles the FCI PEC at all $R$ values, and the absolute error decreases with
increasing $R$. The height of the PEC is significantly underestimated by B3LYP
at all $R$ values.

The $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ methods accurately model the
PEC near $\theta=70^{\circ}$ and 90∘, however there is an obvious cusp at
$\theta\sim 85^{\circ}$ and $\theta\sim 95^{\circ}$. This is due to an
inadequate description of the 2-RDM upon transitioning from a $D_{4h}$ to
$D_{2h}$ arrangement, which is analysed, and elaborated upon, later.

The barrier height, defined as the difference in energy at $\theta=90^{\circ}$
and 70∘, predicted by each method at each $R$ value is presented in Table 4.

Table 4:  Barrier height555Barrier height is defined as $E[90^{\circ}]-E[70^{\circ}]$ (in kJ mol-1) for $D_{2h}\to D_{4h}\to D_{2h}$ transition of H4 cluster at varying $R$. |  $R$ (Å)
---|---
method | 0.8 | 1.2 | 1.7
$\Delta\text{NO}$ | 305 | 212 | 76
$\Delta\text{NO-OF}$ | 287 | 200 | 76
$\Delta\text{NO-CS}$ | 284 | 200 | 78
B3LYP | 238 | 157 | 50
CR-CCL | 292 | 186 | 51
MRSCF | 281 | 190 | 72
MRMP2 | 286 | 204 | 80
FCI | 293 | 208 | 81

At $R=0.8$ Å, $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$ underestimate the
barrier by 6 and 9 kJ mol-1, respectively. Compare that to $\Delta\text{NO}$,
which overestimates the barrier by 12 kJ mol-1. MRMP2 provides accuracy
similar to $\Delta\text{NO-OF}$ and $\Delta\text{NO-CS}$, while CR-CCL only
differs from the FCI barrier by 1 kJ mol-1. As $R$ is increased, the agreement
of $\Delta\text{NO}$ with FCI improves, which highlights the importance of
static correlation to these PECs. At $R=1.2$, $\Delta\text{NO-OF}$ and
$\Delta\text{NO-CS}$ both differ from FCI by 8 kJ mol-1, and differ by 5 and 3
kJ mol-1, respectively, at $R=1.7$ Å. At the same time, the agreement of MRMP2
with FCI also improves, to within 4 kJ mol-1 at $R=1.2$ Å and within 1 kJ
mol-1 at $R=1.7$ Å. In contrast, due to the cusp, the CR-CCL values worsen.

Although the barrier heights predicted by $\Delta\text{NO-OF}$ and
$\Delta\text{NO-CS}$ are in good agreement with FCI values, the deviation of
the PECs from FCI approximately 5∘ from $\theta=90^{\circ}$ is concerning. An
initial analysis of the correlation captured by $\Delta\text{NO}$ is provided
by the electron transfer variables, $\\{\Delta_{ia}\\}$. The values of
$\Delta_{14}$ and $\Delta_{23}$ from $\theta=70^{\circ}$ to
$\theta=110^{\circ}$, at $R=1.7$ Å, are presented in Figure 10.

Figure 10:  $\Delta$ values for rectangle-to-square H4 cluster as a function
of $\theta$ at $R=1.7$ Å.

From $\theta=70^{\circ}$ to $\theta=77^{\circ}$, $\Delta_{14}$ and
$\Delta_{23}$ are equivalent, indicating that the motion of both electron
pairs is correlated in the same manner. This points to the rectangular cluster
behaving as two H2 molecules, which is not surprising at particularly large
separations like $R=1.7$ Å. Beyond $\theta=77^{\circ}$, the values of
$\Delta_{14}$ and $\Delta_{23}$ diverge, and there is an obvious transition
that occurs from $\theta=84^{\circ}$ to $\theta=87^{\circ}$, and from
$\theta=93^{\circ}$ to $\theta=96^{\circ}$ (grey rectangles). At the
$\theta=90^{\circ}$ side of the transition, the electron pairs have
significantly different amounts of static correlation. The $\Delta_{23}$
electron pair is very near the strong correlation limit, $\Delta\to 1/2$,
while $\Delta_{14}\approx 0.28$. Essentially, there are two regimes of
correlation. For $\theta$ near 70∘ (and 110∘) there is the 2H2 regime, and for
$\theta$ near 90∘ there is the H4 regime.

The presence of two different regimes, or models, for the correlation in H4 by
$\Delta\text{NO}$, and their transition, is supported further by an analysis
of the $\Delta\text{NO}$ orbitals (Figure 11).

| $\theta$ (degrees)
---|---
orbital | 84 | 85 | 86 | 87
1 |  |  |  |
2 |  |  |  |
3 |  |  |  |
4 |  |  |  |
Figure 11:  Change in $\Delta\text{NO}$ active orbitals of H4 with $\theta$.

From Figure 11, it is seen that from $\theta=84^{\circ}$ to $85^{\circ}$ the
$\Delta\text{NO}$ orbitals resemble those of two separate H2 molecules, albeit
with small tails on the opposite pair of protons. At $\theta=86^{\circ}$, the
orbitals begin to extend significantly to the neighbouring pair of protons,
and by $\theta=87^{\circ}$ they are almost evenly distributed. This transition
occurs quickly with respect to the PEC, in just 3∘.

Analysis of the $\Delta$ values and the orbitals of $\Delta\text{NO}$
partially explains the behaviour of the model for the this challenging H4 PEC,
but it is does not explain why it differs from the FCI result. A comparison of
the $\Delta\text{NO}$ and FCI two-electron densities,
$\Gamma^{\Delta\text{NO}}(\mathbf{r}_{1},\mathbf{r}_{2})$ and
$\Gamma(\mathbf{r}_{1},\mathbf{r}_{2})$, can show conclusively how the
$\Delta\text{NO}$ model of electron correlation differs from how the electrons
actually behave. The $\Delta\text{NO}$ and FCI opposite-spin ($\alpha\beta$)
and parallel-spin ($\alpha\alpha$) two-electron densities for H4 at
$\theta=70^{\circ},85^{\circ}$ and 90∘, with $R=1.7$ Å, are shown in Figure
12.

$\Gamma^{\alpha\beta}(\mathbf{r},\mathbf{r}_{0})$

$\Gamma^{\Delta\text{NO},\alpha\beta}(\mathbf{r},\mathbf{r}_{0})$

$\Gamma^{\alpha\alpha}(\mathbf{r},\mathbf{r}_{0})$

$\Gamma^{\Delta\text{NO},\alpha\alpha}(\mathbf{r},\mathbf{r}_{0})$

Figure 12:  FCI and $\Delta\text{NO}$ two-electron densities of H4 cluster
when a test electron is placed at a H nucleus, $\mathbf{r}_{0}$ (orange dot).
Densities are shown for H4 at $R=1.7$ Å and $\theta=$ $70^{\circ}$ (left),
$85^{\circ}$ (middle), and $90^{\circ}$ (right).

The two-electron densities are plotted for a test-electron placed at the upper
right H atom (indicated by an orange dot). Therefore, the plots show the
probability of finding the other electron(s), opposite or parallel-spin,
within the H4 cluster.

At $\theta=70^{\circ}$ (left column of Figure 12), the $\Delta\text{NO}$ and
FCI two-electron densities closely resemble each other. For opposite-spin
electrons (top half of Figure 12), if an $\alpha$ electron is found at the
upper-right H atom, then it is most likely to find a $\beta$ electron on the
closest H atom, to which it is paired. It would then be equally likely to find
the other $\beta$ electron on either of the two more distant H atoms. In the
case of parallel-spin electrons (bottom half of Figure 12), if an $\alpha$
electron is found at the upper-right H atom then the other $\alpha$ electron
must be found in the “other” H2 molecule, with an equal probability of being
found at either proton.

At $\theta=90^{\circ}$ (right column of Figure 12), the behaviour is much
different. For electrons of opposite-spin, when an $\alpha$ electron is found
at the upper-right H atom, it is equally likely to find the $\beta$ electron
at either of the nearest protons, with a reduced probability of finding a
$\beta$ electron at the furthest proton. Comparing the $\Delta\text{NO}$ and
FCI two-electron densities, the same qualitative behaviour is seen but the
$\Delta\text{NO}$ 2-RDM results in too high a probability of finding the
$\beta$ electron on the furthest proton compared to the probability at the
nearest protons. In the case of parallel-spin electrons, if an $\alpha$
electron is found at the upper-right H atom then it is most likely that
another $\alpha$ electron is on the furthest proton, with only a slight chance
of an $\alpha$ electron on the nearest protons. Again, $\Delta\text{NO}$
models this behaviour qualitatively, but compared to FCI, $\Delta\text{NO}$
predicts the probability of finding an $\alpha$ electron on the nearest
protons to be too high relative to the furthest proton. This suggests that the
amount of exchange, or Fermi correlation, which occurs between parallel-spin
electrons, in the $\Delta\text{NO}$ 2-RDM is slightly inadequate at
$\theta=90^{\circ}$.

Finally, at $\theta=85^{\circ}$ (middle column of Figure 12), the difference
between $\Delta\text{NO}$ and FCI two-electron densities is the most
pronounced. It is evident that the electrons of the $\Delta\text{NO}$ 2-RDM at
$\theta=85^{\circ}$ behave more like the electrons at $\theta=70^{\circ}$ (2H2
regime) than the actual behaviour given by the FCI two-electron density. In
the case of the $\Delta\text{NO}$ two-electron densities, both opposite and
parallel-spin, the probabilities of finding electrons on the left-side protons
(other H2) is near equivalent. Whereas, the FCI two-electron density shows
uneven probabilities for the left-side protons, which more resembles the
behaviour at $\theta=90^{\circ}$ (H4 regime).

The inability of $\Delta\text{NO}$ to accurately model the PEC from
$\theta=70^{\circ}$ to $90^{\circ}$, is due to the energetic preference of
$\Delta\text{NO}$ to model H4 as 2H2, at $\theta$ values close to
$90^{\circ}$. The error observed in the parallel-spin two-electron density
suggests that $\Delta\text{NO}$ is incorrectly describing the correlation
between different electron pairs (interpair correlation), which is similar to
the conclusion reached in a study with PNOF6.Ramos-Cordoba  _et al._ (2015)
For $\Delta\text{NO}$, this could mean an error in the exchange component
(Fermi correlation), or possibly something else, but it implies that the high-
spin correction (HSC) is not adequate. The HSC provides the correct
description in the strong correlation limit (both pairs), but it has
deficiencies for moderate amounts of static correlation and must be improved
upon.

The inadequacy of the HSC is highlighted further by examining the
$N$-representabiltity of the $\Delta\text{NO}$ 2-RDM; in particular, the
violation of the $PQG$-conditions (Subsection II.6). The sums of the negative
eigenvalues of the $\bm{P}$, $\bm{Q}$, and $\bm{G}$ matrices corresponding to
the $\Delta\text{NO}$ 2-RDM for rectangle-to-square H4 are shown in Figure 13.

Figure 13:  Sum of the negative eigenvalues of $\bm{P}$, $\bm{Q}$ and $\bm{G}$
for the rectangle-to-square H4 cluster as a function of $\theta$ at $R=1.7$ Å.

The sums are separated into contributions from different spin-blocks of the
$\bm{P}$, $\bm{Q}$, and $\bm{G}$ matrices. The $\bm{P}$ and $\bm{Q}$ matrices
are separated into spin-blocks for each specific spin-pair, whereas the
$\bm{G}$ matrix is separated into only three blocks,
$\bm{G}^{\alpha\beta\alpha\beta}$, $\bm{G}^{\beta\alpha\beta\alpha}$, and
$\bm{G}^{\sigma\sigma\sigma\sigma}$, due to the coupling of the parallel-spin
blocks [Equation (70c)]. For all values of $\theta$, the parallel-spin blocks
of the $\bm{P}$ and $\bm{Q}$ matrices have no negative eigenvalues. This is
completely expected considering that all parallel-spin terms of the
$\Delta\text{NO}$ 2-RDM are derived from the $\Delta\text{NO}$ wave function.
The only terms not originating from the $\Delta\text{NO}$ wave function are
those corresponding to the HSC, which are applied to the opposite-spin block
of the 2-RDM only [Equation (35)]. The off-diagonal,
$\bm{G}^{\alpha\alpha\beta\beta}$ and $\bm{G}^{\beta\beta\alpha\alpha}$, spin-
blocks of $\bm{G}^{\sigma\sigma\sigma\sigma}$ contain 2-RDM elements
corresponding to the HSC and they are responsible for the contribution to
$G^{\sigma\sigma\sigma\sigma}_{\text{neg}}$. The most severe violation of the
$PQG$-condition is seen for $\bm{G}^{\alpha\beta\alpha\beta}$, also known as
the particle-hole density matrix. The sum of the negative eigenvalues of
$\bm{G}^{\alpha\beta\alpha\beta}$ reach a minimum at $\theta=90^{\circ}$,
where $G^{\alpha\beta\alpha\beta}_{\text{neg}}=-0.229$. Due to the fact that
the HSC is applied to the exchange-like terms (i.e.
$\Gamma^{\alpha\beta}_{pqqp}$) of the $\Delta\text{NO}$ 2-RDM, the sum of the
negative eigenvalues of the $\bm{P}^{\alpha\beta\alpha\beta}$ and
$\bm{Q}^{\alpha\beta\alpha\beta}$ matrices are equivalent for all values of
$\theta$. In the case of $P_{\text{neg}}$ and $Q_{\text{neg}}$, the minimum
value observed is -0.0818, which occurs at $\theta=86^{\circ}$ and
$\theta=94^{\circ}$. As seen for the energy and orbital occupancies, the
violation of the $PQG$-conditions exhibits different behaviour, with respect
to $\theta$, in the neighbourhood of $\theta=90^{\circ}$ (between the grey
rectangles) compared to $\theta$ closer to $70^{\circ}$ or $110^{\circ}$
(outside the grey rectangles).

##  V Conclusions

The $\Delta\text{NO}$ method is one of many that combines different models for
different aspects of electron correlation in an effort to efficiently capture
the total. The $\Delta\text{NO}$-ODF approach is unique in the manner that it
combines an approximate 2-RDM for static correlation with an ODF for dynamic
correlation. With the exception of the high-spin correction, the
$\Delta\text{NO}$ 2-RDM and energy can be derived from a multideterminantal
wave function with restricted $2n$-tuple excitations. Combining the
$\Delta\text{NO}$ 2-RDM for static correlation with an ODF for dynamic
correlation requires the inclusion of a double-counting correction, which can
be done by replacing the short-range correlation counted by the
$\Delta\text{NO}$ 2-RDM.

Previously, $\Delta\text{NO}$, and other NOF methods, optimized the 2-RDM
through alternating optimization of the occupation numbers and the orbitals
until convergence was reached. The orbitals themselves were optimized by
iterative diagonalization of a pseudo-Fock matrix. The trust-region Newton
method introduced here, which allows for simultaneous optimization of the
orbitals and occupancies, drastically reduces the number of optimization steps
required and has the ability to identify and avoid saddle points.

Even simple hydrogen clusters can prove challenging for electronic structure
methods, particularly single-reference methods. Overall, the
$\Delta\text{NO}$-ODF approach performs well when describing the potential
energy surfaces of H2, H3 and H4, relative to the more accurate, but costly,
multireference wave function method, MRMP2.

Analysis of the two-electron density of H4 reveals that $\Delta\text{NO}$
incorrectly describes interorbital correlation. Despite being accurate in the
strong-correlation limit, the ad hoc HSC does have deficiencies, particularly
near the $D_{2h}\to D_{4h}$ transition. A more complete description of
interorbital correlation within $\Delta\text{NO}$ is a target of future
development.

###### Acknowledgements.

JWH thanks the Natural Sciences and Engineering Research Council of Canada
(NSERC) for a Discovery Grant, Compute/Calcul Canada for computing resources
and the Discovery Institute for Computation and Synthesis for useful
consultations.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon request.

##  Appendix A $\Delta\text{NO}$ energy derivatives

The optimization of the $\Delta\text{NO}$ 2-RDM is acheived by simultaneously
finding the optimal $\Delta\text{NO}$ orbitals, $\\{\phi_{p}\\}$, and electron
transfer variables, $\\{\Delta_{pq}\\}$. The optimization is performed using a
trust-region Newton algorithm, which requires the analytical calculation of
the energy gradient, $\mathbf{g}$, and hessian, $\mathbf{H}$, with respect to
the skew-matrix parameters, $\mathbf{y}$, and the variables used to
parameterize the electron-transfer variables, $\mathbf{\theta}$.

The gradient, with respect to $\mathbf{y}$, is evaluated at $\mathbf{y}=0$.
The first and second derivatives of the individual orbitals,

$\left.\frac{\partial\phi_{a}}{\partial
y_{pq}}\right|_{\mathbf{y}=0}=\delta_{aq}\phi_{p}-\delta_{ap}\phi_{q}$ (73)

and

$\displaystyle\left.\frac{\partial^{2}\phi_{a}}{\partial y_{pq}\partial
y_{rs}}\right|_{\mathbf{y}=0}$
$\displaystyle=\left(\delta_{qr}\delta_{as}-\delta_{qs}\delta_{ar}\right)\phi_{p}$
$\displaystyle-\left(\delta_{pr}\delta_{as}-\delta_{ps}\delta_{ar}\right)\phi_{q}$
$\displaystyle+\left(\delta_{ps}\delta_{aq}-\delta_{qs}\delta_{ap}\right)\phi_{r}$
$\displaystyle-\left(\delta_{pr}\delta_{aq}-\delta_{qr}\delta_{ap}\right)\phi_{s},$
(74)

can be used to derive expressions for derivatives of the total energy,
$E^{\Delta\text{NO}}$. The total gradient with respect to the orbital-mixing
parameters, is given by

$g^{\phi}_{pq}=2\left(\lambda_{qp}-\lambda_{pq}\right),$ (75)

where

$\displaystyle\lambda_{pq}$
$\displaystyle=(2-W_{p})n_{p}\left(h_{pq}+G_{pq}^{\text{cl}}\right)$
$\displaystyle+2(1-W_{p})n_{p}G_{pq}^{Opt.~{}Commun.~{}}+W_{p}n_{p}G_{pq}^{\text{op}}$
$\displaystyle+2\sum_{r}\eta_{pr}\left(2\langle pr|qr\rangle-\langle
pr|rq\rangle\right)$
$\displaystyle+2\sum_{r}\left(\zeta_{pr}-\xi_{pr}\right)\langle pq|rr\rangle$
$\displaystyle-2\sum_{r}\kappa_{pr}\langle pr|rq\rangle$ (76)

and the two-electron integrals are defined as follows

$\langle
pq|rs\rangle=\int\frac{\phi^{*}_{p}(\mathbf{r}_{1})\phi^{*}_{q}(\mathbf{r}_{2})\phi_{r}(\mathbf{r}_{1})\phi_{s}(\mathbf{r}_{2})}{r_{12}}d\mathbf{r}_{1}d\mathbf{r}_{2},$
(77)

and

$G_{pq}^{\text{cl}}=\sum_{\mu\nu}P_{\nu\mu}^{\text{cl}}\left(\langle
p\mu|q\nu\rangle-\frac{1}{2}\langle p\mu|\nu q\rangle\right),$ (78)
$G_{pq}^{\text{op}}=\sum_{\mu\nu}P_{\nu\mu}^{\text{op}}\left(\langle
p\mu|q\nu\rangle-\langle p\mu|\nu q\rangle\right),$ (79)
$G_{pq}^{Opt.~{}Commun.~{}}=\sum_{\mu\nu}P_{\nu\mu}^{\text{op}}\left(\langle
p\mu|q\nu\rangle-\frac{1}{2}\langle p\mu|\nu q\rangle\right),$ (80)

where $\mu$ and $\nu$ are indices of atomic basis functions. The density
matrices are defined as

$P_{\mu\nu}^{\text{cl}}=\sum^{\text{cl}}_{p}2n_{p}C_{\mu p}C_{\nu p}^{*}$ (81)

and

$P_{\mu\nu}^{\text{op}}=\sum^{\text{op}}_{p}n_{p}C_{\mu p}C_{\nu p}^{*},$ (82)

where $C_{\mu p}$ is a $\Delta\text{NO}$ orbital coefficient.

For the second derivatives with respect to the orbital-mixng parameters, it is
useful to decompose $E^{\Delta\text{NO}}$. The one-electron energy is given by

$E^{\Delta\text{NO}}_{\text{1e}}=\sum^{\text{cl}}_{p}2n_{p}h_{pp}+\sum^{\text{op}}_{p}n_{p}h_{pp},$
(83)

and the contribution of $E^{\Delta\text{NO}}_{\text{1e}}$ to the orbital
hessian is given by

$\displaystyle H^{\phi\phi,\text{1e}}_{pq,rs}=$
$\displaystyle\delta_{pr}\left(n^{\alpha}_{p}+n^{\beta}_{p}+n^{\alpha}_{r}+n^{\beta}_{r}-n^{\alpha}_{q}-n^{\beta}_{q}-n^{\alpha}_{s}-n^{\beta}_{s}\right)h_{qs}$
$\displaystyle+$
$\displaystyle\delta_{qs}\left(n^{\alpha}_{q}+n^{\beta}_{q}+n^{\alpha}_{s}+n^{\beta}_{s}-n^{\alpha}_{p}-n^{\beta}_{p}-n^{\alpha}_{r}-n^{\beta}_{r}\right)h_{pr}$
$\displaystyle+$
$\displaystyle\delta_{ps}\left(n^{\alpha}_{q}+n^{\beta}_{q}+n^{\alpha}_{r}+n^{\beta}_{r}-n^{\alpha}_{p}-n^{\beta}_{p}-n^{\alpha}_{s}-n^{\beta}_{s}\right)h_{qr}$
$\displaystyle+$
$\displaystyle\delta_{qr}\left(n^{\alpha}_{p}+n^{\beta}_{p}+n^{\alpha}_{s}+n^{\beta}_{s}-n^{\alpha}_{q}-n^{\beta}_{q}-n^{\alpha}_{r}-n^{\beta}_{r}\right)h_{ps}.$
(84)

The two-electron energy associated with the 0-1RDM component of the energy is
given by

$E^{\Delta\text{NO}}_{\text{0-1RDM-2e}}=E^{\Delta\text{NO}}_{\text{0-1RDM}}-E^{\Delta\text{NO}}_{\text{1e}},$
(85)

and the corresponding contribution to the orbital hessian is

$\displaystyle H^{\phi\phi,\text{0-1RDM-2e}}_{pq,rs}=$
$\displaystyle\delta_{pr}\Big{(}\left[n^{\alpha}_{p}+n^{\beta}_{p}+n^{\alpha}_{r}+n^{\beta}_{r}-n^{\alpha}_{q}-n^{\beta}_{q}-n^{\alpha}_{s}-n^{\beta}_{s}\right]G^{\text{cl}}_{qs}$
$\displaystyle+\left[n^{\alpha}_{p}+n^{\alpha}_{r}-n^{\alpha}_{q}-n^{\alpha}_{s}\right]G^{\text{op}}_{qs}$
$\displaystyle+\left[n^{\beta}_{p}+n^{\beta}_{r}-n^{\beta}_{q}-n^{\beta}_{s}\right]\left[2G^{Opt.~{}Commun.~{}}_{qs}-G^{\text{op}}_{qs}\right]\Big{)}$
$\displaystyle+\delta_{qs}\Big{(}\left[n^{\alpha}_{q}+n^{\beta}_{q}+n^{\alpha}_{s}+n^{\beta}_{s}-n^{\alpha}_{p}-n^{\beta}_{p}-n^{\alpha}_{r}-n^{\beta}_{r}\right]G^{\text{cl}}_{pr}$
$\displaystyle+\left[n^{\alpha}_{q}+n^{\alpha}_{s}-n^{\alpha}_{p}-n^{\alpha}_{r}\right]G^{\text{op}}_{pr}$
$\displaystyle+\left[n^{\beta}_{q}+n^{\beta}_{s}-n^{\beta}_{p}-n^{\beta}_{r}\right]\left[2G^{Opt.~{}Commun.~{}}_{pr}-G^{\text{op}}_{pr}\right]\Big{)}$
$\displaystyle+\delta_{ps}\Big{(}\left[n^{\alpha}_{q}+n^{\beta}_{q}+n^{\alpha}_{r}+n^{\beta}_{r}-n^{\alpha}_{p}-n^{\beta}_{p}-n^{\alpha}_{s}-n^{\beta}_{s}\right]G^{\text{cl}}_{qr}$
$\displaystyle+\left[n^{\alpha}_{q}+n^{\alpha}_{r}-n^{\alpha}_{p}-n^{\alpha}_{s}\right]G^{\text{op}}_{qr}$
$\displaystyle+\left[n^{\beta}_{q}+n^{\beta}_{r}-n^{\beta}_{p}-n^{\beta}_{s}\right]\left[2G^{Opt.~{}Commun.~{}}_{qr}-G^{\text{op}}_{qr}\right]\Big{)}$
$\displaystyle+\delta_{qr}\Big{(}\left[n^{\alpha}_{p}+n^{\beta}_{p}+n^{\alpha}_{s}+n^{\beta}_{s}-n^{\alpha}_{q}-n^{\beta}_{q}-n^{\alpha}_{r}-n^{\beta}_{r}\right]G^{\text{cl}}_{ps}$
$\displaystyle+\left[n^{\alpha}_{p}+n^{\alpha}_{s}-n^{\alpha}_{q}-n^{\alpha}_{r}\right]G^{\text{op}}_{qr}$
$\displaystyle+\left[n^{\beta}_{p}+n^{\beta}_{s}-n^{\beta}_{q}-n^{\beta}_{r}\right]\left[2G^{Opt.~{}Commun.~{}}_{qr}-G^{\text{op}}_{qr}\right]\Big{)}$
$\displaystyle+4\left[\left(n^{\alpha}_{p}-n^{\alpha}_{q}\right)\left(n^{\alpha}_{r}-n^{\alpha}_{s}\right)+\left(n^{\beta}_{p}-n^{\beta}_{q}\right)\left(n^{\beta}_{r}-n^{\beta}_{s}\right)\right]\langle
pr||qs\rangle$
$\displaystyle+4\left[\left(n^{\alpha}_{p}-n^{\alpha}_{q}\right)\left(n^{\beta}_{r}-n^{\beta}_{s}\right)+\left(n^{\beta}_{p}-n^{\beta}_{q}\right)\left(n^{\alpha}_{r}-n^{\alpha}_{s}\right)\right]\langle
pr|qs\rangle,$ (86)

where the antisymmetrized two-electron integrals are defined as

$\langle pr||qs\rangle=\langle pr|qs\rangle-\langle pr|sq\rangle.$ (87)

For the contributions to the orbital hessian from the remaining components of
$E^{\Delta\text{NO}}$, it is useful to consider a general expression for
component $X$,

$E^{\Delta\text{NO}}_{X}=\sum_{abcd}\Gamma_{abcd}^{X}\langle ab|cd\rangle$
(88)

where $\Gamma_{abcd}^{X}$ is a 2-RDM element, and in this case
$X=\text{pair},\text{stat}$ or HSC. The second derivative with respect to the
skew-matrix parameters for this general expression is given by

$\displaystyle H^{\phi\phi,X}_{pq,rs}=2\sum_{bcd}\Big{[}\big{(}$
$\displaystyle\delta_{ps}\Gamma^{X}_{rbcd}-\delta_{pr}\Gamma^{X}_{sbcd}\big{)}\langle
qb|cd\rangle$ $\displaystyle+\big{(}$
$\displaystyle\delta_{qr}\Gamma^{X}_{pbcd}-\delta_{pr}\Gamma^{X}_{qbcd}\big{)}\langle
sb|cd\rangle$ $\displaystyle+\big{(}$
$\displaystyle\delta_{ps}\Gamma^{X}_{qbcd}-\delta_{qs}\Gamma^{X}_{pbcd}\big{)}\langle
rb|cd\rangle$ $\displaystyle+\big{(}$
$\displaystyle\delta_{qr}\Gamma^{X}_{sbcd}-\delta_{qs}\Gamma^{X}_{rbcd}\big{)}\langle
pb|cd\rangle\Big{]}$ $\displaystyle+4\sum_{ab}\Big{[}$
$\displaystyle\Gamma^{X}_{prab}\langle qs|ab\rangle-\Gamma^{X}_{qrab}\langle
ps|ab\rangle$ $\displaystyle-$ $\displaystyle\Gamma^{X}_{psab}\langle
qr|ab\rangle+\Gamma^{X}_{qsab}\langle pr|ab\rangle$ $\displaystyle+$
$\displaystyle\Gamma^{X}_{parb}\langle qa|sb\rangle-\Gamma^{X}_{pasb}\langle
qa|rb\rangle$ $\displaystyle-$ $\displaystyle\Gamma^{X}_{qarb}\langle
pa|sb\rangle+\Gamma^{X}_{qasb}\langle pa|rb\rangle$ $\displaystyle+$
$\displaystyle\Gamma^{X}_{pabr}\langle qa|bs\rangle-\Gamma^{X}_{qabr}\langle
pa|bs\rangle$ $\displaystyle-$ $\displaystyle\Gamma^{X}_{pabs}\langle
qa|br\rangle+\Gamma^{X}_{qabs}\langle pa|br\rangle\Big{]}$ (89)

and can be found elsewhere.Bozkaya  _et al._ (2011) In the case of
$\Delta\text{NO}$, for which the 2-RDM is a JKL-functional, the above
expression can be simplified due to the fact that each term of the 2-RDM is
described by only two indices, rather than four. The simplification can be
applied specifically to terms involving $J$, $K$ or $L$ integrals. In the case
of Coulomb integrals, $J$, the 2-RDM simplifies as
$\Gamma_{abcd}=\delta_{ac}\delta_{bd}\gamma_{ab}$, giving

$E^{\Delta\text{NO}}_{X_{J}}=\sum_{ab}\gamma_{ab}J_{ab},$ (90)

and leading to the following orbital hessian contribution,

$\displaystyle H^{\phi\phi,{X_{J}}}_{pq,rs}=4\big{(}\gamma_{pr}$
$\displaystyle+\gamma_{qs}-\gamma_{ps}-\gamma_{qr}\big{)}\left(\langle
pr|qs\rangle+\langle ps|qr\rangle\right)$ $\displaystyle+2\sum_{a}\Big{[}$
$\displaystyle\delta_{ps}\left(\gamma_{qa}+\gamma_{ra}-\gamma_{pa}-\gamma_{sa}\right)\langle
qa|ra\rangle$ $\displaystyle+$
$\displaystyle\delta_{pr}\left(\gamma_{pa}+\gamma_{ra}-\gamma_{qa}-\gamma_{sa}\right)\langle
qa|sa\rangle$ $\displaystyle+$
$\displaystyle\delta_{qr}\left(\gamma_{pa}+\gamma_{sa}-\gamma_{qa}-\gamma_{ra}\right)\langle
pa|sa\rangle$ $\displaystyle+$
$\displaystyle\delta_{qs}\left(\gamma_{qa}+\gamma_{sa}-\gamma_{pa}-\gamma_{ra}\right)\langle
pa|ra\rangle\Big{]}$ (91)

In the case of exchange integrals, $K$, the 2-RDM simplifies as
$\Gamma_{abcd}=\delta_{ad}\delta_{bc}\gamma_{ab}$, giving

$E^{\Delta\text{NO}}_{X_{K}}=\sum_{ab}\gamma_{ab}K_{ab}$ (92)

and leading to the following orbital hessian contribution,

$\displaystyle H^{\phi\phi,{X_{K}}}_{pq,rs}=4\big{(}\gamma_{pr}$
$\displaystyle+\gamma_{qs}-\gamma_{ps}-\gamma_{qr}\big{)}\left(\langle
pr|sq\rangle+\langle ps|rq\rangle\right)$ $\displaystyle+2\sum_{a}\Big{[}$
$\displaystyle\delta_{ps}\left(\gamma_{qa}+\gamma_{ra}-\gamma_{pa}-\gamma_{sa}\right)\langle
qa|ar\rangle$ $\displaystyle+$
$\displaystyle\delta_{pr}\left(\gamma_{pa}+\gamma_{ra}-\gamma_{qa}-\gamma_{sa}\right)\langle
qa|as\rangle$ $\displaystyle+$
$\displaystyle\delta_{qr}\left(\gamma_{pa}+\gamma_{sa}-\gamma_{qa}-\gamma_{ra}\right)\langle
pa|as\rangle$ $\displaystyle+$
$\displaystyle\delta_{qs}\left(\gamma_{qa}+\gamma_{sa}-\gamma_{pa}-\gamma_{ra}\right)\langle
pa|ar\rangle\Big{]}$ (93)

Similarly, for time-inversion eschange integrals, $L$, the 2-RDM simplifies as
$\Gamma_{abcd}=\delta_{ab}\delta_{cd}\gamma_{ac}$ and the orbital hessian
contribution simplifies in the same manner.

Simultaneous optimization of the electron-transfer variables,
$\\{\Delta_{pq}\\}$, and consequently the orbital occupancies, requires
including the gradient and hessian contributions from $\mathbf{\theta}$, where
$\theta_{pq}$ is related to $\Delta_{pq}$ via Equation (58). The gradient
contributions are given by,

$g^{\theta}_{pq}=-\frac{\sin(2\theta_{pq})}{2}\frac{\partial
E^{\Delta\text{NO}}}{\partial\Delta_{pq}}$ (94)

The derivative of $E^{\Delta\text{NO}}_{\text{0-1RDM}}$ with respect to
$\Delta_{pq}$ is derived using the following identity,

$\frac{\partial n_{t}}{\partial\Delta_{pq}}=\delta_{qt}-\delta_{pt},$ (95)

which leads to

$\frac{\partial
E^{\Delta\text{NO}}_{\text{0-1RDM}}}{\partial\Delta_{pq}}=2\left(h_{qq}-h_{pp}+G_{qq}^{\text{cl}}-G_{pp}^{\text{cl}}+G_{qq}^{Opt.~{}Commun.~{}}-G_{pp}^{Opt.~{}Commun.~{}}\right).$
(96)

The contribution to the gradient from the remaining components of the energy
can be calculated by replacing the coefficients, $\eta_{tu}$, $\zeta_{tu}$,
$\xi_{tu}$ and $\kappa_{tu}$ [Equations (15), (20) and (25)] with their
corresponding derivatives. The derivatives of the coefficients are given by

$\displaystyle\frac{\partial\eta_{tu}}{\partial\Delta_{pq}}=\delta_{tu}\Bigg{[}$
$\displaystyle\delta_{qt}\left(1-2n_{q}\right)-\delta_{pt}\left(1-2n_{p}\right)\Bigg{]}$
$\displaystyle+\left(1-\delta_{tu}\right)\Bigg{[}$ $\displaystyle
O_{t}V_{u}\Big{\\{}\delta_{pt}\delta_{qu}\left(n_{u}-n_{t}-\Delta_{tu}\right)$
$\displaystyle\hphantom{O_{t}V_{t}}+\Delta_{tu}\left(\delta_{qu}+\delta_{pt}-\delta_{pt}\delta_{qu}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
V_{t}O_{u}\Big{\\{}\delta_{pu}\delta_{qt}\left(n_{t}-n_{u}-\Delta_{ut}\right)$
$\displaystyle\hphantom{O_{t}V_{t}}+\Delta_{ut}\left(\delta_{qt}+\delta_{pu}-\delta_{pu}\delta_{qt}\right)\Big{\\}}$
$\displaystyle+$ $\displaystyle
V_{t}V_{u}\left(\delta_{qt}\Delta_{pu}+\delta_{qu}\Delta_{pt}\right)\Bigg{]},$
(97)
$\frac{\partial\zeta_{tu}}{\partial\Delta_{pq}}=\frac{V_{t}V_{u}}{2}\left(\delta_{qt}\sqrt{\frac{\Delta_{pu}}{\Delta_{pt}}}+\delta_{qu}\sqrt{\frac{\Delta_{pt}}{\Delta_{pu}}}\right),$
(98) $\displaystyle\frac{\partial\xi_{tu}}{\partial\Delta_{pq}}$
$\displaystyle=\frac{O_{t}V_{u}}{2}\left(\delta_{pt}\delta_{qu}\sqrt{\frac{n_{t}}{\Delta_{tu}}}-\delta_{pt}\sqrt{\frac{\Delta_{tu}}{n_{t}}}\right)$
$\displaystyle+\frac{O_{u}V_{t}}{2}\left(\delta_{pu}\delta_{qt}\sqrt{\frac{n_{u}}{\Delta_{ut}}}-\delta_{pu}\sqrt{\frac{\Delta_{ut}}{n_{u}}}\right),$
(99)

and

$\displaystyle\frac{\partial\kappa_{tu}}{\partial\Delta_{pq}}$
$\displaystyle=\left(1-\delta_{tu}\right)\Bigg{[}\sum_{v\neq
w}\left(\xi_{uw}\frac{\partial\xi_{tv}}{\partial\Delta_{pq}}+\xi_{tv}\frac{\partial\xi_{uw}}{\partial\Delta_{pq}}\right)$
$\displaystyle+W_{t}\frac{n_{t}}{2\sqrt{2}}\sum_{v}\frac{\partial\xi_{uv}}{\partial\Delta_{pq}}+W_{u}\frac{n_{u}}{2\sqrt{2}}\sum_{v}\frac{\partial\xi_{tv}}{\partial\Delta_{pq}}\Bigg{]}.$
(100)

Using the chain rule, the $\mathbf{\theta}$ contribution to the hessian can be
expressed in terms of derivatives with respect to the $\\{\Delta_{pq}\\}$,

$\displaystyle H_{pq,rs}^{\theta\theta}$
$\displaystyle=\frac{\sin(2\theta_{pq})\sin(2\theta_{rs})}{4}\frac{\partial^{2}E^{\Delta\text{NO}}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle-\delta_{pr}\delta_{qs}\cos(2\theta_{pq})\frac{\partial
E^{\Delta\text{NO}}}{\partial\Delta_{pq}}$ (101)

In the case of the second derivative, it is useful to express the 0-1RDM and
pair contributions together, where

$E^{\Delta\text{NO}}_{\text{0-1RDM-
pair}}=E^{\Delta\text{NO}}_{\text{0-1RDM}}+E^{\Delta\text{NO}}_{\text{pair}},$
(102)

and

$\displaystyle\frac{\partial^{2}E^{\Delta\text{NO}}_{\text{0-1RDM-
pair}}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle=4\left(J_{pr}+J_{qs}-J_{ps}-J_{qr}\right)$
$\displaystyle-2\left(K_{pr}+K_{qs}-K_{ps}-K_{qr}\right)$ (103)

For the remainder of the contributions to the hessian, as with the gradient,
the second derivatives can be calculated by replacing the coefficients,
$\zeta_{tu}$, $\xi_{tu}$ and $\kappa_{tu}$, with their corresponding second
derivative,

$\displaystyle\frac{\partial^{2}\zeta_{tu}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle=\delta_{pr}\frac{V_{t}V_{u}}{4}\Bigg{(}\frac{\delta_{pt}\delta_{su}+\delta_{qu}\delta_{st}}{\sqrt{\Delta_{pt}\Delta_{pu}}}$
$\displaystyle-\delta_{qt}\delta_{st}\sqrt{\frac{\Delta_{pu}}{\Delta_{pt}^{3}}}-\delta_{qu}\delta_{su}\sqrt{\frac{\Delta_{pt}}{\Delta_{pu}^{3}}}\Bigg{)},$
(104)
$\displaystyle\frac{\partial^{2}\xi_{tu}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle=\delta_{pr}\frac{O_{t}V_{u}}{4}\Bigg{(}\frac{-\delta_{pt}\delta_{qu}-\delta_{pt}\delta_{su}}{\sqrt{n_{t}\Delta_{tu}}}$
$\displaystyle-\delta_{pt}\delta_{qu}\delta_{su}\sqrt{\frac{n_{t}}{\Delta_{tu}^{3}}}-\delta_{pt}\sqrt{\frac{\Delta_{tu}}{n_{t}^{3}}}\Bigg{)}$
$\displaystyle+\delta_{pr}\frac{O_{u}V_{t}}{4}\Bigg{(}\frac{-\delta_{pu}\delta_{qt}-\delta_{pu}\delta_{st}}{\sqrt{n_{u}\Delta_{ut}}}$
$\displaystyle-\delta_{pu}\delta_{qt}\delta_{st}\sqrt{\frac{n_{u}}{\Delta_{ut}^{3}}}-\delta_{pu}\sqrt{\frac{\Delta_{ut}}{n_{u}^{3}}}\Bigg{)}$
(105)

and

$\displaystyle\frac{\partial^{2}\kappa_{tu}}{\partial\Delta_{pq}\partial\Delta_{rs}}=\left(1-\delta_{tu}\right)\Bigg{[}\sum_{v\neq
w}\Big{(}$
$\displaystyle\xi_{uw}\frac{\partial^{2}\xi_{tv}}{\partial\Delta_{pq}\partial\Delta_{rs}}+\xi_{tv}\frac{\partial^{2}\xi_{uw}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle+$
$\displaystyle\frac{\partial\xi_{tv}}{\partial\Delta_{pq}}\frac{\partial\xi_{uw}}{\partial\Delta_{rs}}+\frac{\partial\xi_{uw}}{\partial\Delta_{pq}}\frac{\partial\xi_{tv}}{\partial\Delta_{rs}}\Big{)}$
$\displaystyle+W_{t}\frac{n_{t}}{2\sqrt{2}}\sum_{v}\frac{\partial^{2}\xi_{uv}}{\partial\Delta_{pq}\partial\Delta_{rs}}$
$\displaystyle+W_{u}\frac{n_{u}}{2\sqrt{2}}\sum_{v}\frac{\partial^{2}\xi_{tv}}{\partial\Delta_{pq}\partial\Delta_{rs}}\Bigg{]}.$
(106)

The hessian terms involving derivatives with respect to both $\mathbf{y}$ and
$\mathbf{\theta}$, can be expressed in terms of derivatives of
$\mathbf{\lambda}$,

$H^{\phi\theta}_{pq,rs}=-\sin(2\theta_{rs})\left(\frac{\partial\lambda_{qp}}{\partial\Delta_{rs}}-\frac{\partial\lambda_{pq}}{\partial\Delta_{rs}}\right).$
(107)

The derivative of the $\lambda_{pq}$ component corresponding to
$E^{\Delta\text{NO}}_{\text{0-1RDM}}$ [first three terms of Equation (A)] is
given by

$\displaystyle\frac{\partial\lambda^{\text{0-1RDM}}_{pq}}{\partial\Delta_{rs}}=\big{(}$
$\displaystyle 2-W_{p}\big{)}$ $\displaystyle\times\Big{[}$
$\displaystyle\left(\delta_{ps}-\delta_{pr}\right)\left(h_{pq}+G_{pq}^{\text{cl}}+(1-W_{p})G_{pq}^{Opt.~{}Commun.~{}}+W_{p}G_{pq}^{\text{op}}\right)$
$\displaystyle+2\left(\langle ps|qs\rangle-\langle
pr|qr\rangle\right)-\left(\langle ps|sq\rangle-\langle
pr|rq\rangle\right)\Big{]}.$ (108)

The remaining components of $\frac{\partial\lambda_{pq}}{\partial\Delta_{rs}}$
are derived by substituting the coefficients, $\zeta_{tu}$, $\xi_{tu}$ and
$\kappa_{tu}$, with their corresponding derivative.

## References

  * Cizek (1969) J. Cizek, Adv. Chem. Phys. 14, 35 (1969).
  * Jeziorski and Monkhorst (1981) B. Jeziorski and H. J. Monkhorst, Phys. Rev. A 24, 1668 (1981).
  * Andersson  _et al._ (1990) K. Andersson, P. Malmqvist, B. O. Roos, A. J. Sadlej, and K. Wolinski, J. Phys. Chem. 94, 5483 (1990).
  * Mok, Neumann, and Handy (1996) D. K. W. Mok, R. Neumann, and N. C. Handy, J. Phys. Chem. 100, 6225 (1996).
  * Handy and Cohen (2001) N. C. Handy and A. J. Cohen, Mol. Phys. 99, 403 (2001).
  * Cremer (2001) D. Cremer, Mol. Phys. 99, 1899 (2001).
  * Pulay (2011) P. Pulay, Int. J. Quantum Chem. 111, 3273 (2011).
  * Becke (2013) A. D. Becke, J. Chem. Phys. 138, 074109 (2013).
  * Crittenden (2013) D. L. Crittenden, J. Phys. Chem. A 117, 3852 (2013).
  * Tsuchimochi and Voorhis (2014) T. Tsuchimochi and T. V. Voorhis, J. Chem. Phys. 141, 164117 (2014).
  * Wallace and Crittenden (2014) A. J. Wallace and D. L. Crittenden, J. Phys. Chem. A 118, 2138 (2014).
  * Hollett, Hosseini, and Menzies (2016) J. W. Hollett, H. Hosseini, and C. Menzies, J. Chem. Phys. 145, 084106 (2016).
  * Ramos-Cordoba, Salvador, and Matito (2016) E. Ramos-Cordoba, P. Salvador, and E. Matito, Phys. Chem. Chem. Phys. 18, 24015 (2016).
  * Benavides-Riveros, Lathiotakis, and Marques (2017) C. L. Benavides-Riveros, N. N. Lathiotakis, and M. A. L. Marques, Phys. Chem. Chem. Phys. 19, 12655 (2017).
  * Evangelista (2018) F. A. Evangelista, J. Chem. Phys. 149, 030901 (2018).
  * Via-Nadal  _et al._ (2019) M. Via-Nadal, M. Rodriguez-Mayorga, E. Ramos-Cordoba, and E. Matito, J. Phys. Chem. Lett. 10, 4032 (2019).
  * Fromager, Toulouse, and Jørgen Aa. Jensen (2007) E. Fromager, J. Toulouse, and H. Jørgen Aa. Jensen, J. Chem. Phys. 126, 074111 (2007).
  * Weimer, Sala, and Görling (2008) M. Weimer, F. D. Sala, and A. Görling, J. Chem. Phys. 128, 144109 (2008).
  * Toulouse  _et al._ (2009) J. Toulouse, I. C. Gerber, G. Jansen, A. Savin, and J. G. Angyán, Phys. Rev. Lett. 102, 096404 (2009).
  * Janesko, Henderson, and Scuseria (2009) B. G. Janesko, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 131, 034110 (2009).
  * Chai and Head-Gordon (2009) J.-D. Chai and M. Head-Gordon, J. Chem. Phys. 131, 174105 (2009).
  * Stoyanova  _et al._ (2013) A. Stoyanova, A. M. Teale, J. Toulouse, T. Helgaker, and E. Fromager, J. Chem. Phys. 139, 134113 (2013).
  * Yamanaka  _et al._ (2006) S. Yamanaka, K. Nakata, T. Ukai, T. Takada, and K. Yamaguchi, Int. J. Quantum Chem. 106, 3312 (2006).
  * Rohr, Toulouse, and Pernal (2010) D. R. Rohr, J. Toulouse, and K. Pernal, Phys. Rev. A 82, 052502 (2010).
  * Piris (2017) M. Piris, Phys. Rev. Lett. 119, 063002 (2017).
  * Piris (2018) M. Piris, Phys. Rev. A 98, 022504 (2018).
  * Pérez-Jiménez, Pérez-Jordá, and Sancho-Garcia (2007) Á. J. Pérez-Jiménez, J. M. Pérez-Jordá, and J. C. Sancho-Garcia, J. Chem. Phys. 127, 104102 (2007).
  * Manni  _et al._ (2014) G. L. Manni, R. K. Carlson, S. Luo, D. Ma, J. Olsen, D. G. Truhlar, and L. Gagliardi, J. Chem. Theory Comput. 10, 3669 (2014).
  * Carlson, Truhlar, and Gagliardi (2015) R. K. Carlson, D. G. Truhlar, and L. Gagliardi, J. Chem. Theory Comput. 11, 4077 (2015).
  * Gagliardi  _et al._ (2017) L. Gagliardi, D. G. Truhlar, G. L. Manni, R. K. Carlson, C. E. Hoyer, and J. L. Bao, Acc. Chem. Res. 50, 66 (2017).
  * Malcolm and McDouall (1998) N. O. J. Malcolm and J. J. W. McDouall, Chem. Phys. Lett. 282, 121 (1998).
  * Köhn  _et al._ (2013) A. Köhn, M. Hanauer, L. A. Mück, T.-C. Jagau, and J. Gauss, WIREs Comput. Mol. Sci. 3, 176 (2013).
  * Hubert, Hedegård, and H. J (2016) M. Hubert, E. D. Hedegård, and A. J. H. J, J. Chem. Theory Comput. 12, 2203 (2016).
  * Hollett and Loos (2019) J. W. Hollett and P.-F. Loos, J. Chem. Phys. in press (2019).
  * Hirao (1992a) K. Hirao, Chem. Phys. Lett. 190, 374 (1992a).
  * Grimme and Waletzke (1999) S. Grimme and M. Waletzke, J. Chem. Phys. 111, 5645 (1999).
  * Gräfenstein and Cremer (2000) J. Gräfenstein and D. Cremer, Chem. Phys. Lett. 316, 569 (2000).
  * Angeli, Cimiraglia, and Malrieu (2001) C. Angeli, R. Cimiraglia, and J. Malrieu, Chem. Phys. Lett. 350, 297 (2001).
  * Angeli, Cimiraglia, and Malrieu (2002) C. Angeli, R. Cimiraglia, and J. Malrieu, J. Chem. Phys. (2002).
  * Mostafanejad, Liebenthal, and DePrince III (2020) M. Mostafanejad, M. D. Liebenthal, and A. E. DePrince III, Journal of chemical theory and computation 16, 2274 (2020).
  * Furche (2008) F. Furche, The Journal of Chemical Physics 129, 114105 (2008).
  * Paldus  _et al._ (1993) J. Paldus, P. Piecuch, L. Pylypow, and B. Jeziorski, Physical Review A 47, 2738 (1993).
  * Kowalski and Jankowski (1998a) K. Kowalski and K. Jankowski, Physical Review Letters 81, 1195 (1998a).
  * Kowalski and Jankowski (1998b) K. Kowalski and K. Jankowski, Chemical physics letters 290, 180 (1998b).
  * Schipper, Gritsenko, and Baerends (1999) P. Schipper, O. Gritsenko, and E. Baerends, The Journal of chemical physics 111, 4056 (1999).
  * Jankowski and Kowalski (1999) K. Jankowski and K. Kowalski, The Journal of chemical physics 110, 3714 (1999).
  * Van Voorhis and Head-Gordon (2000) T. Van Voorhis and M. Head-Gordon, The Journal of Chemical Physics 113, 8873 (2000).
  * Hachmann, Cardoen, and Chan (2006) J. Hachmann, W. Cardoen, and G. K.-L. Chan, The Journal of chemical physics 125, 144101 (2006).
  * Al-Saidi, Zhang, and Krakauer (2007) W. A. Al-Saidi, S. Zhang, and H. Krakauer, The Journal of chemical physics 127, 144101 (2007).
  * Rassolov and Xu (2007) V. A. Rassolov and F. Xu, The Journal of chemical physics 127, 044104 (2007).
  * Sinitskiy, Greenman, and Mazziotti (2010) A. V. Sinitskiy, L. Greenman, and D. A. Mazziotti, The Journal of chemical physics 133, 014104 (2010).
  * Stella  _et al._ (2011) L. Stella, C. Attaccalite, S. Sorella, and A. Rubio, Physical Review B 84, 245117 (2011).
  * Bytautas  _et al._ (2011) L. Bytautas, T. M. Henderson, C. A. Jiménez-Hoyos, J. K. Ellis, and G. E. Scuseria, J. Chem. Phys. 135, 044119 (2011).
  * Ramos-Cordoba  _et al._ (2015) E. Ramos-Cordoba, X. Lopez, M. Piris, and E. Matito, The Journal of chemical physics 143, 164112 (2015).
  * Bulik, Henderson, and Scuseria (2015) I. W. Bulik, T. M. Henderson, and G. E. Scuseria, Journal of chemical theory and computation 11, 3171 (2015).
  * Robinson and Knowles (2012a) J. B. Robinson and P. J. Knowles, The Journal of chemical physics 136, 054114 (2012a).
  * Robinson and Knowles (2012b) J. B. Robinson and P. J. Knowles, The Journal of chemical physics 137, 054301 (2012b).
  * Robinson and Knowles (2012c) J. B. Robinson and P. J. Knowles, Journal of chemical theory and computation 8, 2653 (2012c).
  * Robinson and Knowles (2012d) J. B. Robinson and P. J. Knowles, Physical Chemistry Chemical Physics 14, 6729 (2012d).
  * Sand and Mazziotti (2013) A. M. Sand and D. A. Mazziotti, Computational and Theoretical Chemistry 1003, 44 (2013).
  * Limacher  _et al._ (2013) P. A. Limacher, P. W. Ayers, P. A. Johnson, S. De Baerdemacker, D. Van Neck, and P. Bultinck, Journal of chemical theory and computation 9, 1394 (2013).
  * Kats and Manby (2013) D. Kats and F. R. Manby, The Journal of chemical physics 139, 021102 (2013).
  * Henderson  _et al._ (2014) T. M. Henderson, I. W. Bulik, T. Stein, and G. E. Scuseria, The Journal of chemical physics 141, 244104 (2014).
  * Motta  _et al._ (2017) M. Motta, D. M. Ceperley, G. K.-L. Chan, J. A. Gomez, E. Gull, S. Guo, C. A. Jiménez-Hoyos, T. N. Lan, J. Li, F. Ma, _et al._ , Physical Review X 7, 031059 (2017).
  * Motta  _et al._ (2020) M. Motta, C. Genovese, F. Ma, Z.-H. Cui, R. Sawaya, G. K.-L. Chan, N. Chepiga, P. Helms, C. Jiménez-Hoyos, A. J. Millis, _et al._ , Physical Review X 10, 031058 (2020).
  * Giner  _et al._ (2020) E. Giner, A. Scemama, P.-F. Loos, and J. Toulouse, The Journal of chemical physics 152, 174104 (2020).
  * Marie, Kossoski, and Loos (2021) A. Marie, F. Kossoski, and P.-F. Loos, The Journal of Chemical Physics 155, 104105 (2021).
  * Sand, Truhlar, and Gagliardi (2017) A. M. Sand, D. G. Truhlar, and L. Gagliardi, J. Chem. Phys. 146, 034101 (2017).
  * Perdew, Burke, and Ernzerhof (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Garza  _et al._ (2015) A. J. Garza, I. W. Bulik, T. M. Henderson, and G. E. Scuseria, J. Chem. Phys. 142, 044109 (2015).
  * Stein, Henderson, and Scuseria (2014) T. Stein, T. M. Henderson, and G. E. Scuseria, The Journal of chemical physics 140, 214113 (2014).
  * Poelmans  _et al._ (2015) W. Poelmans, M. Van Raemdonck, B. Verstichel, S. De Baerdemacker, A. Torre, L. Lain, G. E. Massaccesi, D. R. Alcoba, P. Bultinck, and D. Van Neck, Journal of chemical theory and computation 11, 4064 (2015).
  * Chen, Zhou, and Wu (2015) Z. Chen, C. Zhou, and W. Wu, Journal of chemical theory and computation 11, 4102 (2015).
  * Limacher (2015) P. A. Limacher, Journal of chemical theory and computation 11, 3629 (2015).
  * Van Meer, Gritsenko, and Baerends (2018) R. Van Meer, O. Gritsenko, and E. Baerends, The Journal of chemical physics 148, 104102 (2018).
  * Brzek  _et al._ (2019) F. Brzek, K. Boguslawski, P. Tecmer, and P. S. Zuchowski, Journal of chemical theory and computation 15, 4021 (2019).
  * Vu and DePrince III (2020) N. Vu and A. E. DePrince III, The Journal of Chemical Physics 152, 244103 (2020).
  * Kutzelnigg and Mukherjee (1999) W. Kutzelnigg and D. Mukherjee, J. Chem. Phys. 110, 2800 (1999).
  * Szabo and Ostlund (1996) A. Szabo and N. S. Ostlund,  _Modern quantum chemistry_ (Dover, New York, 1996).
  * Bytautas, Scuseria, and Ruedenberg (2015) L. Bytautas, G. E. Scuseria, and K. Ruedenberg, The Journal of chemical physics 143, 094105 (2015).
  * Hollett and Pegoretti (2018) J. W. Hollett and N. Pegoretti, J. Chem. Phys. 148, 164111 (2018).
  * Colle and Salvetti (1975) R. Colle and O. Salvetti, Theor. Chim. Acta 37, 329 (1975).
  * Colle and Salvetti (1979) R. Colle and O. Salvetti, Theoret. Chim. Acta 53, 55 (1979).
  * Note (1) This value differs from the value determined in reference 81, $\kappa=1.73$, because it was found to overestimate the correlation energy between electrons of parallel-spin in these H clusters.
  * Helgaker, Jorgensen, and Olsen (2014) T. Helgaker, P. Jorgensen, and J. Olsen,  _Molecular electronic-structure theory_ (John Wiley & Sons, 2014).
  * Bozkaya  _et al._ (2011) U. Bozkaya, J. M. Turney, Y. Yamaguchi, H. F. Schaefer III, and C. D. Sherrill, The Journal of chemical physics 135, 104103 (2011).
  * Bozkaya and Sherrill (2013) U. Bozkaya and C. D. Sherrill, The Journal of chemical physics 139, 054104 (2013).
  * Bozkaya and Sherrill (2014) U. Bozkaya and C. D. Sherrill, The Journal of chemical physics 141, 204105 (2014).
  * Bozkaya (2016) U. Bozkaya, Journal of chemical theory and computation 12, 1179 (2016).
  * Nocedal and Wright (2006) J. Nocedal and S. J. Wright,  _Numerical Optimization_ (Springer, New York, New York, 2006).
  * Garrod and Percus (1964) C. Garrod and J. K. Percus, Journal of Mathematical Physics 5, 1756 (1964).
  * Mazziotti (2006) D. A. Mazziotti, Accounts of chemical research 39, 207 (2006).
  * Mazziotti (2012) D. A. Mazziotti, Phys. Rev. Lett. 108, 263002 (2012).
  * Li, Liebenthal, and DePrince (2021) R. R. Li, M. Liebenthal, and E. DePrince, The Journal of Chemical Physics (2021).
  * Mazziotti (2002) D. A. Mazziotti, Physical Review A 65, 062511 (2002).
  * Rodríguez-Mayorga  _et al._ (2017) M. Rodríguez-Mayorga, E. Ramos-Cordoba, M. Via-Nadal, M. Piris, and E. Matito, Physical Chemistry Chemical Physics 19, 24029 (2017).
  * Piecuch  _et al._ (2002) P. Piecuch, S. A. Kucharski, K. Kowalski, and M. Musiał, Computer Physics Communications 149, 71 (2002).
  * Piecuch and Włoch (2005) P. Piecuch and M. Włoch, The Journal of chemical physics 123, 224105 (2005).
  * Becke (1993) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
  * Stephens  _et al._ (1994) P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, The Journal of physical chemistry 98, 11623 (1994).
  * Hertwig and Koch (1997) R. H. Hertwig and W. Koch, Chemical Physics Letters 268, 345 (1997).
  * Chaban, Schmidt, and Gordon (1997) G. Chaban, M. W. Schmidt, and M. S. Gordon, Theoretical Chemistry Accounts 97, 88 (1997).
  * Hirao (1992b) K. Hirao, Chemical physics letters 190, 374 (1992b).
  * Hirao (1992c) K. Hirao, Chemical physics letters 196, 397 (1992c).
  * Hirao (1992d) K. Hirao, International Journal of Quantum Chemistry 44, 517 (1992d).
  * Hirao (1993) K. Hirao, Chemical physics letters 201, 59 (1993).
  * Barca  _et al._ (2020) G. M. J. Barca, C. Bertoni, L. Carrington, D. Datta, N. De Silva, J. E. Deustua, D. G. Fedorov, J. R. Gour, A. O. Gunina, E. Guidez, T. Harville, S. Irle, J. Ivanic, K. Kowalski, S. S. Leang, H. Li, W. Li, J. J. Lutz, I. Magoulas, J. Mato, V. Mironov, H. Nakata, B. Q. Pham, P. Piecuch, D. Poole, S. R. Pruitt, A. P. Rendell, L. B. Roskop, K. Ruedenberg, T. Sattasathuchana, M. W. Schmidt, J. Shen, L. Slipchenko, M. Sosonkina, V. Sundriyal, A. Tiwari, J. L. Galvez Vallejo, B. Westheimer, M. Wloch, P. Xu, F. Zahariev, and M. S. Gordon, The Journal of Chemical Physics 152, 154102 (2020).
  * Garniron  _et al._ (2019) Y. Garniron, K. Gasperich, T. Applencourt, A. Benali, A. Ferte, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J. P. Malrieu, M. Veril, M. Caffarel, P. F. Loos, E. Giner, and A. Scemama, J. Chem. Theory Comput. 15, 3591 (2019).
  * Pritchard  _et al._ (2019) B. P. Pritchard, D. Altarawy, B. Didier, T. D. Gibson, and T. L. Windus, J. Chem. Inf. Model. 59, 4814 (2019).
  * Dunning Jr. (1989) T. H. Dunning Jr., J. Chem. Phys. 90, 1007 (1989).
  * Poirier, Hollett, and Warburton (2015) R. A. Poirier, J. W. Hollett, and P. L. Warburton, MUNgauss. Memorial University, Chemistry Department, St. John’s, NL A1B 3X7 (2015), with contributions from A. Alrawashdeh, J.-P. Becker, J. Besaw, S.D. Bungay, F. Colonna, A. El-Sherbiny, T. Gosse, D. Keefe, A. Kelly, D. Nippard, C.C. Pye, D. Reid, K. Saputantri, M. Shaw, M. Staveley, O. Stueker, Y. Wang, and J. Xidos.
  * Moré and Sorensen (1983) J. J. Moré and D. C. Sorensen, SIAM Journal on scientific and statistical computing 4, 553 (1983).
  * Gill, Johnson, and Pople (1993) P. M. W. Gill, B. G. Johnson, and J. A. Pople, Chem. Phys. Lett. 209, 506 (1993).
  * Becke (1988) A. D. Becke, J. Chem. Phys. 88, 2547 (1988).
  * Lew-Yee, Piris, and M. del Campo (2021) J. F. H. Lew-Yee, M. Piris, and J. M. del Campo, The Journal of Chemical Physics 154, 064102 (2021).
  * Piris and Ugalde (2009) M. Piris and J. M. Ugalde, J. Comput. Chem. 30, 2078 (2009).