# Tensor Hypercontraction Form of the Perturbative Triples Energy in Coupled-
Cluster Theory

Andy Jiang  Center for Computational Quantum Chemistry, Department of
Chemistry, University of Georgia, Athens, GA 30602  Center for Computational
Molecular Science and Technology, School of Chemistry and Biochemistry, School
of Computational Science and Engineering, Georgia Institute of Technology,
Atlanta, GA 30332-0400  Justin M. Turney<EMAIL_ADDRESS>Center for
Computational Quantum Chemistry, Department of Chemistry, University of
Georgia, Athens, GA 30602  Henry F. Schaefer III<EMAIL_ADDRESS>Center for
Computational Quantum Chemistry, Department of Chemistry, University of
Georgia, Athens, GA 30602

###### Abstract

We present the working equations for a reduced-scaling method of evaluating
the perturbative triples (T) energy in coupled-cluster theory, through the
tensor hypercontraction (THC) of the triples amplitudes ($t_{ijk}^{abc}$).
Through our method we can reduce the scaling of the (T) energy from the
traditional $\mathcal{O}(N^{7})$ to a more modest $\mathcal{O}(N^{5})$. We
also discuss implementation details to aid future research, development, and
software realization of this method. Additionally, we show that this method
yields sub-millihartree (mEh) differences from CCSD(T) when evaluating
absolute energies, and sub-0.1 kcal/mol energy differences when evaluating
relative energies. Finally, we demonstrate that this method converges to the
true CCSD(T) energy through the systematic increasing of the rank or
eigenvalue tolerance of the orthogonal projector, as well as exhibiting sub-
linear to linear error growth with respect to system size.

##  I Introduction

Coupled-cluster (CC) theory Crawford and H. F. Schaefer (2007); Bartlett and
Musial (2007) is one of the most important advances of modern quantum
chemistry, allowing for a polynomial-time evaluation of the electronic
energies and wavefunction of a molecule, as a size-extensive alternative to
truncated configuration interaction (CI) methods Cramer (2002); Sherrill and
H. F. Schaefer (1999). Truncated CC methods also avoid the intractable super-
exponential scaling of full configuration interaction (FCI), yielding
reasonable and chemically accurate relative energies compared to both the FCI
limit and to experimental results, especially in the context of CCSD(T), also
known as the “gold standard” method in computational quantum chemistry
Raghavachari  _et al._ (1989). The tractability and accuracy of CC methods
make the development of efficient CC methods crucial for the future of quantum
chemistry, as evaluation of accurate energies and wavefunctions is made
possible for larger and more complex systems through hardware advances such as
massively parallel computing Hirata (2003); Auer  _et al._ (2006); Janowski
_et al._ (2007); Janowski and Pulay (2008); van Dam  _et al._ (2011); Deumens
_et al._ (2011); Kobayashi and Rendell (1997); Anisimov  _et al._ (2014);
Solomonik  _et al._ (2014); Peng  _et al._ (2016); Lyakh (2019); Gyevi-Nagy
_et al._ (2020); Peng  _et al._ (2020); Datta and Gordon (2021); Gyevi-Nagy
_et al._ (2021); Kowalski  _et al._ (2021); Calvin  _et al._ (2021) and GPUs
Seritan  _et al._ (2020); Wang  _et al._ (2020); Peng  _et al._ (2019);
Kaliman and Krylov (2017); DePrince  _et al._ (2014); Ma  _et al._ (2011);
DePrince and Hammond (2011).

However, there is still a tremendous gap in applicability between coupled-
cluster theories (formally scaling at least $\mathcal{O}(N^{6})$) and lower-
scaling methods like Møller–Plesset perturbation theory (MP2) Møller and
Plesset (1934); Cremer (2011) and density functional theory (DFT) Hohenberg
and Kohn (1964); Kohn and Sham (1965) (scaling $\mathcal{O}(N^{5})$ or
better). Because of this, DFT and MP2 can be run on system tens or even
hundreds of times the size of a system typically evaluated with CC methods
Maurer  _et al._ (2014); Dawson  _et al._ (2022). To close the gap between CC
and less reliable electron correlation methods, it is useful to devise
approximation schemes to CC which reduce the scaling, but also allow a means
to systematically control the error compared to the non-approximated CC
method. One such approach involves local-correlation Li  _et al._ (2002, 2006,
2009); Neese  _et al._ (2009); Li and Piecuch (2010a, b); Rolik and Kállay
(2011); Rolik  _et al._ (2013); Riplinger and Neese (2013); Riplinger  _et
al._ (2013); Liakos  _et al._ (2015); Schwilk  _et al._ (2017), such as used
in the DLPNO methods Pinski  _et al._ (2015); Riplinger  _et al._ (2016). With
large enough molecules, these methods achieve asymptotic linear-scaling.

Another approach is the rank reduction of the coupled-cluster amplitudes
Parrish  _et al._ (2019), using orthogonal projectors that transform the
single and double cluster-amplitudes into a smaller basis

$\displaystyle T^{V}$ $\displaystyle=U^{V}_{ia}t_{i}^{a}$ (1) $\displaystyle
T^{VW}$ $\displaystyle=U^{V}_{ia}t_{ij}^{ab}U^{W}_{jb}\ .$ (2)

Because of the orthogonal nature of the projectors, getting the full
amplitudes from the rank-reduced form is trivial

$\displaystyle t_{i}^{a}$ $\displaystyle=U^{V}_{ia}T^{V}$ (3) $\displaystyle
t_{ij}^{ab}$ $\displaystyle=U^{V}_{ia}T^{VW}U^{W}_{jb}\ .$ (4)

As shown by Parrish and co-workers, the size of the $V$ and $W$ indices, also
known as the projector rank, can be made directly proportional to the system
size, while maintaining a set relative error from the absolute energy of a
molecule Parrish  _et al._ (2019). More recently, Hohenstein et. al. have
shown how to create a tensor hypercontracted (THC) form of the $t_{ij}^{ab}$
amplitudes, through the CANCENCOMP/PARAFAC (CP) decomposition Kolda and Bader
(2009) of the orthogonal projectors Hohenstein  _et al._ (2022).

$\displaystyle U^{V}_{ia}$ $\displaystyle=\sum_{X}y_{i}^{X}y_{a}^{X}\tau_{VX}$
(5) $\displaystyle t_{ij}^{ab}$
$\displaystyle=\sum_{XY}y_{i}^{X}y_{a}^{X}\widetilde{T}^{XY}y_{j}^{Y}y_{b}^{Y}$
(6) $\displaystyle\widetilde{T}^{XY}$
$\displaystyle=\sum_{VW}\tau_{VX}T^{VW}\tau_{WY}$ (7)

Parrish and Hohenstein have also shown that, in the context of CCSD, the size
of the $X$ index can be made proportional to the system size to maintain a set
relative error. Rank-reduction methods have also been applied to coupled-
cluster theories involving higher levels of excitation, recently by Lesiuk
with the SVD-CCSDT method Lesiuk (2020), where the concept of orthogonal
projectors is used to approximate the triples amplitude in CCSDT theory

$t_{ijk}^{abc}=U^{U}_{ia}U^{V}_{jb}U^{W}_{kc}T^{UVW}\ .$ (8)

In the following sections, we will combine the concepts of orthogonal
projectors and THC to develop working equations for a reduced-scaling variant
of the non-iterative perturbative triples correction to the CCSD energy
Raghavachari  _et al._ (1989). Recently, Lesiuk derived an $O(N^{6})$ approach
to the (T) energy with orthogonal projectors which he calls RR-CCSD(T) Lesiuk
(2022). In the current paper, we will improve upon the work of Lesiuk’s
approach utilizing tensor hypercontraction. Similar to how the THC-CCSD method
Hohenstein  _et al._ (2022) improves upon the RR-CCSD method Parrish  _et al._
(2019); Lesiuk (2022), our new approach, which we name THC-CCSD(T), will
commensurately enhance RR-CCSD(T), reducing the scaling of Lesiuk’s from
$\mathcal{O}(N^{6})$ to $\mathcal{O}(N^{5})$. For consistency, we will use
many of the same formalisms as Lesiuk Lesiuk (2020) and Hohenstein Hohenstein
_et al._ (2022).

