# Externally corrected CCSD with renormalized perturbative triples
(R-ecCCSD(T)) and density matrix renormalization group and selected
configuration interaction external sources

Seunghoon Lee  Division of Chemistry and Chemical Engineering, California
Institute of Technology, Pasadena, California 91125, USA<EMAIL_ADDRESS>Huanchen Zhai  Division of Chemistry and Chemical Engineering, California
Institute of Technology, Pasadena, California 91125, USA  Sandeep Sharma
Department of Chemistry, The University of Colorado at Boulder, Boulder,
Colorado 80302, USA<EMAIL_ADDRESS>C. J. Umrigar  Laboratory of
Atomic and Solid State Physics, Cornell University, Ithaca, New York 14853,
USA<EMAIL_ADDRESS>Garnet Kin-Lic Chan  Division of Chemistry and
Chemical Engineering, California Institute of Technology, Pasadena, California
91125, USA<EMAIL_ADDRESS>
###### Abstract

We investigate the renormalized perturbative triples correction together with
the externally corrected coupled-cluster singles and doubles (ecCCSD) method.
We take the density matrix renormalization group (DMRG) and heatbath CI (HCI)
as external sources for the ecCCSD equations. The accuracy is assessed for the
potential energy surfaces of H2O, N2, and F2. We find that the triples
correction significantly improves on ecCCSD and we do not see any instability
of the renormalized triples with respect to dissociation. We explore how to
balance the cost of computing the external source amplitudes with respect to
the accuracy of the subsequent CC calculation. In this context, we find that
very approximate wavefunctions (and their large amplitudes) serve as an
efficient and accurate external source. Finally, we characterize the domain of
correlation treatable using the externally corrected method and renormalized
triples combination studied in this work via a well-known wavefunction
diagnostic.

##  1 Introduction

Electronic structure methods that can efficiently handle both static and
dynamic correlation remain an important area of investigation. Because there
is a wide spectrum of strongly correlated problems, ranging from mildly
"quasi"-degenerate scenarios (e.g. in the electronic structure of diradicals
1, 2) to extensively near-degenerate problems (e.g. in the electronic
structure of multi-centre transition metal clusters 3, 4, 5) a variety of
theoretical strategies have been proposed.

For highly-degenerate problems, it is common to combine dynamic correlation
methods with an explicit treatment of a multi-reference state. The
representation of the multi-reference state can range from an exact complete
active space (CAS) representation 6, 7, 8 (for small numbers of orbitals), to
density matrix renormalization group (DMRG) 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23 (e.g. when almost all orbitals are degenerate), to
selected configuration interaction 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39, 40, 41, 42, 43, 44 and Monte Carlo approximations 45, 46,
47, 48, 49 for intermediate cases. On top of these, various flavours of
perturbation theory 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62,
configuration interaction 63, 64, 65, 66, and exponential approximations 67,
68, 69, 58 have been explored. However, the combination of dynamic correlation
with multi-reference representations is not straightforward and usually leads
to added conceptual, implementation, and computational complexity.

For quasi-degenerate problems, an alternative strategy can be used, which
incorporates a limited amount of static correlation on top of an existing SR
method. This has been particularly popular in conjunction with SR coupled
cluster methods. 70, 71, 72, 73, 74, 75 Some examples include variants of
tailored coupled cluster 76, 77, 78, 79, 80, 81 and externally corrected
coupled cluster methods. 82, 83, 84, 85 The computational cost of such SR
static correlation methods is often lower than that of true MR dynamic
correlation methods and thus large active spaces become affordable. However,
the approximations are limited to problems with only a modest amount of
degeneracy.

In this work, we will focus on quasi-degenerate problems, and in particular,
we will investigate the externally corrected coupled cluster method. 82, 83,
84, 85 This extracts static correlation from a MR method by using the MR
wavefunction as an "external" source of higher order coupled cluster
amplitudes. For example, in the ecCCSD approximation, the $T_{3}$ and $T_{4}$
amplitudes are extracted from the external source and a new set of $T_{1}$ and
$T_{2}$ amplitudes are computed in their presence. Should the $T_{3}$ and
$T_{4}$ amplitudes be exact, then the $T_{1}$ and $T_{2}$ amplitudes and the
energy will be exact. A different, spiritually related, approximation is
tailored CCSD 76, 77, 78, which has been of renewed interest of late 79, 80,
81. Here, instead of higher order cluster amplitudes, the large (active space)
$T_{1}$ and $T_{2}$ amplitudes are fixed from an external source.

The ecCCSD method has a long history and external sources, ranging from
unrestricted Hartree-Fock 86, 87, 88, 89 to CASSCF and CASCI 90, 91, 92, 93,
94 and, most recently, full configuration interaction quantum Monte Carlo95,
have been used. One of the more successful applications of ecCCSD is the
reduced multireference (RMR) CCSD method96, 97, which uses a MRCISD wave
function as the external source. RMR CCSD(T), which incorporates some of the
residual dynamic correlation through perturbative triples, has also been
studied.98 Despite the promising performance of RMRCCSD(T) in several studies
99, 100, 101, 102, 103, 104, 105, 106, 107, it suffers from two main
limitations. First, conventional MRCISD can only be applied for modest sizes
of active spaces (typically, up to about $16$ orbitals as limited by the exact
CAS treatment). Second, the (T) correction, although not divergent like its
single-reference counterpart, still overcorrects the dynamical correlation in
the bond-stretched region.108

In this work, we make two modifications to ecCCSD to overcome and ameliorate
the above limitations. First, we utilize variational DMRG and HCI wave
functions as external sources for ecCCSD. This allows for the use of larger
quasidegenerate active spaces, of the size typically treated by DMRG and HCI.
Second, we explore the renormalized perturbative triples correction. This has
been shown in the single reference setting to ameliorate the overcorrection of
standard perturbative triples 109, without affecting the computational
scaling.

We describe the use of DMRG and HCI as external sources for ecCCSD in Section
2.1. It is possible to create a near-exact method by using a near-exact
external source. Since near-exact external sources can be computed by DMRG and
HCI for the small systems we employ as test cases in this paper, the critical
question is not simply the accuracy of the method, but the balance between
cost and accuracy. To this end, we explore a variety of approximate treatments
of the external DMRG and HCI sources as discussed in Sec. 2.2. The various
triples approximations for ecCCSD are discussed in Section 2.3. Computational
details are provided in Section 3 and the accuracy of the renormalized
perturbative triples correction is assessed for three potential energy
surfaces (PESs) in Sections 4.1-4.4. We characterize the range of quasi-
degenerate correlations captured in this work in Section 4.6. Finally we
discuss the limitations of this method in Section 4.7 and summarize our
findings in Section 5.

##  2 Theory

In this work, we use a reference configuration $|\Phi\rangle$ with the same
occupancy as the Hartree-Fock (HF) determinant. It is useful to define a
projection operator onto the space of $k$-tuply excited configurations
relative to the reference; we denote this $Q_{k}$. The external source is used
to provide an important subset of the triply and quadruply excited
configurations; the projector onto this subset is denoted $Q_{k}^{ec}$, and
its complement is $Q_{k}^{c}$. Thus

$\displaystyle Q_{k}$ $\displaystyle=Q^{\rm ec}_{k}+Q^{\rm c}_{k},\ \ $
$\displaystyle k$ $\displaystyle=3,4,$ (1)

###  2.1 Externally corrected CCSD with DMRG and HCI wave functions

In ecCCSD, the coupled cluster operator $T$ is given by

$\displaystyle T=T_{1}+T_{2}+Q^{\rm ec}_{3}T_{3}+Q^{\rm ec}_{4}T_{4}.$ (2)

Here, $T_{n},\ n=1,2,3,4,$ are the $n$-fold cluster operators. In this work,
we extract $Q^{\rm ec}_{k}T_{k}$, $k=3,4$ from the DMRG or the HCI variational
wave function. For the DMRG wave function, the triply and quadruply excited
configurations $D_{p}$ that define $Q^{\rm ec}_{k}$, $k=3,4$ are chosen to be
those where the magnitude of the CI coefficient $c_{p}$ is above a threshold,
i.e.

$\displaystyle|c_{p}|>s\sqrt{\omega},$ (3)

where $s$ is an arbitrary scaling factor and $\omega$ is the largest discarded
weight of the density matrix at the maximum bond dimension in two-dot DMRG
sweeps (carried out without noise). An efficient algorithm to convert a
matrix-product state to CI coefficients above a given threshold is described
in Appendix A. For the HCI variational wave function, $D_{p}$ is included in
the projectors $Q^{\rm ec}_{k}$, $k=3,4$ using the heatbath algorithm with a
threshold $\epsilon$,42 i.e., it is included if

$\displaystyle|\langle D_{p}|H|D_{q}\rangle c_{q}|>\epsilon,$ (4)