##  II Theory

###  II.1 Notation

We will use the following conventions to describe the indices appearing in
this work:

  * •

$i,j,k,l$: Occupied molecular orbitals, which ranges from 1 to $n_{occ}$.

  * •

$a,b,c,d$: Virtual molecular orbitals, which ranges from 1 to $n_{virt}$.

  * •

$P,Q$: Auxiliary indices of density-fitted/Cholesky-decomposed ERIs, which
ranges from 1 to $n_{aux}$.

  * •

$w,v$: Laplace denominator weight indices, which ranges from 1 to $n_{w}$.

  * •

$U,V,W$: Rank-reduced dimension of the doubles orthogonal projector, which
ranges from 1 to $n_{proj}$.

  * •

$A,B,C$: Rank-reduced dimension of the triples orthogonal projector, which
ranges from 1 to $n_{proj}$.

  * •

$X,Y,Z$: CP-decomposition rank of the triples orthogonal projector, which
ranges from 1 to $n_{proj}$.

The relative sizes of the indices are as follows:

$n_{occ}<n_{virt}<n_{aux}\approx n_{proj}$ (9)

Note that $n_{w}$ does not grow with increasing molecular system size, and
therefore, run-time analysis of intermediates with $w,v$ indices will only
treat the Laplace index as a prefactor.

The frozen-core approximation was used in all post-Hartree–Fock computations
in this work; i.e., the 1$s$ electrons are not correlated for all first-row
atoms. The occupied space $n_{occ}$ always refers to the number of correlated
occupied orbitals. Einstein summation convention is used throughout – all
indices appearing on the right-hand side but not on the left-hand side of an
expression are summed over.

###  II.2 Perturbative Triples Correction to CCSD

CCSD is often not sufficient to obtain “chemically-reliable” theoretical
predictions, and it has been shown that only after triple excitations are
considered that relative energies of under 1 kcal/mol can be regularly
achieved Riley  _et al._ (2010); Karton  _et al._ (2006); Tajti  _et al._
(2004); Bak  _et al._ (2000); Hopkins and Tschumper (2004); Bartlett  _et al._
(1990). However, an explicit treatment of all triples has a very high cost of
$\mathcal{O}(N^{8})$. Therefore, the triples amplitudes are often determined
in a perturbative manner, based on the work of Raghavachari and co-workers
Raghavachari  _et al._ (1989). In their formalism, the perturbative triples
correction to the CCSD energy is defined as

$E^{(T)}=E^{[4]}_{T}+E^{[5]}_{ST}$ (10)

where

$\displaystyle E^{[4]}_{T}$ $\displaystyle=\left<T_{2}\mid[W,T_{3}]\right>$
(11) $\displaystyle E^{[5]}_{ST}$
$\displaystyle=\left<T_{1}\mid[W,T_{3}]\right>\ .$ (12)

$T_{1}$, $T_{2}$, and $T_{3}$ are known as the “cluster operators” and, in
second-quantization formalism, are defined as

$\displaystyle T_{1}$ $\displaystyle=t_{i}^{a}E_{ai}$ (13) $\displaystyle
T_{2}$ $\displaystyle=t_{ij}^{ab}E_{ai}E_{bj}$ (14) $\displaystyle T_{3}$
$\displaystyle=t_{ijk}^{abc}E_{ai}E_{bj}E_{ck}$ (15)

$E_{ai}$ represents the singlet, spin-adapted excitation operator, and is
defined as

$E_{ai}=a^{\dagger}_{a}a_{i}+{\bar{a}}^{\dagger}_{a}{\bar{a}}_{i}$ (16)

where the barred creation/annihilation operators refer to the beta spin
orbitals and nonbarred refer to the alpha spin orbitals.

The accuracy of the (T) method stems from a highly favorable error
cancellation between $E^{[4]}_{T}$ and $E^{[5]}_{ST}$. In restricted, single-
reference, closed-shell coupled cluster theory, one can write the equation for
the $(T)$ correction as Gyevi-Nagy  _et al._ (2019)

$E^{(T)}=\frac{1}{3}\frac{(4W_{ijk}^{abc}+W_{ijk}^{bca}+W_{ijk}^{cab})(V_{ijk}^{abc}-V_{ijk}^{cba})}{\epsilon_{i}+\epsilon_{j}+\epsilon_{k}-\epsilon_{a}-\epsilon_{b}-\epsilon_{c}}$
(17)

where

$W_{ijk}^{abc}=P_{L}\left[(ia|bd)t_{kj}^{cd}-(ia|jl)t_{kl}^{cb}\right]$ (18)

and

$V_{ijk}^{abc}=W_{ijk}^{abc}+P_{S}\left[t_{i}^{a}(jb|kc)\right]$ (19)

Following the formalism of Lesiuk Lesiuk (2020), we define $P_{L}$ and
$P_{S}$, or the “long” and “short” permutation operations as

$\displaystyle P_{L}(A_{ijk}^{abc})$
$\displaystyle=A_{ijk}^{abc}+A_{ikj}^{acb}+A_{jik}^{bac}+A_{jki}^{bca}+A_{kij}^{cab}+A_{kji}^{cba}$
(20) $\displaystyle P_{S}(A_{ijk}^{abc})$
$\displaystyle=A_{ijk}^{abc}+A_{jik}^{bac}+A_{kij}^{cab}$ (21)

The perturbative triples amplitude ($t_{ijk}^{abc}$), is defined as

$t_{ijk}^{abc}=\frac{W_{ijk}^{abc}}{\epsilon_{i}+\epsilon_{j}+\epsilon_{k}-\epsilon_{a}-\epsilon_{b}-\epsilon_{c}}$
(22)

Using the perturbative triples amplitude, as well as the permutational
symmetry of the Laplace denominator, one can rewrite Equation 17 as:

$E^{(T)}=t_{ijk}^{abc}\cdot(\frac{4}{3}V_{ijk}^{abc}-2V_{ijk}^{cba}+\frac{2}{3}V_{ijk}^{cab})$
(23)

We will use this equation when deriving the formulas for the THC-CCSD(T)
energy.

The cost of evaluating expression 23 scales as $\mathcal{O}(N^{6})$. However,
the cost of evaluating expression 18 scales as $\mathcal{O}(N^{7})$, leading
to an overall unfavorable $\mathcal{O}(N^{7})$ scaling of the CCSD(T) method.

###  II.3 Orthogonal Projectors

One crucial step of rank-reduced coupled cluster methods is the the formation
of the orthogonal projectors to reduce the dimensionality of the amplitudes,
as given in equations 1-4 and 8. There are a variety of methods that can be
used to compute orthogonal projectors. One such method for the CCSD doubles
amplitude is to form them from the definition of the MP2 $t_{ij}^{ab}$
amplitudes Parrish  _et al._ (2019).

$t_{ij}^{ab}=\frac{(ia|jb)}{\epsilon_{i}+\epsilon_{j}-\epsilon_{a}-\epsilon_{b}}$
(24)

Using density-fitting (DF) Dunlap  _et al._ (1979); Weigend  _et al._ (1998),
also known as resolution-of-the-identity (RI), or Cholesky Decomposition (CD)
Røeggen and Johansen (2008), the set of electron-repulsion integrals (ERIs) in
the molecular orbital (MO) basis $(ia|jb)$ can be written as follows DePrince
and Sherrill (2013):