for at least one determinant $D_{q}$ which is already in the variational
space. The extracted CI coefficients are then converted into cluster
amplitudes.

The $T_{1}$ and $T_{2}$ amplitudes are obtained by solving the ecCCSD
equations using fixed $Q^{\rm ec}_{3}T_{3}$ and $Q^{\rm ec}_{4}T_{4}$,

$0=(Q_{1}+Q_{2})(H_{N}e^{T_{1}+T_{2}+Q^{\rm ec}_{3}T_{3}+Q^{\rm
ec}_{4}T_{4}})_{C}|\Phi\rangle,$ (5)

where $H_{N}$ is the Hamiltonian in normal-ordered form, and the subscript $C$
denotes the connected part of the corresponding operator expression. With the
relaxed $T_{1}$ and $T_{2}$, the ecCCSD correlation energy of the ground state
is obtained as

$\displaystyle\Delta E^{\rm
ecCCSD}_{0}=\langle\Phi|(H_{N}e^{T_{1}+T_{2}+Q^{\rm ec}_{3}T_{3}+Q^{\rm
ec}_{4}T_{4}})_{C}|\Phi\rangle.$ (6)

###  2.2 Approximations in the external source

While DMRG and HCI can, in small molecules, be a source of nearly exact
$T_{3}$ and $T_{4}$ amplitudes even in the full orbital space, this provides
no computational advantage as such calculations are more expensive than the
subsequent coupled cluster calculation. Consequently, it is important to
balance the cost of the external source calculation and that of the subsequent
coupled cluster calculation by making approximations in the external source.
This introduces the additional complication that one must ensure that errors
introduced into the external source do not lead to unacceptable errors in the
final coupled cluster calculation.

In this work, we consider six different types of approximate external sources
with different sizes of active spaces and different values of parameters
summarized in Table 1. Type I uses CASSCF-like external sources in minimal
active spaces. Since the minimal active space of F2, i.e., two electrons in
two orbitals (2e,2o), does not contain $T_{3}$ and $T_{4}$, we perform minimal
active space ecCC calculations for only H2O and N2 (with (4e,4o) and (6e,6o),
respectively). These provide amplitudes that are very close to the exact
CASSCF amplitudes. However, these amplitudes lack the relaxation that comes
from allowing excitations within a larger space of orbitals, and of course the
amplitudes outside the minimal active space are completely absent. Type III
uses CASCI-like external sources (the orbitals are not optimized to save
computer time) in larger active spaces ((8e,18o), (10e,16o), and (14e,16o) for
H2O, N2, and F2, respectively).

In either case, one can potentially introduce bad external amplitudes if the
effect of relaxation on the amplitude upon going to a larger space is large
relative to the size of the amplitude (e.g. changes its sign). Thus, we study
also Type II and Type IV external sources which are similar to Type I and Type
II external sources respectively except that they employ an additional
threshold to screen out all except the largest $T_{3}$ and $T_{4}$ amplitudes.
The absolute values of the $T_{3}$ and $T_{4}$ elements at the most stretched
geometry of each molecule are sorted in a single large vector. The norm of the
vector is computed. Only the largest elements of the vector are retained such
that the resulting norm is more than 80% of the norm of the full vector. Along
PESs of each molecule, we used the same set of elements of $T_{3}$ and $T_{4}$
(but with the appropriate values for each geometry) as the external sources,
to maintain the smoothness of the PESs.

As discussed in Section 4.4, the type-III and type-IV sources improve upon the
PESs obtained from the type-I and type-II sources, but the DMRG calculations
to obtain the sources incur a higher computational cost than the subsequent CC
calculations. To reduce the cost, we have also tried type-V sources, which
employ loosely converged DMRG wave functions with small bond dimensions of
$M=25,50,$ and $100$. Finally, type-VI sources employ large thresholds
$\epsilon=0.01,$ and $0.003$ to obtain loosely converged HCI wave functions in
the full orbital spaces. This combination has the advantage that it can be
considered a black-box method wherein a single parameter $\epsilon$ controls
the tradeoff between accuracy and cost.

Table 1: Six types of approximate external sources. Detailed explanation to be
found in the main text.

Type | Method | Active space | Cut-off parameters
---|---|---|---
Wave function | Cluster amplitude
I | DMRG | Minimala | $M=2000$ |  $100$% with $s=0.1$
HCI | $\epsilon=10^{-7}$ |  $100$%
II | DMRG | Minimala | $M=2000$ |  $80$% with $s=0.1$
HCI | $\epsilon=10^{-7}$ |  $80$%
III | DMRG | Largerb | $M=2000$ |  $100$% with $s=0.1$
IV | DMRG | Largerb | $M=2000$ |  $80$% with $s=0.1$
V | DMRG | Largerb | $M=25,50,100$ |  $80$% for H2O and N2,
$100$% for F2 with $s=0.01$
VI | HCI | Fullc | $\epsilon=0.01,0.003$ |  $100$%

$\displaystyle{}^{a}\mathrm{(4e,4o)}\ $ $\displaystyle\mathrm{for}\
\mathrm{H}{}_{2}{\mathrm{O}},\ \ $ $\displaystyle\mathrm{(6e,6o)}\ $
$\displaystyle\mathrm{for}\ \mathrm{N}{}_{2}\ \ $
$\displaystyle{}^{b}\mathrm{(8e,18o)}\ $ $\displaystyle\mathrm{for}\
\mathrm{H}{}_{2}{\mathrm{O}},\ \ $ $\displaystyle\mathrm{(10e,16o)}\ $
$\displaystyle\mathrm{for}\ \mathrm{N}{}_{2},\ \ $
$\displaystyle\mathrm{(14e,16o)}\ $ $\displaystyle\mathrm{for}\ \
\mathrm{F}{}_{2}$ $\displaystyle{}^{c}\mathrm{(10e,58o)}\ $
$\displaystyle\mathrm{for}\ \mathrm{H}{}_{2}{\mathrm{O}},\ \ $
$\displaystyle\mathrm{(14e,60o)}\ $ $\displaystyle\mathrm{for}\
\mathrm{N}{}_{2},\ \ $ $\displaystyle\mathrm{(14e,58o)}\ $
$\displaystyle\mathrm{for}\ \ \mathrm{F}{}_{2}$

###  2.3 Perturbative triples corrections

Expressions for the standard, renormalized, and completely renormalized
perturbative triples corrections can be written down in analogy with their
single reference definitions 109. We first define the completely renormalized
(CR)-ecCCSD(T) correction. We use the state $|\Psi^{\rm ecCCSD(T)}\rangle$
defined as

$\displaystyle|\Psi^{\rm ecCCSD(T)}\rangle$
$\displaystyle=(1+T_{1}+T_{2}+Q^{\rm ec}_{3}T_{3}+Q^{\rm
c}_{3}T^{[2]}_{3}+Q^{\rm c}_{3}Z_{3})|\Phi\rangle,$ (7) $\displaystyle
T^{[2]}_{3}|\Phi\rangle$
$\displaystyle=R^{(3)}_{0}(V_{N}T_{2})_{C}|\Phi\rangle,$ (8) $\displaystyle
Z_{3}|\Phi\rangle$ $\displaystyle=R^{(3)}_{0}V_{N}T_{1}|\Phi\rangle.$ (9)

where $V_{N}$ is the two-body part of the Hamiltonian in normal-ordered form
and $R^{(3)}_{0}$ denotes the three-body component of the reduced resolvent
operator in many-body perturbation theory, given by differences of orbital
energies in the denominator.75 The resulting formula for the CR-ecCCSD(T)
energy correction is

$\displaystyle\delta_{0}^{\rm CR-ecCCSD(T)}$
$\displaystyle=\frac{\langle\Psi^{\rm ecCCSD(T)}|Q^{\rm
c}_{3}(H_{N}e^{T})_{C}|\Phi\rangle}{\langle\Psi^{\rm
ecCCSD(T)}|e^{T}|\Phi\rangle}.$ (10)

The energy corrections for renormalized (R)-ecCCSD(T) and ecCCSD(T) can be
obtained by taking the lowest-order estimates of the correction and by
assuming the denominator to be one, i.e.,

$\displaystyle\delta_{0}^{\rm R-ecCCSD(T)}$
$\displaystyle=\frac{\langle\Psi^{\rm ecCCSD(T)}|Q^{\rm
c}_{3}(V_{N}T_{2})_{C}|\Phi\rangle}{\langle\Psi^{\rm
ecCCSD(T)}|e^{T}|\Phi\rangle},$ (11) $\displaystyle\delta_{0}^{\rm ecCCSD(T)}$
$\displaystyle=\langle\Psi^{\rm ecCCSD(T)}|Q^{\rm
c}_{3}(V_{N}T_{2})_{C}|\Phi\rangle.$ (12)

Unlike perturbative triples without external correction (as in CCSD(T)), the
approximate perturbative expression $T_{3}^{[2]}$ is only evaluated for
determinants omitted in the external source. Thus it is not expected to
diverge as long as the external source includes all degeneracies. Nonetheless,
it can still overestimate the triples correlation. The role of the denominator
in the "renormalized" triples approximations is to rescale this correction,
which can be expected to reduce the overestimation.

##  3 Computational Details

All CC calculations were performed using cc-pVTZ basis sets. 110 The ecCC
calculations were performed with the six different types of external sources
summarized in Table 1 and discussed in Section 2.2. The CASSCF-like external
sources (type I and II) used natural orbitals in the minimal active spaces and
CASSCF orbitals for the core and external spaces, while the CASCI-like
external sources (type III-VI) used HF orbitals.

All CC calculations were carried out using a local version of PySCF 111
interfaced with StackBLOCK 112, 113, 114, 115 for DMRG and Arrow 44, 116, 117
for HCI. We used Dice 116, 117, 118 to get accurate SHCI PESs.

##  4 Numerical results

###  4.1 PESs with type-I external sources

The dissociation PESs of H2O and N2, shown in Fig. 1, were obtained by the
ecCC methods using the type-I external sources in Table 1 (i.e. near exact
wave functions in the minimal active spaces and all amplitudes of $T_{3}$ and
$T_{4}$). As expected, the ecCC curves with tightly converged HCI and DMRG
external sources are almost identical. Thus, we show the PESs of ecCC using
only one of these external sources (colored solid lines) and the PESs of CC
(colored dotted lines) in Fig. 1. These are compared against accurate PESs
represented as black lines obtained by SHCI in the full space. The accurate
energies from SHCI are given in the Supporting Information. The mean absolute
errors (MAE) and the non-parallelity errors (NPE) of the PESs are listed in
Table 2.

The CCSD curves (blue dotted lines) have an unphysical dip due to an
inadequate treatment of static correlation. This problem is completely
eliminated within ecCCSD (blue solid line). However, there are significant
errors with respect to the black curve at large distances, giving a (MAE,NPE)
of ($24.0,65.9$) m$E_{\rm H}$ for H2O and ($64.1,149.9$) m$E_{\rm H}$ for N2.

We see that the (T) correction captures much of the missing dynamic
correlation, such that the green ecCCSD(T) curves reach an accuracy of 1 mEH
around the equilibrium geometry. However, when the bond is stretched, the (T)
correction overestimates the dynamic correlation, leading to another
unphysical dip in the PESs. Although this overestimation can easily be reduced
by increasing the size of the active space for the small systems treated here,
this would be very expensive for large systems. Alternatively, we can use the
renormalized triples formula to damp the (T) correction in Eq. (11).
R-ecCCSD(T) (the solid red curves) completely removes the unphysical dips in
the PESs. These attain (MAE,NPE) = ($3.9,23.1$) and ($12.8,41.4$) m$E_{\rm H}$
for H2O and N2, respectively.

Figure 1: PESs of H2O and N2 in the top and bottom panels, respectively,
obtained with the CC methods and the ecCC methods using the type-I external
sources. The type-I sources correspond to near exact wave functions from the
minimal active space, and use all the amplitudes of $T_{3}$ and $T_{4}$. The
blue, green, and red solid lines are the ecCCSD, ecCCSD(T), and R-ecCCSD(T)
PESs, respectively. These are to be compared with the PESs of CC represented
as dotted lines. The black lines are accurate PESs obtained by SHCI.

###  4.2 PESs with type-II external sources

The difference between the type-I and type-II external sources is that the
former use all the amplitudes of $T_{3}$ and $T_{4}$ while the latter use only
a small number of the largest amplitudes, which contribute about 80% of the
total weight of $T_{3}$ and $T_{4}$ (see Table 1 and Section 2.2). For the
minimal active space of H2O, i.e. (4e,4o), the external source contains only
one non-zero $T_{4}$ amplitude, and all elements of $T_{3}$ are zero to within
numerical noise. Thus, the PESs of ecCC using 80% of the external amplitudes
(type-II) are almost identical to those using 100% of the external amplitudes
(type-I). On the other hand, for the minimal active space of N2, i.e. (6e,6o),
the external source contains several large $T_{3}$ and $T_{4}$ amplitudes. At
the stretched geometry corresponding to a bond length of $10a_{0}$, 80% of the
total $T_{3}$ and $T_{4}$ amplitude weight is recovered by the nine largest
elements of $T_{4}$. Figure 2 shows that the type-II external sources,
although using fewer amplitudes, improve the PESs of ecCCSD and R-ecCCSD(T)
relative to the type-I sources. The PES of R-ecCCSD(T) displays a (MAE,NPE) =
($5.4,13.9$) m$E_{\rm H}$.

One concern with partial use of the amplitudes in the minimal active space is
the possibility of divergence in (T). Assuming a HF occupation of the natural
orbitals, the recomputed orbital energies (diagonal parts of the Fock matrix)
of the holes and those of the particles are such that the system retains a
sizable gap. N2 at the $10.0a_{0}$ bond length, for example, has a gap of
around $0.1$ Hartrees between the three occupied orbitals and the three
unoccupied orbitals. This prevents the divergence of the (T) and the
renormalized (T) corrections at the bond lengths we considered (although (T)
largely overestimates dynamic correlation).

Figure 2: PESs of N2 obtained by the ecCC methods using all the external
amplitudes (type-I) and some of the largest external amplitudes (type-II) from
the minimal active space (solid and dotted lines, respectively). Other
descriptions are the same as in Fig.1

###  4.3 PESs with type-III and type-IV external sources

In addition, we investigated PESs of R-ecCCSD(T) using larger active spaces
with all amplitudes (type-III) and 80% of the amplitudes (type-IV). We used
active spaces of (8e,18o), (10e,16o), and (14e,16o) for H2O, N2, and F2,
respectively, and obtained near exact wave functions in the spaces. The PESs
of ecCC using the type-III and type-IV external sources are shown as colored
solid and dotted lines, respectively, in Figure 3. The MAE and NPE of the PESs
are given in Table 2.

For H2O and N2, the resulting PESs of R-ecCCSD(T) with all amplitudes (red
solid lines in the top and middle panels) achieve (MAE,NPE) = ($2.9,18.0$) and
($12.1,44.3$) m$E_{\rm H}$, respectively. These are minor improvements
compared to the results using the minimal active space external sources.
However, the PESs obtained using only 80% of the external amplitudes (red
dotted lines) are much closer to the black, accurate PES, and reach a
(MAE,NPE) = ($1.1,8.1$) and ($4.6,8.4$) m$E_{\rm H}$. Similarly to in the
minimal active space, 80% of the amplitude weight in the larger active space
corresponds to only one element of $T_{4}$ for H2O and the nine largest
elements of $T_{4}$ for N2. When we can find such large elements, the partial
usage of the amplitudes clearly has advantages in both accuracy and
efficiency. It furthermore points to the importance of accounting for the
error from missing relaxation in the external source amplitudes.

Unlike in H2O and N2, there are no particularly large elements in the $T_{3}$
and $T_{4}$ amplitudes of F2. To truncate the amplitudes to 80% of their
weight, we used the approximately 400 largest elements of $T_{3}$ and $T_{4}$
summing to 80% of the total $T_{3}$ and $T_{4}$ weights at the bond length of
$5.0$ Å. The red solid and dotted lines in the bottom panel of Figure 3 show
the PESs of R-ecCCSD(T) using all and 80% of the amplitudes, respectively.
These two PESs are very close to the accurate black curve and reach a
(MAE,NPE) = ($0.9,1.4$) and ($1.7,1.6$) m$E_{\rm H}$ for 100% and 80% of the
amplitudes, respectively. Although both are accurate, the partial use of the
external source in this case (where there is not a small number of large
elements) leads to slightly worse accuracy.

Figure 3: PESs of H2O, N2, and F2 in the top, middle, and bottom panels,
respectively, for the ecCC methods using the larger than minimal active spaces
and the near exact external sources. The solid and dotted lines correspond to
PESs obtained using all external amplitudes of $T_{3}$ and $T_{4}$ (type-III)
and only the largest amplitudes of $T_{4}$ (type-IV), respectively. The
descriptions are otherwise the same as those in Fig.2. Table 2: MAE and NPE
(m$E_{\rm H}$) of H2O, N2, and F2 PESs in a range of geometries
$R\in[0.68,7.80]$ Å, $R\in[1.5,10.0]$ $a_{0}$, and $R\in[1.14,5.00]$ Å,
respectively. ecCC methods using the near exact external sources with all
elements of $T_{3}$ and $T_{4}$ amplitudes ("100%") or only the largest
elements of the $T_{4}$ amplitudes ("80%") from the external source.