$(ia|jb)\approx B^{Q}_{ia}B^{Q}_{jb}$ (25)

The energy denominator can be factored with a constant-sized index $w$ (with
growing molecular system size) through the Laplace denominator approach Häser
and Almlöf (1992)

$\frac{1}{\epsilon_{i}+\epsilon_{j}-\epsilon_{a}-\epsilon_{b}}=-D_{iw}D_{jw}D_{aw}D_{bw}$
(26)

Combining these techniques, and the following intermediates, as defined by
Parrish and co-workers Parrish  _et al._ (2019),

$\displaystyle L^{Qw}_{ia}$ $\displaystyle=D_{iw}D_{aw}B^{Q}_{ia}$ (27)
$\displaystyle M_{Pw,Qv}$ $\displaystyle=L^{Pw}_{ia}L^{Qv}_{ia}$ (28)

allows us to diagonalize $M$ and form the MP2 projector ($U^{V}_{ia}$) as

$\displaystyle M_{Pw,Qv}$ $\displaystyle=V^{V}_{Pw}\tau^{V}V^{V}_{Qv}$ (29)
$\displaystyle U^{V}_{ia}$
$\displaystyle=L^{Qw}_{ia}V^{V}_{Qw}\frac{1}{\sqrt{\tau^{V}}}\ .$ (30)

Note that the size of the index V can be truncated based on the magnitude of
the corresponding eigenvalue $\tau^{V}$. Even though the diagonalization of
$M$ is technically cubic-scaling, the size of the $w$ index can provide a
large prefactor. In the case of larger molecules, the size of the $V$ index is
often much smaller than the size of the $[Qw]$ index, and thus truncated
diagonalization approaches like the one given in reference 67 may be used.
Overall, this approach scales $\mathcal{O}(N^{4})$. Similarly, projectors can
be derived from MP3, albeit the equations are more complex Parrish  _et al._
(2019); Lesiuk (2022),

For triples amplitudes, we present two approaches devised by Lesiuk. In his
SVD-CCSDT algorithm Lesiuk (2020), he took guess $t_{ijk}^{abc}$ amplitudes,
such as from CC3, and applied either a TUCKER-3 decomposition (scaling
$\mathcal{O}(N^{8})$) or an iterative SVD approach (scaling
$\mathcal{O}(N^{6})$), yielding the form of equation 8.

In his RR-CCSD(T) paper, Lesiuk devised an $\mathcal{O}(N^{5})$ scheme to
compute projectors from the form of the perturbative triples amplitudes
(Equation 22), in a variant of HO-OI (Higher Order-Orthogonal Iteration)
Lesiuk (2022). The steps of the algorithm are as follows:

  * •

Start with the a guess of the triples projector $V_{ia}^{A}$. This can be done
naively by setting $V_{ia}^{A}=U_{ia}^{A}$ from the doubles amplitudes.

  * •

Evaluate $t_{ia,BC}$ from the current guess of the triples amplitudes, where

$t_{ia,BC}=t_{ijk}^{abc}V_{jb}^{B}V_{kc}^{C}$ (31)

By using the explicit expression for $t_{ijk}^{abc}$ and $W_{ijk}^{abc}$, this
can be evaluated in $\mathcal{O}(N^{5})$. The working equations are presented
in reference 54.

  * •

Compute the SVD of $t_{ia,BC}$, and take the largest $n_{proj}$ left singular
vectors as the next $V^{A}_{ia}$. This can be done in $\mathcal{O}(N^{5})$
time using a modified variant of truncated SVD, given in reference 67. In this
algorithm, we save the singular values of this step ($\sigma_{A}$), when we
perform the CP decomposition of the triples projector. Pseudocode for this
will be presented in Section IV.

  * •

Iterate until convergence. Convergence is defined when the difference between
the Frobenius norm of the rank-reduced triples amplitudes $t_{ABC}$, defined
as

$t_{ABC}=V^{A}_{ia}t_{ia,BC}$ (32)

between two successive iterations, falls below $10^{-5}$.

Since the source of the orthogonal projectors is not relevant to the scope of
this paper, we will only present results from computations utilizing the MP2
projector for the doubles amplitudes, and Lesiuk’s HO-OI approach for the
perturbative triples amplitudes.

###  II.4 Tensor Hypercontraction (THC)

Tensor hypercontraction (THC) can be viewed as a “double approximation,” where
two auxiliary indices are introduced to fit a high-dimensional tensor instead
of just one. The THC form of electron repulsion integrals is defined as
Hohenstein  _et al._ (2012a):

$(pq|rs)\approx x_{p}^{I}x_{q}^{I}Z^{IJ}x_{r}^{J}x_{s}^{J}$ (33)

This can be derived from the CP decomposition of $B^{Q}_{ia}$, (Equation 25)

$\displaystyle B^{Q}_{ia}\approx x_{i}^{I}x_{j}^{I}\eta^{QI}$ (34)
$\displaystyle Z^{IJ}=\eta^{QI}\eta^{QJ}\ .$ (35)

Similarly, the THC form of coupled-cluster amplitudes can be derived from the
tensor hypercontraction of the orthogonal projectors, given by, in the case of
the doubles projector:Hohenstein  _et al._ (2022)

$U^{V}_{ia}=y_{i}^{X}y_{a}^{X}\tau_{VX}\ .$ (36)

For the triples projector, it assumes a very similar form,

$V^{A}_{ia}=z_{i}^{X}z_{a}^{X}\theta_{AX}$ (37)

A PARAFAC/CANDENCOMP (CP) decomposition approach on $V^{A}_{ia}$ may be used.
This approach is not dependent on the source of the projectors, and any of the
projector building approaches from Section C may be used. Here we use the
variant of CP decomposition, first introduced by Hohenstein et. al. for the
doubles projector Hohenstein  _et al._ (2022), where the eigenvalues of the
doubles projector are in the CP decomposition, into the alternating least-
squares (ALS) iterations.

In our algorithm, for the decomposition of the triples amplitude, instead of
using the eigenvalues of the doubles projector, we use the singular values of
the $t_{ia,BC}$ intermediate ($\sigma_{A}$). The functional to minimize is
hence:

$L_{CP}=\sum_{ia}(\sigma_{A}[V^{A}_{ia}-z_{i}^{X}z_{a}^{X}\theta_{AX}])^{2}$
(38)

And the update rule for each intermediate is given as

$\displaystyle z_{i}^{X}$
$\displaystyle=\sum_{aA}\sigma_{A}^{2}V^{A}_{ia}\sum_{Y}z_{a}^{Y}\theta_{AY}[\sum_{b}z_{b}^{X}z_{b}^{Y}\sum_{B}\sigma_{B}^{2}\theta_{BX}\theta_{BY}]^{-1}$
(39) $\displaystyle z_{a}^{X}$
$\displaystyle=\sum_{iA}\sigma_{A}^{2}V^{A}_{ia}\sum_{Y}z_{i}^{Y}\theta_{AY}[\sum_{j}z_{j}^{X}z_{j}^{Y}\sum_{B}\sigma_{B}^{2}\theta_{BX}\theta_{BY}]^{-1}$
(40) $\displaystyle\theta_{AX}$
$\displaystyle=\sum_{ia}V^{A}_{ia}\sum_{Y}z_{i}^{Y}z_{a}^{Y}[\sum_{j}z_{j}^{X}z_{j}^{Y}\sum_{b}z_{b}^{X}z_{b}^{Y}]^{-1}$
(41)

Note that the update rule for $\theta$ is the same as in traditional CP
decomposition.

Since a CP decomposition does not exactly recreate the original projector, the
projectors lose their orthogonal property Hohenstein  _et al._ (2022).
Therefore, we have to re-create the projectors after the CP decomposition:

$\displaystyle S_{AB}$ $\displaystyle=V^{A}_{ia}V^{B}_{ia}$ (42)
$\displaystyle\theta_{AX}$ $\displaystyle=\theta_{BX}S_{AB}^{-1/2}$ (43)
$\displaystyle V^{A}_{ia}$ $\displaystyle=z_{i}^{X}z_{a}^{X}\theta_{AX}$ (44)

The $t_{ijk}^{abc}$ amplitudes can now be rewritten as, from equation 8:

$\displaystyle t_{ijk}^{abc}$
$\displaystyle=z_{i}^{X}z_{a}^{X}z_{j}^{Y}z_{b}^{Y}z_{k}^{Z}z_{c}^{Z}t_{XYZ}$
(45) $\displaystyle t_{XYZ}$
$\displaystyle=\theta_{AX}\theta_{BY}\theta_{CZ}t_{ABC}$ (46)

Recently, Hohenstein et. al. have devised an algorithm that takes advantage of
the THC form of the $t_{ij}^{ab}$ amplitudes to develop an $O(N^{4})$ scaling
implementation of CCSD Hohenstein  _et al._ (2022). In the next section, we
will show how to extend this to the (T) correction with the THC form of the
$t_{ijk}^{abc}$ amplitudes.

##  III Derivation of Working Equations

We first define a couple of intermediates. From Lesiuk Lesiuk (2022), we
define:

$D^{QV}_{jb}=(B^{Q}_{bd}U^{W}_{jd}-B^{Q}_{lj}U^{W}_{lb})T^{VW}$ (47)

Next, we define the following chain of intermediates from contracting the
polyadic vectors ($z_{i}^{X}$ and $z_{a}^{X}$) of the triples projector with
the the doubles projector, the DF/RI or CD decomposed ERIs, the D intermediate
from equation 47, as well as the $T_{1}$ amplitudes.

$\displaystyle\overline{U}^{VX}$ $\displaystyle=U^{V}_{ia}z_{i}^{X}z_{a}^{X}$
(48) $\displaystyle\widetilde{U}^{VXY}$
$\displaystyle=U^{V}_{ia}z_{i}^{X}z_{a}^{Y}$ (49)
$\displaystyle\overline{B}^{QX}$ $\displaystyle=B^{Q}_{ia}z_{i}^{X}z_{a}^{X}$
(50) $\displaystyle\widetilde{B}^{QXY}$
$\displaystyle=B^{Q}_{ia}z_{i}^{X}z_{a}^{Y}$ (51)
$\displaystyle\overline{D}^{QVX}$
$\displaystyle=D^{QV}_{ia}z_{i}^{X}z_{a}^{X}$ (52)
$\displaystyle\widetilde{D}^{QVXY}$
$\displaystyle=D^{QV}_{ia}z_{i}^{X}z_{a}^{Y}$ (53)
$\displaystyle\overline{t_{1}}^{X}$
$\displaystyle=t_{i}^{a}z_{i}^{X}z_{a}^{X}$ (54)
$\displaystyle\widetilde{t_{1}}^{XY}$
$\displaystyle=t_{i}^{a}z_{i}^{X}z_{a}^{Y}$ (55)

We then take Equation 23, Equation 19, Equation 45, and the previously defined
intermediates, to arrive at the THC form of the triples energy correction:

$\displaystyle
E^{(T)}\mathrel{+}=8\cdot\overline{U}^{VX}\overline{D}^{VQY}\overline{B}^{QZ}t_{XYZ}$
(56) $\displaystyle
E^{(T)}\mathrel{+}=4\cdot\overline{t_{1}}^{X}\overline{B}^{QY}\overline{B}^{QZ}t_{XYZ}$
(57) $\displaystyle
E^{(T)}\mathrel{-}=4\cdot\widetilde{U}^{VXZ}\overline{D}^{VQY}\widetilde{B}^{QZX}t_{XYZ}$
(58) $\displaystyle
E^{(T)}\mathrel{-}=4\cdot\widetilde{U}^{VXZ}\widetilde{D}^{VQZX}\overline{B}^{QY}t_{XYZ}$
(59) $\displaystyle
E^{(T)}\mathrel{-}=4\cdot\overline{U}^{VX}\widetilde{D}^{VQXZ}\widetilde{B}^{QZX}t_{XYZ}$
(60) $\displaystyle
E^{(T)}\mathrel{-}=4\cdot\widetilde{t_{1}}^{XZ}\widetilde{B}^{QZX}\overline{B}^{QY}t_{XYZ}$
(61) $\displaystyle
E^{(T)}\mathrel{-}=2\cdot\overline{t_{1}}^{Y}\widetilde{B}^{QXZ}\widetilde{B}^{QZX}t_{XYZ}$
(62) $\displaystyle
E^{(T)}\mathrel{+}=2\cdot\widetilde{U}^{VXZ}\widetilde{D}^{VQYX}\widetilde{B}^{QZY}t_{XYZ}$
(63) $\displaystyle
E^{(T)}\mathrel{+}=2\cdot\widetilde{U}^{VXZ}\widetilde{D}^{VQZY}\widetilde{B}^{QYX}t_{XYZ}$
(64) $\displaystyle
E^{(T)}\mathrel{+}=2\cdot\widetilde{t_{1}}^{XZ}\widetilde{B}^{QYX}\widetilde{B}^{QZY}t_{XYZ}$
(65)

Equations 56 and 57 correspond to the first term in equation 23, equations 58
through 62 the second term, and equations 63 to 65 the third term. All of the
contractions can be determined in $\mathcal{O}(N^{5})$ time or less.

##  IV Implementation Details

To aid future research and development, we present pseudocode for some of the
algorithms we use for the optimal contraction of intermediate terms to
evaluate the THC-CCSD(T) energy. We first present our non-iterative SVD
algorithm to factorize the $t_{ia,BC}$ intermediate, inspired by the truncated
SVD and diagonalization algorithms given in Ref. 67. In Algorithm 1, we
present a non-iterative truncated SVD algorithm to avoid the $O(N^{6})$
scaling of a traditional SVD of the $t_{ia,BC}$ intermediate. In Algorithms
2-4, we present suggested contraction orders, as well as tensor slicings, for
each term of the THC-CCSD(T) energy expression. We try to make the
contractions such that highly-efficient level 3 BLAS matrix multiplication
calls are utilized as much as possible. For each step of each algorithm, the
runtime is given, and if a level 3 BLAS matrix multiplication call is
possible, then the term (GEMM) is added. Additionally, the
$\widetilde{D}^{QVXY}$ intermediate is never fully built to help with memory
costs. The runtime of this algorithm is $\mathcal{O}(N^{5})$, with
$\mathcal{O}(N^{4})$ storage costs, the only quartic memory requirements
involve the storage of the $t_{ia,BC}$ and $D^{QV}_{jb}$ intermediates. It may
be possible to reduce the memory cost in future implementations of this
method, but that is beyond the scope of this paper.

Algorithm 1 Truncated SVD algorithm for $t_{ia,BC}$

$\Omega_{BC,X}=\textbf{random}(n_{proj}*n_{proj},n_{proj})$ $\triangleright$
$\mathcal{O}({N^{3}})$

$Y_{ia,X}=t_{ia,BC}\Omega_{BC,X}$ $\triangleright$ $\mathcal{O}(N^{5})$, GEMM

$Q_{ia,X},R_{X,Y}=\textbf{QR}(Y_{ia,X})$ $\triangleright$ $\mathcal{O}(N^{4})$

$t^{\prime}_{X,BC}=Q_{ia,X}t_{ia,BC}$ $\triangleright$ $\mathcal{O}(N^{5})$,
GEMM