Method | $T_{3}$, $T_{4}$ | H2O | N2 | F2
---|---|---|---|---
Active | MAE | NPE | Active | MAE | NPE | Active | MAE | NPE
CCSD |  |  | 15.9 | 30.4 |  | 32.4 | 102.3 |  | 34.8 | 57.7
ecCCSD | 100% | 4e,4o | 24.0 | 65.9 | 6e,6o | 64.1 | 149.9 | 2e,2o | a | a
80% | 24.0 | 65.9 | 57.5 | 114.5
100% | 8e,18o | 17.6 | 49.9 | 10e,16o | 48.4 | 113.0 | 14e,16o | 25.2 | 41.4
80% | 15.6 | 39.4 | 41.4 | 74.0 | 27.0 | 40.9
R-CCSD(T) |  |  | 6.8 | 43.0 |  | 32.1 | 162.1 |  | 6.2 | 10.6
R-ecCCSD(T) | 100% | 4e,4o | 3.9 | 23.1 | 6e,6o | 12.8 | 41.4 | 2e,2o | a | a
80% | 3.9 | 23.1 | 5.4 | 13.9
100% | 8e,18o | 2.9 | 18.0 | 10e,16o | 12.1 | 44.3 | 14e,16o | 0.9 | 1.4
80% | 1.1 | 8.1 | 4.6 | 8.4 | 1.7 | 1.6

a no $T_{3}$ and $T_{4}$ in the minimal active space of F2

###  4.4 PESs with type-V external sources

In the previous section, we showed that the use of larger active spaces
significantly improves the ecCC PESs. However, obtaining tightly converged
DMRG wave functions in large active spaces requires more CPU time than the
subsequent CC calculation. For example, optimizing an external wave function
with 16 orbitals requires around a few minutes of CPU time for one DMRG sweep
with $M=2000$. Although a few minutes is not prohibitively large in many
applications, it is large compared to the subsequent ecCC calculation which
only takes tens of seconds, at least for the small molecules considered here.
In addition, the fact that in some cases only a partial use of the amplitudes
led to better results in the last section suggests that it is not a good use
of computational time to tightly converge the external source.

We thus now consider loosely converged DMRG sources (type-V) in the larger
active spaces. We re-computed PESs of R-ecCCSD(T) shown in Figure 3 using a
DMRG source with bond dimensions 25, 50, and 100. We used the truncation to
80% amplitude weight for H2O and N2, while we used all the external amplitudes
for F2. Table 3 shows the corresponding MAE and NPE in the same range of
geometries in Table 2.

For all cases, when we reduced the bond dimension to $M=100$, the MAE
increased by $0.1\sim 0.5$ m$E_{H}$, and the NPE increased by $0.3\sim 1.0$
m$E_{H}$, compared to using $M=2000$. However for $M=100$, one DMRG sweep with
16 orbitals took only a few seconds of CPU time at the 10$a_{0}$ bond length
of N2, giving a better computational balance between the DMRG calculation and
the subsequent ecCC calculations. When we further reduced the bond dimension
to $M=25$, the MAE increased by $0.9$ m$E_{H}$ and the NPE increased by $0.6$
and $2.7$ m$E_{H}$ for H2O and N2. In the case of F2, which does not have a
small number of large $T_{3}$ and $T_{4}$ elements, the MAE and NPE increased
more, by $3.4$ and $5.4$ m$E_{H}$.

###  4.5 PESs with type-VI external sources

Type VI sources use the full orbital space and employ the single HCI
parameter, $\epsilon$, to select the large $T_{3}$ and $T_{4}$ amplitudes.
They do not require a choice for the CAS space, and we did not further
threshold the $T_{3}$ and $T_{4}$ amplitudes. However, for systems where the
molecule has a low-spin ground state, but the dissociated fragments have high-
spin ground states, the HCI wave functions at stretched geometries are
strongly spin-contaminated when $\epsilon$ is large, even when time-reversal
symmetry is employed to reduce the spin contamination. (This type of spin-
contamination is avoided in the small bond dimension DMRG calculations via the
full use of spin symmetry). Although this has little effect on the HCI energy,
it makes the amplitudes unsuitable for ecCC calculations.

In the second part of Table 3 the MAE and NPE from type-VI sources using
$\epsilon=0.01$ and $0.003$ are shown. These are evaluated over a smaller
range of geometries than those used in the first part of Table 3 for the
reason mentioned above. Since N2 is a singlet that dissociates into atoms that
are quartets, it has particularly large errors. The errors improve upon going
from $\epsilon=0.01$ to $\epsilon=0.003$, in particular the PEC has an
unphysical dip at $\epsilon=0.01$ which disappears at $\epsilon=0.003$, but
the errors are still substantial. Similarly to how using all the amplitudes
could lead to larger errors than only using some of the amplitudes in the
previous sections, the errors incurred from the type VI sources here emphasize
that the quality of the external amplitudes cannot be judged solely from the
energy obtained by the external method.

Table 3: MAE and NPE (m$E_{\rm H}$) of R-ecCCSD(T) PESs for H2O, N2, and F2
using the approximate DMRG and HCI sources with a small bond dimension ($M$)
or large threshold ($\epsilon$). Detailed explanation to be found in the main
text.

$M$ | H2O | N2 | F2
---|---|---|---
Active | MAE a | NPE a | Active | MAE a | NPE a | Active | MAE a | NPE a
25 | (8e,18o) | 3.8 | 8.7 | (10e,16o) | 5.5 | 11.1 | (14e,16o) | 4.3 | 6.8
50 | 3.4 | 8.2 | 5.2 | 9.6 | 3.5 | 3.4
100 | 3.0 | 8.4 | 4.9 | 9.4 | 1.4 | 2.1
2000 | 2.9 | 8.1 | 4.6 | 8.4 | 0.9 | 1.4
$\epsilon$ | H2O | N2 | F2
Active | MAE b | NPE b | Active | MAE b | NPE b | Active | MAE b | NPE b
$0.01$ | (10e,58o) | 1.6 | 7.1 | (14e,60o) | 8.1 | 61.5 | (14e,58o) | 3.6 | 4.7
$0.003$ | 1.3 | 4.4 | 3.6 | 15.8 | 3.0 | 5.8

a $R\in[0.68,7.80]$ Å, $[1.5,10.0]$ $a_{0}$, and $[1.14,5.00]$ Å for H2O, N2,
and F2 PESs, respectively.
b $R\in[0.68,3.41]$ Å, $[1.5,5.4]$ $a_{0}$, and $[1.14,2.80]$ Å for H2O, N2,
and F2 PESs, respectively.

###  4.6 Error analysis

Figure 4: Absolute errors of R-CCSD(T) and ecCC methods on a log scale plotted
against the $D_{2}$ diagnostic. Detailed explanation in the main text.

In this section, we present the errors of R-ecCCSD(T) for the systems in this
work and analyze them using the well known CC error diagnostic $D_{2}$,
defined by the matrix 2-norm of the $T_{2}$ amplitudes. 119 The magnitude of
$D_{2}$ can be used to distinguish between the SR and MR character of the
different geometries on the PESs. Organic molecules are sometimes considered
to have MR character when $D_{2}$ is larger than $0.18$. 119 Figure 4 shows
the absolute errors on a log scale versus the $D_{2}$ diagnostic. Each symbol
represents a geometry on the PESs of one of the molecules, H2O, N2, and F2.
Square symbols denote R-CCSD(T) and circle symbols denote the various
externally corrected theories, using the type-V and the type-VI external
sources. For $D_{2}$ ranging from $0$ to $1.7$, the absolute errors of
R-ecCCSD(T) (red circles) are less than $0.015$ EH and most of the errors are
smaller than those of R-CCSD(T) (black squares) and ecCCSD (blue circles). The
absolute errors of R-ecCCSD(T) are less than those of ecCCSD(T) (green
circles) in the MR region where $D_{2}>0.4$, while they are mostly greater
than those of ecCCSD(T) in the range $0.0<D_{2}<0.4$. Overall, it is clear to
see that R-ecCCSD(T) offers the most balanced treatment of errors across a
wide range of SR and MR character. However, for very weakly correlated
systems, the original (T) correction is slightly more accurate than the
renormalized (T) correction.

###  4.7 Limitations

Although the work above shows that it is possible to obtain quantitative
accuracy across the full potential energy surface, at a reasonable
computational cost, using the combination of external sources and the
renormalized (T) correction, SR ecCC approaches have a fundamental limitation
when the CI coefficient of the reference configuration (the HF configuration
in this work) in the external source is exactly or numerically zero. This
leads to overflow due to the assumption of intermediate normalization. The
smallest reference coefficient value we encountered in this work was $0.2$ at
the stretched $10a_{0}$ geometry of N2. Although this coefficient is not
numerically zero, it is difficult to converge the ecCCSD energy. We extracted
initial amplitudes of $T_{1}$ and $T_{2}$ from the external source and then
iteratively converged the ecCCSD energy using a damping parameter of $0.05$ to
update the amplitudes. (We did not use the direct inversion in the iterative
subspace (DIIS) algorithm). At this geometry, the energy could be converged
monotonically up to a threshold of $10^{-5}$ Hartrees, although the norm of
the amplitudes could not be converged, and we simply used amplitudes from the
last iteration in the set of monotonically decreasing energies.

##  5 Conclusion

In this work, we explored externally corrected coupled cluster with a
renormalized triples correction (R-ecCCSD(T)) using DMRG and HCI and external
sources. The critical question is how to best balance the accuracy and cost of
computing the external source with the cost of the overall method. To this
end, we considered multiple types of external sources: "exact" external
sources, where the DMRG and HCI wavefunctions were tightly converged, and
"approximate" external sources where they were loosely converged. We also
considered both "full" usage of the $T_{3}$ and $T_{4}$ amplitudes, and
"partial" usage, wherein we retained only the largest elements.

For all systems considered here, we found that R-ecCCSD(T) can significantly
improve on the description of the potential energy surface given either by the
external source alone or CC alone. For example, the unphysical dips in the PES
of H2O and N2 in the bond-stretched region can be completely eliminated. The
use of approximate external sources, possibly with truncation to only the
large $T_{3}$ and $T_{4}$ amplitudes, appears to be a practical way to balance
the cost of the external calculation and the coupled cluster calculation in
small molecules. Using the $D_{2}$ diagnostic to characterize the different
points on the potential energy surfaces, we find that R-ecCCSD(T) gives
absolute errors of less than 15 m$E_{H}$ in the range of $D_{2}$ from $0.0$ to
$1.7$. In fact, the errors of R-ecCCSD(T) are less than those of ecCCSD and
R-CCSD(T) in almost all cases, except when $D_{2}$ is very small, where the
renormalized (T) correction appears to be slightly less accurate than the
simple (T) correction.

There are several interesting questions remaining which lie beyond what we
have considered in this work. For example, while R-ecCCSD(T) appears quite
stable up to large values of the $D_{2}$ diagnostic, what is the largest
amount of multi-reference character which can be handled? Here, the difficulty
in solving the CC equations, and the divergence of the amplitudes reflecting
the problems of intermediate normalization, cannot be ignored. In addition, in
the realm of quasidegenerate problems, we can ask whether other non-iterative
corrections such as the "completely renormalized" triples and quadruples
corrections, 109 corresponding to CR-ecCCSD(T) and CR-ecCCSD(T,Q), would
further improve on the present R-ecCCSD(T) method. Finally, the current
approximation, with its modest computational requirements on the external
source, is applicable to the same scale of systems that can be handled by
single reference coupled cluster methods. Thus understanding the performance
of this method in larger correlated systems is of interest.

Note: As this work was being prepared for submission, we were made aware of a
related recent submission to the arxiv 120, that also discusses selected
configuration interaction as an external source and perturbative triples
corrections in the context of externally corrected coupled cluster methods.

Work by SL, HZ and GKC was supported by the US National Science Foundation via
Award CHE-1655333. GKC is a Simons Investigator. SS was supported by the US
National Science Foundation grant CHE-1800584 and the Sloan research
fellowship. CJU was supported in part by the AFOSR under grant
FA9550-18-1-0095.