$X_{XY}=t^{\prime}_{X,BC}t^{\prime}_{Y,BC}$ $\triangleright$
$\mathcal{O}(N^{4})$, GEMM

$V^{\prime}_{XY},\epsilon_{Y}=\textbf{diagonalize}(X_{XY})$ $\triangleright$
$\mathcal{O}(N^{3})$

$V^{A}_{ia}=Q_{ia,B}V^{\prime}_{BA}$ $\triangleright$ $\mathcal{O}(N^{4})$,
GEMM

$\sigma_{A}=\sqrt{\epsilon_{A}}$ $\triangleright$ $\mathcal{O}(N)$

$\textbf{return }V^{A}_{ia},\sigma_{A}$

Algorithm 2 $E_{1}^{(T)}$ Contractions (Equations 56 \- 57)

$A_{VYZ}=\overline{D}^{QVY}\overline{B}^{QZ}$ $\triangleright$
$\mathcal{O}(N^{4})$, GEMM

$B_{XYZ}=\overline{U}^{VX}\overline{A}^{VYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$, GEMM

$E^{(T)}\mathrel{+}=8\cdot B_{XYZ}t_{XYZ}$ $\triangleright$
$\mathcal{O}(N^{3})$

$C_{YZ}=\overline{B}^{QY}\overline{B}^{QZ}$ $\triangleright$
$\mathcal{O}(N^{3})$, GEMM

$D_{X}=C_{YZ}t_{XYZ}$ $\triangleright$ $\mathcal{O}(N^{3})$

$E^{(T)}\mathrel{+}=4\cdot D_{X}\overline{t_{1}}^{X}$ $\triangleright$
$\mathcal{O}(N)$

Algorithm 3 $E_{2}^{(T)}$ Contractions (Equations 58 \- 62)

for $V\texttt{ in }[0,n_{proj})$ do $\triangleright$ parallelize

$\widetilde{D}^{QVXY}=D^{QV}_{ia}z_{i}^{X}z_{a}^{Y}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM, built on the fly to save storage

$A^{YZX}=\overline{D}^{QVY}\widetilde{B}^{QZX}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM

$B^{YZX}=A^{YZX}\widetilde{U}^{VXZ}$ $\triangleright$ $\mathcal{O}(N^{4})$

$E^{(T)}\mathrel{-}=4\cdot B^{XYZ}t_{XYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$

$C^{ZXY}=\widetilde{D}^{QVZX}\overline{B}^{QY}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM

$D^{ZXY}=\widetilde{U}^{VXZ}C^{ZXY}$ $\triangleright$ $\mathcal{O}(N^{4})$

$E^{(T)}\mathrel{-}=4\cdot D^{XYZ}t_{XYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$

end for

for $Q\texttt{ in }[0,n_{aux})$ do $\triangleright$ parallelize

$\widetilde{D}^{QVXY}=D^{QV}_{ia}z_{i}^{X}z_{a}^{Y}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM, built on the fly to save storage

$F^{YXZ}=\overline{U}^{VY}\widetilde{D}^{QVXZ}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM

$G^{YXZ}=F^{YXZ}\widetilde{B}^{QZX}$ $\triangleright$ $\mathcal{O}(N^{4})$

$E^{(T)}\mathrel{-}=4\cdot G^{XYZ}t_{XYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$

end for

$H^{ZXY}=\widetilde{B}^{QZX}\overline{B}^{QY}$ $\triangleright$
$\mathcal{O}(N^{4})$, GEMM

$I^{ZXY}=H^{ZXY}\widetilde{t_{1}}^{XZ}$ $\triangleright$ $\mathcal{O}(N^{3})$

$E^{(T)}\mathrel{-}=4\cdot I^{XYZ}t_{XYZ}$ $\triangleright$
$\mathcal{O}(N^{3})$

$J^{XY}=t_{XYZ}\overline{t_{1}}^{Z}$ $\triangleright$ $\mathcal{O}(N^{3})$

$K^{XY}=\widetilde{B}^{QXY}\widetilde{B}^{QYX}$ $\triangleright$
$\mathcal{O}(N^{3})$

$E^{(T)}\mathrel{-}=2\cdot J^{XY}K^{XY}$ $\triangleright$ $\mathcal{O}(N^{2})$

Algorithm 4 $E_{3}^{(T)}$ Contractions (Equations 63 \- 65)

for $V\texttt{ in }[0,n_{proj})$ do $\triangleright$ parallelize

$\widetilde{D}^{QVXY}=D^{QV}_{ia}z_{i}^{X}z_{a}^{Y}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM, built on the fly to save storage

for $Y\texttt{ in }[0,n_{proj})$ do

$A^{XZ}=\widetilde{D}^{QVYX}\widetilde{B}^{QZY}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM

$B^{XZ}=A^{XZ}U^{VXZ}$ $\triangleright$ $\mathcal{O}(N^{4})$

$E^{(T)}\mathrel{+}=2\cdot B^{XZ}T^{XYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$

$C^{XZ}=\widetilde{D}^{QVZY}\widetilde{B}^{QYX}$ $\triangleright$
$\mathcal{O}(N^{5})$, GEMM

$D^{XZ}=C^{XZ}U^{VXZ}$ $\triangleright$ $\mathcal{O}(N^{4})$

$E^{(T)}\mathrel{+}=2\cdot D^{XZ}T^{XYZ}$ $\triangleright$
$\mathcal{O}(N^{4})$

end for

end for

for $Y\texttt{ in }[0,n_{proj})$ do $\triangleright$ parallelize

$F^{XZ}=\widetilde{B}^{QYX}\widetilde{B}^{QZY}$ $\triangleright$
$\mathcal{O}(N^{4})$, GEMM

$G^{XZ}=F^{XZ}\widetilde{t_{1}}^{XZ}$ $\triangleright$ $\mathcal{O}(N^{3})$

$E^{(T)}\mathrel{+}=2\cdot G^{XZ}T^{XYZ}$ $\triangleright$
$\mathcal{O}(N^{3})$

end for

The code is implemented in a developmental plugin version of the Psi4 Quantum
Chemistry code Smith  _et al._ (2020), following the completion of an exact
CCSD computation. Tensor contractions are performed with the help of the
EinsumsInCpp software (public on GitHub). The compressed doubles amplitudes
$T^{VW}$ used to build the triples projector are formed by transforming the
exact CCSD amplitudes from the preceding computation by the MP2 projector
amplitudes. This method is designed to be fully compatible and used with
Hohenstein’s THC-CCSD method Hohenstein  _et al._ (2022). Future studies of
using THC-CCSD(T) in conjunction with THC-CCSD is encouraged.

##  V Results

###  V.1 Conformation Energies

We first evaluate our new THC-CCSD(T) method on the CYCONF Wilke  _et al._
(2009) Goerigk and Grimme (2010) data set, a set containing 11 different
conformations of gaseous cysteine, with 10 corresponding conformation
energies, relative to the lowest conformer. We evaluate conformation energies
for each of the 10 conformations in CCSD, CCSD(T), and THC-CCSD(T), and for
each system, and we use the exact CCSD(T) conformation energy as the
reference. We do this using the cc-pVDZ and jun-cc-pVDZ Dunning correlation-
consistent basis sets Dunning (1989); Woon and Dunning (1993, 1994); Papajak
and Truhlar (2011). The basis set jun-cc-pVDZ consists of diffuse functions
added to all heavy atoms, except for the basis functions with the highest
angular momentum. For the THC-CCSD(T) computations, we set the eigenvalue
tolerance of the MP2 projector to be $10^{-4}$. In other words, the ranks
($n_{proj}$) of the doubles and triples projectors are determined from how
many eigenvalues of the MP2 $t_{ij}^{ab}$ amplitudes are greater than
$10^{-4}$, defined as $\tau$ from Equation 29 in our work. For these
computations, $n_{proj}$ is around 400, compared to the max possible rank of
2205 ($n_{occ}n_{virt}$) in the cc-pVDZ basis, yielding a compression ratio of
around 18%. Similarly, in the jun-cc-pVDZ basis, the ratio is 440/2793, which
is around 16%.

The summary statistics are presented in Table I, and the results for each
individual conformation are presented in Figure 1. In the table, for the THC-
CCSD(T) algorithms, the eigenvalue tolerance is given in parentheses. To
summarize the findings, THC-CCSD(T) consistently gives lower errors compared
to CCSD, for both basis sets, and the errors are on the order of less than
$0.1$ kcal/mol. It is further encouraging to note that the absolute energy
errors for these sets of computations hover around $0.3-0.4$ kcal/mol, such
that the evaluation of relative energies benefits from favorable error
cancellation. The error also does not significantly grow with the addition of
diffuse functions, from cc-pVDZ to jun-cc-pVDZ.

Table 1: Errors in conformation energy compared to the exact CCSD(T) reference (kcal/mol). The number in parenthesis is the eigenvalue tolerance used to determine projector rank. Test Set |  Mean Error |  MAE |  RMSE |  Std Dev
---|---|---|---|---
CCSD/cc-pVDZ |  $-0.343$ |  0.343 |  0.384 |  0.173
THC-CCSD(T)/cc-pVDZ ($10^{-4}$) |  $-0.072$ |  0.072 |  0.075 |  0.023
CCSD/jun-cc-pVDZ |  $-0.291$ |  0.291 |  0.323 |  0.141
THC-CCSD(T)/jun-cc-pVDZ ($10^{-4}$) |  $-0.076$ |  0.076 |  0.082 |  0.031
Figure 1: Errors in conformation energies for CCSD and THC-CCSD(T) evaluated
on the CYCONF data set, compared to the exact CCSD(T) reference, evaluated in
with the cc-pVDZ and jun-cc-pVDZ basis sets, with a $10^{-4}$ eigenvalue
tolerance.

###  V.2 Potential Energy Surface

We perform next, a potential energy surface scan on the benzene-HCN dimer
system (compound 19 from the on S22 data set Jurečka  _et al._ (2006)), with
the hydrogen atom of HCN pointing towards the $\pi$-bonds in the benzene. We
measured the energy of the system at five different inter-atomic distances,
relative to the equilibrium geometry, ranging from 0.9 to 2.0 times the
equilibrium geometry length, with the geometries coming from the S22x5 data
set Gráfová  _et al._ (2010). In Figure 2, we plot the shape of the potential
energy surface of the THC-CCSD(T) method at an eigenvalue tolerance of
$10^{-4}$, as well as using predetermined projector ranks of 400 and 500. For
all systems, an eigenvalue tolerance of $10^{-4}$ corresponds to a projector
rank between 420-430. All THC-CCSD(T) computations better capture the
potential energy surface than the reference CCSD computations, with the
computations with the predetermined projector ranks better capturing the shape
of the surface than the one with a set eigenvalue tolerance. The THC-CCSD(T)
potential energy surface with $n_{proj}$ set to 500 exactly matches the
CCSD(T) potential energy surface, for practical purposes, with a max error of
0.027 kcal/mol, and a RMSE of 0.014 kcal/mol. The shape of the potential
energy surface, for each method, is shown in Figure 2, while the error
statistics are presented in Table II. The errors are especially encouraging
for the case of $n_{proj}$ set to 500, as the absolute energy error of each
system compared to CCSD(T) hover around 0.4 kcal/mol.

Table 2: Errors in relative energies compared to the exact CCSD(T) reference (kcal/mol), for a reference CCSD computation, as well as THC-CCSD(T) computations with varying parameters. Test Set |  Mean Error |  MAE |  RMSE |  Std Dev
---|---|---|---|---
CCSD |  $-0.138$ |  0.200 |  0.236 |  0.191
THC-CCSD(T), tol = $10^{-4}$ |  $-0.100$ |  0.103 |  0.132 |  0.086
THC-CCSD(T), $n_{proj}=400$ |  $-0.098$ |  0.098 |  0.128 |  0.082
THC-CCSD(T), $n_{proj}=500$ |  $-0.001$ |  0.010 |  0.014 |  0.014

Figure 2: The relative energies of benzene-HCN dimer (S22 system 19) evaluated
with each method at five different dimer separation distances relative to the
equilibrium geometry.

###  V.3 Rank Convergence

Next, to demonstrate the convergence of the THC-CCSD(T) method, compared to
the exact CCSD(T) energy, we ran a series of computations of the water dimer
from the S22 set Jurečka  _et al._ (2006), at eigenvalue tolerances from
$10^{-3}$ to $10^{-11}$. An eigenvalue tolerance of $10^{-11}$ corresponds to
no rank compression for this system. The errors with respect to eigenvalue
tolerance and compression ranks are plotted in Figure 3, and it is encouraging
to see the errors decrease smoothly to the true CCSD(T) energy, within the
DF/RI approximation of the ERIs. We attribute the “kink” in the graph from
$10^{-4}$ to $10^{-6}$ as an artifact of the CP decomposition of the triples
projector, with the CP error increasing slightly between the projector ranks
of 122 - 156, before going back down. This artifact is well known on studies
of the CP decomposition algorithm Kolda and Bader (2009), where medium CP
decomposition ranks suffer larger losses in accuracy compared to small or
large ranks. Further studies and work are encouraged to look for ways to
mitigate this phenomenon in the context of decomposing CC amplitudes.

Figure 3: The convergence of the absolute energy of a water dimer system
(S22), with respect to eigenvalue tolerance and rank.

###  V.4 Scaling

To establish the $\mathcal{O}(N^{5})$ scaling of the THC-CCSD(T) method, it
must be shown that the projector rank, or $n_{proj}$ must scale linearly with
respect to system size. Hohenstein and Parrish have previously established the
linear scaling of $n_{proj}$ for doubles amplitudes in their previous work
Parrish  _et al._ (2019); Hohenstein  _et al._ (2022). However, to verify this
in our algorithm, we must show that the error does not grow more than linearly
with linear increases in system size. Below, we present THC-CCSD(T)
computations on systematically larger systems of waterclusters and linear
alkanes, from 1-8 heavy atoms, in the cc-pVDZ and jun-cc-pVDZ basis sets,
evaluated at an eigenvalue tolerance of $10^{-4}$. As shown in Figures 4-7,
sub-linear to linear error growth are shown, with respect to projector rank
and system size, with virtually no loss in accuracy from cc-pVDZ to jun-cc-
pVDZ in both systems.

Figure 4: The growth of the absolute energy error, and projector rank, in a
growing series of water clusters $(H_{2}O)_{n}$, cc-pVDZ basis. Figure 5: The
growth of the absolute energy error, and projector rank, in a growing series
of water clusters $(H_{2}O)_{n}$, jun-cc-pVDZ basis. Figure 6: The growth of
the absolute energy error, and projector rank, in a growing series of linear
alkanes $(C_{n}H_{2n+2})$, cc-pVDZ basis. Figure 7: The growth of the absolute
energy error, and projector rank, in a growing series of linear alkanes
$(C_{n}H_{2n+2})$, jun-cc-pVDZ basis.

##  VI Conclusions

In this paper, we present the working equations for the THC-CCSD(T) method, a
$\mathcal{O}(N^{5})$ scaling approximation to CCSD(T), that allows for
systematic control of errors. In our pilot implementation, we show the errors
are controllable to the point of maintaining chemical accuracy of less than
0.1 kcal/mol for relative energies, and 1 mEh for absolute energies, while
maintaining size extensivity. We also showed that the method yields continuous
potential energy surfaces that closely matches the CCSD(T) surfaces with
sufficient projector rank. In the future, we hope to consider ways to improve
the errors of the method at a given eigenvalue tolerance, such as through
using other sources for the orthogonal projector. We would also like to look
into alternative approaches to the THC factorization of orthogonal projectors.
Though a CP decomposition is generally applicable, and relatively easy to
implement, it does not assume any underlying form about the amplitudes. One
avenue is the extension of the quadrature-based approach of Parrish,
Hohenstein, Martinez, and Sherrill with Least-Squares Tensor Hypercontraction
(LS-THC) to the triples amplitudes Parrish  _et al._ (2012); Hohenstein  _et
al._ (2012b).

## Acknowledgements

The authors gratefully acknowledge financial support from the U.S. Department
of Energy, Basic Energy Sciences Division, Computational and Theoretical
Chemistry (CTC) Grant DE-SC0018164.

## Data Availability

The data that supports the findings of this study are available with the
article and its supplementary material.

## References

  * Crawford and H. F. Schaefer (2007) T. D. Crawford and H. F. Schaefer, Reviews in Computational Chemistry , 33 (2007).
  * Bartlett and Musial (2007) R. J. Bartlett and M. Musial, Reviews of Modern Physics 79 (2007).
  * Cramer (2002) C. J. Cramer,  _Essentials of Computational Chemistry_ (2002) pp. 191–232.
  * Sherrill and H. F. Schaefer (1999) C. D. Sherrill and H. F. Schaefer, Advances in Quantum Chemistry 34, 143 (1999).
  * Raghavachari  _et al._ (1989) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chemical Physics Letters 157, 479 (1989).
  * Hirata (2003) S. Hirata, The Journal of Physical Chemistry A 107, 9887 (2003).
  * Auer  _et al._ (2006) A. A. Auer, G. Baumgartner, D. E. Bernholdt, A. Bibireata, V. Choppella, D. Cociorva, X. Gao, R. Harrison, S. Krishnamoorthy, S. Krishnan, C.-C. Lam, Q. Lu, M. Nooijen, R. Pitzer, J. Ramanujam, P. Sadayappan, and A. Sibiryakov, Molecular Physics 104, 211 (2006).
  * Janowski  _et al._ (2007) T. Janowski, A. R. Ford, and P. Pulay, Journal of Chemical Theory and Computation 3, 1368 (2007).
  * Janowski and Pulay (2008) T. Janowski and P. Pulay, Journal of Chemical Theory and Computation 4, 1585 (2008).
  * van Dam  _et al._ (2011) H. van Dam, W. de Jong, E. Bylaska, N. Govind, K. Kowalski, T. Straatsma, and M. Valiev, Wiley Interdisciplinary Reviews: Computational Molecular Science 1, 888 (2011).
  * Deumens  _et al._ (2011) E. Deumens, V. F. Lotrich, A. Perera, M. J. Ponton, B. A. Sanders, and R. J. Bartlett, Wiley Interdisciplinary Reviews: Computational Molecular Science 1, 895 (2011).
  * Kobayashi and Rendell (1997) R. Kobayashi and A. P. Rendell, Chemical Physics Letters 265, 1 (1997).
  * Anisimov  _et al._ (2014) V. M. Anisimov, G. H. Bauer, K. Chadalavada, R. M. Olson, J. W. Glenski, W. T. C. Kramer, E. Aprà, and K. Kowalski, Journal of Chemical Theory and Computation 10, 4307 (2014).
  * Solomonik  _et al._ (2014) E. Solomonik, D. Matthews, J. R. Hammond, J. F. Stanton, and J. Demmel, Journal of Parallel and Distributed Computing 74, 3176 (2014).
  * Peng  _et al._ (2016) C. Peng, J. A. Calvin, F. Pavošević, J. Zhang, and E. F. Valeev, The Journal of Physical Chemistry A 120, 10231 (2016).
  * Lyakh (2019) D. I. Lyakh, International Journal of Quantum Chemistry 119, e25926 (2019).
  * Gyevi-Nagy  _et al._ (2020) L. Gyevi-Nagy, M. Kállay, and P. R. Nagy, Journal of Chemical Theory and Computation 16, 366 (2020).
  * Peng  _et al._ (2020) C. Peng, C. A. Lewis, X. Wang, M. C. Clement, K. Pierce, V. Rishi, F. Pavošević, S. Slattery, J. Zhang, N. Teke, A. Kumar, C. Masteran, A. Asadchev, J. A. Calvin, and E. F. Valeev, The Journal of Chemical Physics 153, 044120 (2020).
  * Datta and Gordon (2021) D. Datta and M. S. Gordon, Journal of Chemical Theory and Computation 17, 4799 (2021).
  * Gyevi-Nagy  _et al._ (2021) L. Gyevi-Nagy, M. Kállay, and P. R. Nagy, Journal of Chemical Theory and Computation 17, 860 (2021).
  * Kowalski  _et al._ (2021) K. Kowalski, R. Bair, N. P. Bauman, J. S. Boschen, E. J. Bylaska, J. Daily, W. A. de Jong, T. Dunning, N. Govind, R. J. Harrison, M. Keçeli, K. Keipert, S. Krishnamoorthy, S. Kumar, E. Mutlu, B. Palmer, A. Panyala, B. Peng, R. M. Richard, T. P. Straatsma, P. Sushko, E. F. Valeev, M. Valiev, H. J. J. van Dam, J. M. Waldrop, D. B. Williams-Young, C. Yang, M. Zalewski, and T. L. Windus, Chemical Reviews 121, 4962 (2021).
  * Calvin  _et al._ (2021) J. A. Calvin, C. Peng, V. Rishi, A. Kumar, and E. F. Valeev, Chemical Reviews 121, 1203 (2021).
  * Seritan  _et al._ (2020) S. Seritan, C. Bannwarth, B. S. Fales, E. G. Hohenstein, S. I. L. Kokkila-Schumacher, N. Luehr, J. W. Snyder, C. Song, A. V. Titov, I. S. Ufimtsev, and T. J. Martínez, The Journal of Chemical Physics 152, 224110 (2020).
  * Wang  _et al._ (2020) Z. Wang, M. Guo, and F. Wang, Physical Chemistry Chemical Physics 22, 25103 (2020).
  * Peng  _et al._ (2019) C. Peng, J. A. Calvin, and E. F. Valeev, International Journal of Quantum Chemistry 119 (2019), 10.1002/qua.25894.
  * Kaliman and Krylov (2017) I. A. Kaliman and A. I. Krylov, Journal of Computational Chemistry 38, 842 (2017).
  * DePrince  _et al._ (2014) A. E. DePrince, M. R. Kennedy, B. G. Sumpter, and C. D. Sherrill, Molecular Physics 112, 844 (2014).
  * Ma  _et al._ (2011) W. Ma, S. Krishnamoorthy, O. Villa, and K. Kowalski, Journal of Chemical Theory and Computation 7, 1316 (2011).
  * DePrince and Hammond (2011) A. E. DePrince and J. R. Hammond, Journal of Chemical Theory and Computation 7, 1287 (2011).
  * Møller and Plesset (1934) C. Møller and M. S. Plesset, Physical Review 46, 618 (1934).
  * Cremer (2011) D. Cremer, Wiley Interdisciplinary Reviews: Computational Molecular Science 1, 509 (2011).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Physical Review 136, B864 (1964).
  * Kohn and Sham (1965) W. Kohn and L. J. Sham, Physical Review 140, A1133 (1965).
  * Maurer  _et al._ (2014) S. A. Maurer, L. Clin, and C. Ochsenfeld, The Journal of Chemical Physics 140, 224112 (2014).
  * Dawson  _et al._ (2022) W. Dawson, A. Degomme, M. Stella, T. Nakajima, L. E. Ratcliff, and L. Genovese, WIREs Computational Molecular Science 12 (2022), 10.1002/wcms.1574.
  * Li  _et al._ (2002) S. Li, J. Ma, and Y. Jiang, Journal of Computational Chemistry 23, 237 (2002).
  * Li  _et al._ (2006) S. Li, J. Shen, W. Li, and Y. Jiang, The Journal of Chemical Physics 125, 074109 (2006).
  * Li  _et al._ (2009) W. Li, P. Piecuch, J. R. Gour, and S. Li, The Journal of Chemical Physics 131, 114109 (2009).
  * Neese  _et al._ (2009) F. Neese, F. Wennmohs, and A. Hansen, The Journal of Chemical Physics 130, 114108 (2009).
  * Li and Piecuch (2010a) W. Li and P. Piecuch, The Journal of Physical Chemistry A 114, 6721 (2010a).
  * Li and Piecuch (2010b) W. Li and P. Piecuch, The Journal of Physical Chemistry A 114, 8644 (2010b).
  * Rolik and Kállay (2011) Z. Rolik and M. Kállay, The Journal of Chemical Physics 135, 104111 (2011).
  * Rolik  _et al._ (2013) Z. Rolik, L. Szegedy, I. Ladjánszki, B. Ladóczki, and M. Kállay, The Journal of Chemical Physics 139, 094105 (2013).
  * Riplinger and Neese (2013) C. Riplinger and F. Neese, The Journal of Chemical Physics 138, 034106 (2013).
  * Riplinger  _et al._ (2013) C. Riplinger, B. Sandhoefer, A. Hansen, and F. Neese, The Journal of Chemical Physics 139, 134101 (2013).
  * Liakos  _et al._ (2015) D. G. Liakos, M. Sparta, M. K. Kesharwani, J. M. L. Martin, and F. Neese, Journal of Chemical Theory and Computation 11, 1525 (2015).
  * Schwilk  _et al._ (2017) M. Schwilk, Q. Ma, C. Köppl, and H.-J. Werner, Journal of Chemical Theory and Computation 13, 3650 (2017).
  * Pinski  _et al._ (2015) P. Pinski, C. Riplinger, E. F. Valeev, and F. Neese, The Journal of Chemical Physics 143, 034108 (2015).
  * Riplinger  _et al._ (2016) C. Riplinger, P. Pinski, U. Becker, E. F. Valeev, and F. Neese, The Journal of Chemical Physics 144, 024109 (2016).
  * Parrish  _et al._ (2019) R. M. Parrish, Y. Zhao, E. G. Hohenstein, and T. J. Martínez, The Journal of Chemical Physics 150, 164118 (2019), https://doi.org/10.1063/1.5092505 .
  * Kolda and Bader (2009) T. G. Kolda and B. W. Bader, SIAM Review 51, 455 (2009).
  * Hohenstein  _et al._ (2022) E. G. Hohenstein, B. S. Fales, R. M. Parrish, and T. J. Martínez, The Journal of Chemical Physics 156, 054102 (2022), https://doi.org/10.1063/5.0077770 .
  * Lesiuk (2020) M. Lesiuk, Journal of Chemical Theory and Computation 16, 453 (2020), pMID: 31715103, https://doi.org/10.1021/acs.jctc.9b00985 .
  * Lesiuk (2022) M. Lesiuk, The Journal of Chemical Physics 156, 064103 (2022), https://doi.org/10.1063/5.0071916 .
  * Riley  _et al._ (2010) K. E. Riley, M. Pitoňák, P. Jurečka, and P. Hobza, Chemical Reviews 110, 5023 (2010).
  * Karton  _et al._ (2006) A. Karton, E. Rabinovich, J. M. L. Martin, and B. Ruscic, The Journal of Chemical Physics 125, 144108 (2006).
  * Tajti  _et al._ (2004) A. Tajti, P. G. Szalay, A. G. Császár, M. Kállay, J. Gauss, E. F. Valeev, B. A. Flowers, J. Vázquez, and J. F. Stanton, The Journal of Chemical Physics 121, 11599 (2004).
  * Bak  _et al._ (2000) K. L. Bak, P. Jørgensen, J. Olsen, T. Helgaker, and W. Klopper, The Journal of Chemical Physics 112, 9229 (2000).
  * Hopkins and Tschumper (2004) B. W. Hopkins and G. S. Tschumper, The Journal of Physical Chemistry A 108, 2941 (2004).
  * Bartlett  _et al._ (1990) R. J. Bartlett, J. Watts, S. Kucharski, and J. Noga, Chemical Physics Letters 165, 513 (1990).
  * Gyevi-Nagy  _et al._ (2019) L. Gyevi-Nagy, M. Kállay, and P. R. Nagy, Journal of Chemical Theory and Computation 16, 366 (2019).
  * Dunlap  _et al._ (1979) B. I. Dunlap, J. W. D. Connolly, and J. R. Sabin, The Journal of Chemical Physics 71, 3396 (1979).
  * Weigend  _et al._ (1998) F. Weigend, M. Häser, H. Patzelt, and R. Ahlrichs, Chemical Physics Letters 294, 143 (1998).
  * Røeggen and Johansen (2008) I. Røeggen and T. Johansen, The Journal of Chemical Physics 128, 194107 (2008).
  * DePrince and Sherrill (2013) A. E. DePrince and C. D. Sherrill, Journal of Chemical Theory and Computation 9, 2687 (2013).
  * Häser and Almlöf (1992) M. Häser and J. Almlöf, The Journal of Chemical Physics 96, 489 (1992), https://doi.org/10.1063/1.462485 .
  * Halko  _et al._ (2011) N. Halko, P. G. Martinsson, and J. A. Tropp, SIAM Review 53, 217 (2011), https://doi.org/10.1137/090771806 .
  * Hohenstein  _et al._ (2012a) E. G. Hohenstein, R. M. Parrish, and T. J. Martínez, The Journal of Chemical Physics 137, 044103 (2012a), https://doi.org/10.1063/1.4732310 .
  * Smith  _et al._ (2020) D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. D. Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrince, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, The Journal of Chemical Physics 152, 184108 (2020).
  * Wilke  _et al._ (2009) J. J. Wilke, M. C. Lind, H. F. Schaefer, A. G. Császár, and W. D. Allen, Journal of Chemical Theory and Computation 5, 1511 (2009).
  * Goerigk and Grimme (2010) L. Goerigk and S. Grimme, Journal of Chemical Theory and Computation 6, 107 (2010).
  * Dunning (1989) T. H. Dunning, The Journal of Chemical Physics 90, 1007 (1989).
  * Woon and Dunning (1993) D. E. Woon and T. H. Dunning, The Journal of Chemical Physics 98, 1358 (1993).
  * Woon and Dunning (1994) D. E. Woon and T. H. Dunning, The Journal of Chemical Physics 100, 2975 (1994).
  * Papajak and Truhlar (2011) E. Papajak and D. G. Truhlar, Journal of Chemical Theory and Computation 7, 10 (2011).
  * Jurečka  _et al._ (2006) P. Jurečka, J. Šponer, J. Černý, and P. Hobza, Phys. Chem. Chem. Phys. 8, 1985 (2006).
  * Gráfová  _et al._ (2010) L. Gráfová, M. Pitoňák, J. Řezáč, and P. Hobza, Journal of Chemical Theory and Computation 6, 2365 (2010).
  * Parrish  _et al._ (2012) R. M. Parrish, E. G. Hohenstein, T. J. Martínez, and C. D. Sherrill, The Journal of Chemical Physics 137, 224106 (2012).
  * Hohenstein  _et al._ (2012b) E. G. Hohenstein, R. M. Parrish, C. D. Sherrill, and T. J. Martínez, The Journal of Chemical Physics 137, 221101 (2012b).