## References

  * Abe 2013 Abe, M. Diradicals. _Chem. Rev._ 2013, _113_ , 7011–7088
  * Nakano 2017 Nakano, M. Electronic Structure of Open-Shell Singlet Molecules: Diradical Character Viewpoint. _Top. Curr. Chem._ 2017, _375_
  * Kurashige et al. 2013 Kurashige, Y.; Chan, G. K.-L.; Yanai, T. Entangled quantum electronic wavefunctions of the Mn 4 CaO 5 cluster in photosystem II. _Nature chemistry_ 2013, _5_ , 660–666
  * Sharma et al. 2014 Sharma, S.; Sivalingam, K.; Neese, F.; Chan, G. K.-L. Low-energy spectrum of iron–sulfur clusters directly from many-particle quantum mechanics. _Nat. Chem._ 2014, _6_ , 927–933
  * Li et al. 2019 Li, Z.; Guo, S.; Sun, Q.; Chan, G. K.-L. Electronic landscape of the P-cluster of nitrogenase as revealed through many-electron quantum wavefunction simulations. _Nat. Chem._ 2019, _11_ , 1026–1033
  * Das and Wahl 1972 Das, G.; Wahl, A. C. New Techniques for the Computation of Multiconfiguration Self-Consistent Field (MCSCF) Wavefunctions. _The Journal of Chemical Physics_ 1972, _56_ , 1769–1775
  * Werner and Meyer 1980 Werner, H.-J.; Meyer, W. A quadratically convergent multiconfiguration–self-consistent field method with simultaneous optimization of orbitals and CI coefficients. _The Journal of Chemical Physics_ 1980, _73_ , 2342–2356
  * Werner and Knowles 1985 Werner, H.-J.; Knowles, P. J. A second order multiconfiguration SCF procedure with optimum convergence. _The Journal of chemical physics_ 1985, _82_ , 5053–5063
  * White and Noack 1992 White, S. R.; Noack, R. M. Real-space quantum renormalization groups. _Physical review letters_ 1992, _68_ , 3487
  * White 1992 White, S. R. Density matrix formulation for quantum renormalization groups. _Physical review letters_ 1992, _69_ , 2863
  * White 1993 White, S. R. Density-matrix algorithms for quantum renormalization groups. _Physical Review B_ 1993, _48_ , 10345
  * White and Martin 1999 White, S. R.; Martin, R. L. Ab initio quantum chemistry using the density matrix renormalization group. _The Journal of chemical physics_ 1999, _110_ , 4127–4130
  * Chan and Head-Gordon 2002 Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. _The Journal of chemical physics_ 2002, _116_ , 4462–4476
  * Legeza et al. 2003 Legeza, Ö.; Röder, J.; Hess, B. Controlling the accuracy of the density-matrix renormalization-group method: The dynamical block state selection approach. _Physical Review B_ 2003, _67_ , 125114
  * Chan 2004 Chan, G. K.-L. An algorithm for large scale density matrix renormalization group calculations. _The Journal of chemical physics_ 2004, _120_ , 3172–3178
  * Legeza et al. 2008 Legeza, Ö.; Noack, R.; Sólyom, J.; Tincani, L. _Computational Many-Particle Physics_ ; Springer, 2008; Vol. 739
  * Chan and Zgid 2009 Chan, G.; Zgid, D. The Density Matrix Renormalization Group in Quantum Chemistry, Annu. Rep. _Comput. Chem_ 2009, _5_ , 149–162
  * Marti and Reiher 2010 Marti, K. H.; Reiher, M. The density matrix renormalization group algorithm in quantum chemistry. _Zeitschrift für Physikalische Chemie_ 2010, _224_ , 583–599
  * Chan and Sharma 2011 Chan, G. K.-L.; Sharma, S. The density matrix renormalization group in quantum chemistry. _Annual review of physical chemistry_ 2011, _62_ , 465–481
  * Sharma and Chan 2012 Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithms for quantum chemistry. _The Journal of chemical physics_ 2012, _136_ , 124121
  * Wouters and Van Neck 2014 Wouters, S.; Van Neck, D. The density matrix renormalization group for ab initio quantum chemistry. _The European Physical Journal D_ 2014, _68_ , 1–20
  * Szalay et al. 2015 Szalay, S.; Pfeffer, M.; Murg, V.; Barcza, G.; Verstraete, F.; Schneider, R.; Legeza, Ö. Tensor product methods and entanglement optimization for ab initio quantum chemistry. _International Journal of Quantum Chemistry_ 2015, _115_ , 1342
  * Yanai et al. 2015 Yanai, T.; Kurashige, Y.; Mizukami, W.; Chalupskỳ, J.; Lan, T. N.; Saitow, M. Density matrix renormalization group for ab initio Calculations and associated dynamic correlation methods: A review of theory and applications. _International Journal of Quantum Chemistry_ 2015, _115_ , 283–299
  * Ivanic and Ruedenberg 2001 Ivanic, J.; Ruedenberg, K. Identification of deadwood in configuration spaces through general direct configuration interaction. _Theoretical Chemistry Accounts_ 2001, _106_ , 339–351
  * Huron et al. 1973 Huron, B.; Malrieu, J.; Rancurel, P. Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions. _The Journal of Chemical Physics_ 1973, _58_ , 5745–5759
  * Garniron et al. 2018 Garniron, Y.; Scemama, A.; Giner, E.; Caffarel, M.; Loos, P.-F. Selected configuration interaction dressed by perturbation. _J. Chem. Phys._ 2018, _149_ , 064103
  * Loos et al. 2018 Loos, P.-F.; Scemama, A.; Blondel, A.; Garniron, Y.; Caffarel, M.; Jacquemin, D. A Mountaineering Strategy to Excited States: Highly Accurate Reference Energies and Benchmarks. _J. Chem. Theory Comput._ 2018, _14_ , 4360–4379
  * Garniron et al. 2017 Garniron, Y.; Scemama, A.; Loos, P.-F.; Caffarel, M. Hybrid stochastic-deterministic calculation of the second-order perturbative contribution of multireference perturbation theory. _J. Chem. Phys._ 2017, _147_
  * Buenker and Peyerimhoff 1974 Buenker, R. J.; Peyerimhoff, S. D. Individualized configuration selection in CI calculations with subsequent energy extrapolation. _Theoretica chimica acta_ 1974, _35_ , 33–58
  * Evangelisti et al. 1983 Evangelisti, S.; Daudey, J.-P.; Malrieu, J.-P. Convergence of an improved CIPSI algorithm. _Chemical Physics_ 1983, _75_ , 91–102
  * Harrison 1991 Harrison, R. J. Approximating full configuration interaction with selected configuration interaction and perturbation theory. _The Journal of chemical physics_ 1991, _94_ , 5021–5031
  * Steiner et al. 1994 Steiner, M.; Wenzel, W.; Wilson, K.; Wilkins, J. The efficient treatment of higher excitations in CI calculations: A comparison of exact and approximate results. _Chemical Physics Letters_ 1994, _231_ , 263 – 268
  * Neese 2003 Neese, F. A spectroscopy oriented configuration interaction procedure. _The Journal of Chemical Physics_ 2003, _119_ , 9428–9443
  * Abrams and Sherrill 2005 Abrams, M. L.; Sherrill, C. D. Important configurations in configuration interaction and coupled-cluster wave functions. _Chemical Physics Letters_ 2005, _412_ , 121 – 124
  * Bytautas and Ruedenberg 2009 Bytautas, L.; Ruedenberg, K. A priori identification of configurational deadwood. _Chemical Physics_ 2009, _356_ , 64 – 75, Moving Frontiers in Quantum Chemistry:
  * Evangelista 2014 Evangelista, F. A. A driven similarity renormalization group approach to quantum many-body problems. _The Journal of chemical physics_ 2014, _141_ , 054109
  * Knowles 2015 Knowles, P. J. Compressive sampling in configuration interaction wavefunctions. _Molecular Physics_ 2015, _113_ , 1655–1660
  * Schriber and Evangelista 2016 Schriber, J. B.; Evangelista, F. A. Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy. _The Journal of Chemical Physics_ 2016, _144_ , 161106
  * Tubman et al. 2016 Tubman, N. M.; Lee, J.; Takeshita, T. Y.; Head-Gordon, M.; Whaley, K. B. A deterministic alternative to the full configuration interaction quantum Monte Carlo method. _The Journal of chemical physics_ 2016, _145_ , 044112
  * Liu and Hoffmann 2016 Liu, W.; Hoffmann, M. R. iCI: Iterative CI toward full CI. _Journal of chemical theory and computation_ 2016, _12_ , 1169–1178
  * Caffarel et al. 2016 Caffarel, M.; Applencourt, T.; Giner, E.; Scemama, A. _Recent Progress in Quantum Monte Carlo_ ; ACS Publications, 2016; pp 15–46
  * Holmes et al. 2016 Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-Bath Configuration Interaction: An Efficient Selected Configuration Interaction Algorithm Inspired by Heat-Bath Sampling. _J. Chem. Theory Comput._ 2016, _12_ , 3674–3680
  * Sharma et al. 2017 Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochastic Heat-Bath Configuration Interaction Method: Selected Configuration Interaction with Semistochastic Perturbation Theory. _Journal of Chemical Theory and Computation_ 2017, _13_ , 1595–1604
  * Li et al. 2018 Li, J.; Otten, M.; Holmes, A. A.; Sharma, S.; Umrigar, C. J. Fast Semistochastic Heat-Bath Configuration Interaction. _J. Chem. Phys._ 2018, _148_ , 214110
  * Booth et al. 2009 Booth, G. H.; Thom, A. J.; Alavi, A. Fermion Monte Carlo without fixed nodes: A game of life, death, and annihilation in Slater determinant space. _The Journal of chemical physics_ 2009, _131_ , 054106
  * Cleland et al. 2010 Cleland, D.; Booth, G. H.; Alavi, A. Communications: Survival of the fittest: Accelerating convergence in full configuration-interaction quantum Monte Carlo. 2010
  * Ghanem et al. 2019 Ghanem, K.; Lozovoi, A. Y.; Alavi, A. Unbiasing the initiator approximation in full configuration interaction quantum Monte Carlo. _The Journal of chemical physics_ 2019, _151_ , 224108
  * Sabzevari and Sharma 2018 Sabzevari, I.; Sharma, S. Improved speed and scaling in orbital space variational Monte Carlo. _Journal of chemical theory and computation_ 2018, _14_ , 6276–6286
  * Neuscamman 2016 Neuscamman, E. Subtractive manufacturing with geminal powers: making good use of a bad wave function. _Molecular Physics_ 2016, _114_ , 577–583
  * Andersson et al. 1990 Andersson, K.; Malmqvist, P. A.; Roos, B. O.; Sadlej, A. J.; Wolinski, K. Second-order perturbation theory with a CASSCF reference function. _Journal of Physical Chemistry_ 1990, _94_ , 5483–5488
  * Angeli et al. 2001 Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. Introduction of n-electron valence states for multireference perturbation theory. _The Journal of Chemical Physics_ 2001, _114_ , 10252–10264
  * Angeli et al. 2001 Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. N-electron valence state perturbation theory: a fast implementation of the strongly contracted variant. _Chemical physics letters_ 2001, _350_ , 297–305
  * Angeli et al. 2002 Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. n-electron valence state perturbation theory: A spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. _The Journal of chemical physics_ 2002, _117_ , 9138–9153
  * Neuscamman et al. 2010 Neuscamman, E.; Yanai, T.; Chan, G. K.-L. A review of canonical transformation theory. _International Reviews in Physical Chemistry_ 2010, _29_ , 231–271
  * Yanai et al. 2010 Yanai, T.; Kurashige, Y.; Neuscamman, E.; Chan, G. K.-L. Multireference quantum chemistry through a joint density matrix renormalization group and canonical transformation theory. _The Journal of chemical physics_ 2010, _132_ , 024105
  * Kurashige and Yanai 2011 Kurashige, Y.; Yanai, T. Second-order perturbation theory with a density matrix renormalization group self-consistent field reference function: Theory and application to the study of chromium dimer. _The Journal of chemical physics_ 2011, _135_ , 094104
  * Sharma and Chan 2014 Sharma, S.; Chan, G. K.-L. Communication: A flexible multi-reference perturbation theory by minimizing the Hylleraas functional with matrix product states. 2014
  * Sharma and Alavi 2015 Sharma, S.; Alavi, A. Multireference linearized coupled cluster theory for strongly correlated systems using matrix product states. _The Journal of chemical physics_ 2015, _143_ , 102815
  * Guo et al. 2016 Guo, S.; Watson, M. A.; Hu, W.; Sun, Q.; Chan, G. K.-L. N-electron valence state perturbation theory based on a density matrix renormalization group reference function, with applications to the chromium dimer and a trimer model of poly (p-phenylenevinylene). _Journal of chemical theory and computation_ 2016, _12_ , 1583–1591
  * Sokolov et al. 2017 Sokolov, A. Y.; Guo, S.; Ronca, E.; Chan, G. K.-L. Time-dependent N-electron valence perturbation theory with matrix product state reference wavefunctions for large active spaces and basis sets: Applications to the chromium dimer and all-trans polyenes. _The Journal of chemical physics_ 2017, _146_ , 244102
  * Sharma et al. 2017 Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochastic heat-bath configuration interaction method: Selected configuration interaction with semistochastic perturbation theory. _Journal of chemical theory and computation_ 2017, _13_ , 1595–1604
  * Guo et al. 2018 Guo, S.; Li, Z.; Chan, G. K.-L. Communication: An efficient stochastic algorithm for the perturbative density matrix renormalization group in large active spaces. _The Journal of chemical physics_ 2018, _148_ , 221104
  * Lischka et al. 1981 Lischka, H.; Shepard, R.; Brown, F. B.; Shavitt, I. New implementation of the graphical unitary group approach for multireference direct configuration interaction calculations. _International Journal of Quantum Chemistry_ 1981, _20_ , 91–100
  * Szalay et al. 2012 Szalay, P. G.; Muller, T.; Gidofalvi, G.; Lischka, H.; Shepard, R. Multiconfiguration self-consistent field and multireference configuration interaction methods and applications. _Chemical reviews_ 2012, _112_ , 108–181
  * Lischka et al. 2001 Lischka, H.; Shepard, R.; Pitzer, R. M.; Shavitt, I.; Dallos, M.; Müller, T.; Szalay, P. G.; Seth, M.; Kedziora, G. S.; Yabushita, S., et al. High-level multireference methods in the quantum-chemistry program system COLUMBUS: Analytic MR-CISD and MR-AQCC gradients and MR-AQCC-LRT for excited states, GUGA spin–orbit CI and parallel CI density. _Physical Chemistry Chemical Physics_ 2001, _3_ , 664–673
  * Saitow et al. 2015 Saitow, M.; Kurashige, Y.; Yanai, T. Fully internally contracted multireference configuration interaction theory using density matrix renormalization group: A reduced-scaling implementation derived by computer-aided tensor factorization. _Journal of chemical theory and computation_ 2015, _11_ , 5120–5131
  * Paldus and Li 1999 Paldus, J.; Li, X. Adv. Chem. Phys. 1999, 1–175
  * Kong et al. 2009 Kong, L.; Shamasundar, K.; Demel, O.; Nooijen, M. State specific equation of motion coupled cluster method in general active space. _The Journal of chemical physics_ 2009, _130_ , 114101
  * Jeziorski 2010 Jeziorski, B. Multireference coupled-cluster Ansatz. _Molecular Physics_ 2010, _108_ , 3043–3054
  * Čížek 1966 Čížek, J. On the correlation problem in atomic and molecular systems. Calculation of wavefunction components in Ursell-type expansion using quantum-field theoretical methods. _The Journal of Chemical Physics_ 1966, _45_ , 4256–4266
  * Čížek 1969 Čížek, J. On the use of the cluster expansion and the technique of diagrams in calculations of correlation effects in atoms and molecules. _Advances in chemical physics_ 1969, 35–89
  * Čížek and Paldus 1971 Čížek, J.; Paldus, J. Correlation problems in atomic and molecular systems III. Rederivation of the coupled-pair many-electron theory using the traditional quantum chemical methodst. _International Journal of Quantum Chemistry_ 1971, _5_ , 359–379
  * Paldus et al. 1972 Paldus, J.; Čížek, J.; Shavitt, I. Correlation Problems in Atomic and Molecular Systems. IV. Extended Coupled-Pair Many-Electron Theory and Its Application to the B H 3 Molecule. _Physical Review A_ 1972, _5_ , 50
  * Bartlett 1995 Bartlett, R. J. In _Modern Electronic Structure Theory_ ; Yarkony, D. R., Ed.; World Scientific, 1995; Vol. 1; pp 1047–1131
  * Shavitt and Bartlett 2009 Shavitt, I.; Bartlett, R. J. _Many-body methods in chemistry and physics: MBPT and coupled-cluster theory_ ; Cambridge university press, 2009
  * Kinoshita et al. 2005 Kinoshita, T.; Hino, O.; Bartlett, R. J. Coupled-cluster method tailored by configuration interaction. _The Journal of chemical physics_ 2005, _123_ , 074106
  * Hino et al. 2006 Hino, O.; Kinoshita, T.; Chan, G. K.-L.; Bartlett, R. J. Tailored coupled cluster singles and doubles method applied to calculations on molecular structure and harmonic vibrational frequencies of ozone. _The Journal of chemical physics_ 2006, _124_ , 114311
  * Lyakh et al. 2011 Lyakh, D. I.; Lotrich, V. F.; Bartlett, R. J. The ‘tailored’CCSD (T) description of the automerization of cyclobutadiene. _Chemical Physics Letters_ 2011, _501_ , 166–171
  * Veis et al. 2016 Veis, L.; Antalík, A.; Brabec, J.; Neese, F.; Legeza, O.; Pittner, J. Coupled cluster method with single and double excitations tailored by matrix product state wave functions. _The journal of physical chemistry letters_ 2016, _7_ , 4072–4078
  * Faulstich et al. 2019 Faulstich, F. M.; Máté, M.; Laestadius, A.; Csirik, M. A.; Veis, L.; Antalik, A.; Brabec, J.; Schneider, R.; Pittner, J.; Kvaal, S., et al. Numerical and theoretical aspects of the DMRG-TCC method exemplified by the nitrogen dimer. _Journal of chemical theory and computation_ 2019, _15_ , 2206–2220
  * Vitale et al. 2020 Vitale, E.; Alavi, A.; Kats, D. FCIQMC-tailored distinguishable cluster approach. _Journal of Chemical Theory and Computation_ 2020, _16_ , 5621–5634
  * Paldus and Planelles 1994 Paldus, J.; Planelles, J. Valence bond corrected single reference coupled cluster approach. _Theoretica chimica acta_ 1994, _89_ , 13–31
  * Planelles et al. 1994 Planelles, J.; Paldus, J.; Li, X. Valence bond corrected single reference coupled cluster approach. _Theoretica chimica acta_ 1994, _89_ , 33
  * Planelles et al. 1994 Planelles, J.; Paldus, J.; Li, X. Valence bond corrected single reference coupled cluster approach. _Theoretica chimica acta_ 1994, _89_ , 59
  * Paldus 2017 Paldus, J. Externally and internally corrected coupled cluster approaches: an overview. _Journal of Mathematical Chemistry_ 2017, _55_ , 477–502
  * Paldus and Boyle 1982 Paldus, J.; Boyle, M. Cluster analysis of the full configuration interaction wave functions of cyclic polyene models. _International Journal of Quantum Chemistry_ 1982, _22_ , 1281–1305
  * Paldus et al. 1983 Paldus, J.; Chin, E.; Grey, M. Bond length alternation in cyclic polyenes. II. Unrestricted hartree–fock method. _International journal of quantum chemistry_ 1983, _24_ , 395–409
  * Paldus et al. 1984 Paldus, J.; Čížek, J.; Takahashi, M. Approximate account of the connected quadruply excited clusters in the coupled-pair many-electron theory. _Physical Review A_ 1984, _30_ , 2193
  * Piecuch et al. 1996 Piecuch, P.; Tobol, R.; Paldus, J. Approximate account of connected quadruply excited clusters in single-reference coupled-cluster theory via cluster analysis of the projected unrestricted Hartree-Fock wave function. _Physical Review A_ 1996, _54_ , 1210
  * Peris et al. 1997 Peris, G.; Planelles, J.; Paldus, J. Single-reference CCSD approach employing three-and four-body CAS SCF corrections: A preliminary study of a simple model. _International journal of quantum chemistry_ 1997, _62_ , 137–151
  * Li et al. 1997 Li, X.; Peris, G.; Planelles, J.; Rajadall, F.; Paldus, J. Externally corrected singles and doubles coupled cluster methods for open-shell systems. _The Journal of chemical physics_ 1997, _107_ , 90–98
  * Peris et al. 1998 Peris, G.; Rajadell, F.; Li, X.; Planelles, J.; Paldus, J. Externally corrected singles and doubles coupled cluster methods for open-shell systems. II. Applications to the low lying doublet states of OH, NH2, CH3 and CN radicals. _Molecular Physics_ 1998, _94_ , 235–248
  * Stolarczyk 1994 Stolarczyk, L. Z. Complete active space coupled-cluster method. Extension of single-reference coupled-cluster method using the CASSCF wavefunction. _Chemical physics letters_ 1994, _217_ , 1–6
  * Xu and Li 2015 Xu, E.; Li, S. The externally corrected coupled cluster approach with four-and five-body clusters from the CASSCF wave function. _The Journal of chemical physics_ 2015, _142_ , 094119
  * Deustua et al. 2018 Deustua, J. E.; Magoulas, I.; Shen, J.; Piecuch, P. Communication: Approaching exact quantum chemistry by cluster analysis of full configuration interaction quantum Monte Carlo wave functions. _The Journal of Chemical Physics_ 2018, _149_ , 151101
  * Li and Paldus 1997 Li, X.; Paldus, J. Reduced multireference CCSD method: An effective approach to quasidegenerate states. _The Journal of Chemical Physics_ 1997, _107_ , 6257–6269
  * Li and Paldus 2000 Li, X.; Paldus, J. Truncated version of the reduced multireference coupled-cluster method with perturbation selection of higher than pair clusters. _International Journal of Quantum Chemistry_ 2000, _80_ , 743–756
  * Li and Paldus 2006 Li, X.; Paldus, J. Reduced multireference coupled cluster method with singles and doubles: Perturbative corrections for triples. _The Journal of Chemical Physics_ 2006, _124_ , 174101
  * Li and Paldus 1998 Li, X.; Paldus, J. Dissociation of N2 triple bond: a reduced multireference CCSD study. _Chemical physics letters_ 1998, _286_ , 145–154
  * Li and Paldus 1998 Li, X.; Paldus, J. Singlet-triplet splitting in methylene: An accurate description of dynamic and nondynamic correlation by reduced multireference coupled cluster method. _Collection of Czechoslovak chemical communications_ 1998, _63_ , 1381–1393
  * Li and Paldus 1999 Li, X.; Paldus, J. Simultaneous handling of dynamical and nondynamical correlation via reduced multireference coupled cluster method: Geometry and harmonic force field of ozone. _The Journal of chemical physics_ 1999, _110_ , 2844–2852
  * Li and Paldus 2000 Li, X.; Paldus, J. Reduced multireference coupled cluster method: Ro-vibrational spectra of N 2. _The Journal of Chemical Physics_ 2000, _113_ , 9966–9977
  * Li and Paldus 2000 Li, X.; Paldus, J. Effect of spin contamination on the prediction of barrier heights by coupled-cluster theory: F+ H2→ HF+ H reaction. _International Journal of Quantum Chemistry_ 2000, _77_ , 281–290
  * Li and Paldus 2006 Li, X.; Paldus, J. A truncated version of reduced multireference coupled-cluster method with singles and doubles and noniterative triples: Application to F 2 and Ni (CO) n (n= 1, 2, and 4). _The Journal of chemical physics_ 2006, _125_ , 164107
  * Li and Paldus 2006 Li, X.; Paldus, J. Singlet–triplet separation in BN and C2: simple yet exceptional systems for advanced correlated methods. _Chemical physics letters_ 2006, _431_ , 179–184
  * Li and Paldus 2007 Li, X.; Paldus, J. Reduced multireference coupled-cluster method: Barrier heights for heavy atom transfer, nucleophilic substitution, association, and unimolecular reactions. _The Journal of Physical Chemistry A_ 2007, _111_ , 11189–11197
  * Li and Paldus 2008 Li, X.; Paldus, J. Coupled-cluster approach to spontaneous symmetry breaking in molecules: The linear N3 radical. _International Journal of Quantum Chemistry_ 2008, _108_ , 2117–2127
  * Li and Paldus 2008 Li, X.; Paldus, J. Full potential energy curve for N2 by the reduced multireference coupled-cluster method. _The Journal of Chemical Physics_ 2008, _129_ , 054104
  * Kowalski and Piecuch 2000 Kowalski, K.; Piecuch, P. The method of moments of coupled-cluster equations and the renormalized CCSD[T], CCSD(T), CCSD(TQ), and CCSDT(Q) approaches. _J. Chem. Phys._ 2000, _113_ , 18
  * Dunning 1989 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. _J. Chem. Phys._ 1989, _90_ , 1007–1023
  * Sun et al. 2017 Sun, Q.; Berkelbach, T. C.; Blunt, N. S.; Booth, G. H.; Guo, S.; Li, Z.; Liu, J.; McClain, J. D.; Sayfutyarova, E. R.; Sharma, S.; Wouters, S.; Chan, G. K. PySCF: the Python-based simulations of chemistry framework. 2017; https://onlinelibrary.wiley.com/doi/abs/10.1002/wcms.1340
  * Chan and Head-Gordon 2002 Chan, G. K.-L.; Head-Gordon, M. Highly correlated calculations with a polynomial cost algorithm: A study of the density matrix renormalization group. _The Journal of Chemical Physics_ 2002, _116_ , 4462–4476
  * Chan 2004 Chan, G. K.-L. An algorithm for large scale density matrix renormalization group calculations. _The Journal of Chemical Physics_ 2004, _120_ , 3172–3178
  * Ghosh et al. 2008 Ghosh, D.; Hachmann, J.; Yanai, T.; Chan, G. K.-L. Orbital optimization in the density matrix renormalization group, with applications to polyenes and Î²-carotene. _The Journal of Chemical Physics_ 2008, _128_ , 144117
  * Sharma and Chan 2012 Sharma, S.; Chan, G. K.-L. Spin-adapted density matrix renormalization group algorithms for quantum chemistry. _The Journal of Chemical Physics_ 2012, _136_ , 124121
  * Sharma et al. 2017 Sharma, S.; Holmes, A. A.; Jeanmairet, G.; Alavi, A.; Umrigar, C. J. Semistochastic Heat-Bath Configuration Interaction Method: Selected Configuration Interaction with Semistochastic Perturbation Theory. _J. Chem. Theory Comput._ 2017, _13_ , 1595–1604
  * Holmes et al. 2016 Holmes, A. A.; Tubman, N. M.; Umrigar, C. J. Heat-bath Configuration Interaction: An efficient selected CI algorithm inspired by heat-bath sampling. _J. Chem. Theory Comput._ 2016, _12_ , 3674–3680
  * Smith et al. 2017 Smith, J. E. T.; Mussard, B.; Holmes, A. A.; Sharma, S. Cheap and Near Exact CASSCF with Large Active Spaces. _Journal of Chemical Theory and Computation_ 2017, _13_ , 5468–5478
  * Nielsen and Janssen 1999 Nielsen, I. M. B.; Janssen, C. L. Double-substitution-based diagnostics for coupled-cluster and Møller–Plesset perturbation theory. _Chemical physics letters_ 1999, _310_ , 568–576
  * Magoulas et al. 2021 Magoulas, I.; Gururangan, K.; Piecuch, P.; Deustua, J. E.; Shen, J. Is Externally Corrected Coupled Cluster Always Better than the Underlying Truncated Configuration Interaction? 2021

## Appendix

##  Appendix A A sweep algorithm converting MPS to CI coefficients with a
threshold

In DMRG, the electronic wave function is represented by a matrix-product state
(MPS)

$\displaystyle|\Psi\rangle$
$\displaystyle=\sum_{n_{1},\dots,n_{K}}\sum_{\alpha_{1},\dots,\alpha_{K-1}}A^{n_{1}}_{\alpha_{1}}A^{n_{2}}_{\alpha_{1}\alpha_{2}}\cdots
A^{n_{K}}_{\alpha_{K-1}}|n_{1}n_{2}\cdots n_{K}\rangle,$ (A.1)
$\displaystyle\\{n\\}$ $\displaystyle=\\{{\rm
vac},\uparrow,\downarrow,\uparrow\downarrow\\},$ (A.2)

where $n_{i}$ is the occupation of orbital $i$, $|n_{1}n_{2}\cdots
n_{K}\rangle$ is the occupation-number representation of a determinant, and
$\alpha_{i}$ is an auxiliary index. Here,
$\sum_{\alpha_{i}}A^{n_{i}}_{\alpha_{i-1}\alpha_{i}}A^{n_{i+1}}_{\alpha_{i}\alpha_{i+1}}$
denotes a matrix product and it is assumed that the dimensions of all the
auxiliary indices are the same (the bond dimension $M$). The wave function of
the ground state can be optimized by the efficient DMRG sweep algorithm.

In order to get $Q^{\rm ec}_{3}T_{3}$ and $Q^{\rm ec}_{4}T_{4}$ in Eq. (2), we
extract quadruple- and lower-order CI excitation amplitudes from the MPS by a
sweep algorithm. To avoid repeated or unnecessary computation, we here
describe how to obtain CI coefficients whose values are larger than a
threshold in Eq. (3), with a concomitant reduction in computational cost and
memory usage from a naive approach.

We first start with the MPS in left canonical form. During a sweep to compute
the amplitude, at any given point (e.g. at some site $p$) one has a set of
partial coefficients $c_{\alpha_{p}}(n_{1}n_{2}\ldots
n_{p})=\sum_{\alpha_{1}\ldots\alpha_{p-1}}A^{n_{1}}_{\alpha_{1}}\ldots
A^{n_{p}}_{\alpha_{p-1}\alpha_{p}}$. If
$\sum_{\alpha_{p}}|c_{\alpha_{p}}(n_{1}n_{2}\ldots n_{p})|^{2}<\text{thresh}$,
then this partial coefficient is dropped, as are all determinants involving
the occupancy string $|n_{1}n_{2}\ldots n_{p}\rangle$. This is because if the
MPS is in left canonical form, the above condition on the partial coefficient
guarantees that the coefficient of any determinant generated by the MPS which
contains $|n_{1}n_{2}\ldots n_{p}\rangle$ as a substring is also less than the
threshold in magnitude. In addition, since the orbitals are associated with
definite hole or particle character, we also drop any coefficient associated
with more than 4 holes or 4 particles. Finally, in this process, we can take
advantage of the conserved quantum numbers to only generate symmetry unique
partial coefficients (e.g. if $S_{z}=0$, then the values of $c_{\alpha_{p+1}}$
come in time-reversal pairs and only one needs to be considered). Thus using
the above algorithm we can completely avoid generating any determinants with
coefficients below the threshold.