# Which model density is best in pair natural orbital local correlation
theory?

Réka A. Horváth  Kesha Sorathia  Isabelle Saint  David P. Tew
<EMAIL_ADDRESS>University of Oxford, South Parks Road, Oxford, OX1
3QZ, UK

###### Abstract

Low-scaling electron correlation theory based on the pair natural orbital
approximation, PNO-CCSD(T), has become a powerful computational tool.
Motivated by the recent discovery of large errors for organometallic
molecules, we assess the role of the model density used to discard unimportant
contributions. We find that second-order perturbation theory provides the best
compromise between cost and accuracy, but coupling between localised occupied
orbitals must be accounted for. Errors in the CCSD energy are then well below
1 kcal/mol, even for molecules with moderate multi-reference character, and
the primary remaining source of errors lies in the treatment of the (T) energy
contribution.

##  I Introduction

Low-scaling electron correlation theory has undergone considerable development
over the past decade, and computational approaches have emergedRiplinger  _et
al._ (2016); Schmitz  _et al._ (2014); Ma and Werner (2018a); Nagy and Kállay
(2019) capable of accurately computing heats of reaction and transition state
energies of large molecular systems with modest computational
resources.Franzke  _et al._ Although local approximations that exploit the
short-ranged nature of electron correlation in molecules were introduced many
years ago,Sinanoǧlu (1964); Pulay (1983); Sæbø and Pulay (1985); Schütz  _et
al._ (1999) in practice, the sizable errors resulting from the domain
approximationBoughton and Pulay (1993) prevented widescale application. In
2009 Neese re-introducedNeese  _et al._ (2009) pair natural orbitial (PNO)
theory,Löwdin (1955) which made it possible to use large local domains and
eliminate instabilities arising from the domain approximation. The PNO-CCSD(T)
method (coupled cluster with singles, doubles and perturbative triples), has
emerged as practical alternative to density functional theory for
computational studies, particularly in organometalic chemistry, which requires
accurate description of exchange, delocalisation and dispersion interactions.

In PNO coupled-cluster theory, each pair-correlation function
$|\mu_{ij}\rangle$ is expanded using a pair-specific subset of virtuals, whose
number $v_{ij}$ and character are determined to recover the correlation energy
for that orbital pair to a target accuracy

$\displaystyle|\mu_{ij}\rangle=\frac{1}{2}\sum_{\tilde{a}\tilde{b}}^{v_{ij}}t^{ij}_{\tilde{a}\tilde{b}}\hat{E}_{i\tilde{a}}\hat{E}_{j\tilde{b}}|0\rangle\,.$
(1)

Here the excitation operator
$\hat{E}_{i\tilde{a}}=\hat{a}^{\dagger}_{\tilde{a}\alpha}\hat{a}_{i\alpha}+\hat{a}^{\dagger}_{\tilde{a}\beta}\hat{a}_{i\beta}$
promotes an electron from a localised occupied orbital $i$ to pair natural
virtual orbital $\tilde{a}$. $|0\rangle$ is the Hartree–Fock reference
wavefunction. Pair natural orbitals diagonalise the pair correlation density

$\displaystyle
D_{ab}=2\sum_{c}[\bar{t}^{ij}_{ac}t^{ij}_{bc}+\bar{t}^{ij}_{ca}t^{ij}_{cb}]$
(2)
$\displaystyle\sum_{ab}U^{\dagger}_{a\tilde{a}}D_{ab}U_{b\tilde{b}}=N_{\tilde{a}}\delta_{\tilde{a}\tilde{b}}$
(3)

and provide a maximally convergent expansion for the correlation
function.Löwdin and Shull (1956) Here
$\bar{t}^{ij}_{ac}=2t^{ij}_{ac}-t^{ij}_{ca}$. The size of the PNO space
$v_{ij}$ for each pair is set by retaining PNOs with
$N_{\tilde{a}}\geq\mathcal{T}$, where the PNO threshold $\mathcal{T}$
determines the accuracy of the resulting correlation energy. Formally,
$t^{ij}_{ab}$ is the amplitude for the pair excitation in canonical coupled-
cluster theory. In practice, approximate amplitudes $\tilde{t}^{ij}_{ab}$ are
used to construct a model density to determine the PNOs.Meyer (1973) The
ensuing coupled-cluster calculation is performed at low cost in the retained
subset of PNOs and the contribution to the model amplitudes from the discarded
PNOs is added as an energy correction term.Neese  _et al._ (2009) The PNO
approach is predicated on the assumption that the orbitals obtained from the
approximate amplitudes are also appropriate for representing the coupled
cluster wavefunction. This letter addresses that assumption and asks the title
question: which model density is best in pair natural orbital local
correlation theory?

The standard approach is to use first-order amplitudes from Møller–Plesset
perturbation theory to form the model density. Often a semi-canonical
approximation is used that neglects off-diagonal Fock matrix elements
$f_{i}^{j}$ among the localised occupied orbitals.

$\displaystyle\tilde{t}^{ij}_{ab}=-\frac{g^{ij}_{ab}}{\epsilon_{a}+\epsilon_{b}-f^{i}_{i}-f^{j}_{j}}$
(4)

where $g^{ij}_{ab}=\langle ia|jb\rangle$ are the electron repulsion integrals
and $\epsilon_{a}$ the canonical orbital eigenvalues. Generally, the energy
from a PNO-CCSD(T) calculation $E(\mathcal{T})$ converges to the canonical
result $E$ as the threshold $\mathcal{T}$ is reduced, and the size of the PNO
space increases, according to the power lawSorathia and Tew (2020); Sorathia
_et al._ (2023)

$\displaystyle E(\mathcal{T})=E+A\mathcal{T}^{1/2}$ (5)

To obtain accurate energies for the reliable prediction of chemical phenomena,
it is necessary to ensure that the PNO truncation error $A\mathcal{T}^{1/2}$
is under control. Neese and co-workers demonstratedLiakos  _et al._ (2020)
that mean absolute PNO truncation errors with $\mathcal{T}=10^{-7}$ are around
0.5 kcal/mol or less in PNO-CCSD(T) calculations for the GMTKN55 data
set,Goerigk  _et al._ (2017) which includes reaction energies and
isomerisation energies of small and large molecules, intermolecular and
intramolecular non-covalent interaction energies and reaction barrier heights.

However, for systems with a higher proportion of static correlation, PNO
truncation errors can be several kcal/mol, and $\mathcal{T}=10^{-8}$ or PNO
extrapolation is necessary.Altun  _et al._ (2020); Sorathia  _et al._ (2023)
This is evident from the work of Martin and co-workersSemidalas and Martin
(2022) on metal organic barrier heightsIron and Janes (2019) and expanded
porphyrin macrocycles that can interconvert between Hückel, figure-of-eight,
and Möbius topologies.Sylvetsky  _et al._ (2020)

The large PNO truncation errors are an indication that the MP2 pair
correlation density may fail to accurately model the CCSD pair correlation
density for these systems, reducing the proportion of correlation energy
recovered in the selected PNO subspace and reducing the accuracy of the energy
estimate for the discarded PNOs.

We have therefore implemented the possibility, within the `pnoccsd` module in
the Turbomole program package,TUR  to perform a hierarchy of intermediate
methods between MP2 theory and CCSD theory, to obtain the corresponding model
density, and to perform PNO-CCSD(T) calculations using PNOs obtained from any
of these model densities. The accuracy of each model can be assessed through
the size of the PNO truncation error for computed PNO-CCSD(T) correlation
energies.

If accuracy were the only criteria, then the answer to the title question
would be trivial. The most accurate model density is obviously that obtained
from the canonical CCSD amplitudes, but this is of no practical use when
designing low-cost correlation methods. “Best” here refers to the best
compromise between accuracy and the cost incurred to obtain the model density.

The doubles amplitude equations for our chosen models are

$\displaystyle 0$ $\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\tilde{H}+[\tilde{H},\hat{T}_{2}]+\tfrac{1}{2}[[\tilde{H},\hat{T}_{2}],\hat{T}_{2}]|0\rangle\quad\textrm{CCSD}$
(6) $\displaystyle 0$ $\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\hat{H}+[\hat{H},\hat{T}_{2}]+\tfrac{1}{2}[[\hat{H},\hat{T}_{2}],\hat{T}_{2}]|0\rangle\quad\textrm{CCD}$
(7) $\displaystyle 0$ $\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\tilde{H}+[\tilde{H},\hat{T}_{2}]|0\rangle\qquad\qquad\qquad\qquad\textrm{LCCSD}$
(8) $\displaystyle 0$ $\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\hat{H}+[\hat{H},\hat{T}_{2}]|0\rangle\qquad\qquad\qquad\qquad\textrm{LCCD}$
(9) $\displaystyle 0$ $\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\tilde{H}+[\hat{F},\hat{T}_{2}]|0\rangle\qquad\qquad\qquad\qquad\,\textrm{CC2}$
(10) $\displaystyle 0$
$\displaystyle=\langle\overline{\begin{subarray}{c}ab\\\
ij\end{subarray}}|\hat{H}+[\hat{F},\hat{T}_{2}]|0\rangle\qquad\qquad\qquad\qquad\,\textrm{MP2}$
(11)

where $\tilde{H}$ is the $T_{1}$ transformed Hamiltonian
$e^{-\hat{T}_{1}}\hat{H}e^{\hat{T}_{1}}$ and
$\langle\overline{\begin{subarray}{c}ab\\\ ij\end{subarray}}|$ are the
biorthogonal bra states for the doubles projection manifold. The methods CCD
and LCCD are obtained from CCSD and LCCSD, respectively, by setting all
$T_{1}$ amplitudes to zero so that $\tilde{H}=\hat{H}$. The singles amplitude
equations for CC2 and LCCSD are the same as for CCSD

$\displaystyle 0$ $\displaystyle=\langle\begin{subarray}{c}a\\\
i\end{subarray}|\tilde{H}+[\tilde{H},\hat{T}_{2}]|0\rangle$ (12)

To assess the accuracy of each model, we chose to perform calculations on the
MOBH35 data set of barrier heights for organometallic reactionsIron and Janes
(2019) and the POLYPYR21 data set of Hükel-Möbius porphyrin
macrocycles,Sylvetsky  _et al._ (2020) which contain species with moderate
amounts of static correlation and large PNO truncation errors. Although
triple-zeta quality basis sets are the minimum required to obtain reliable
CCSD(T) energies for comparison with experiment, the behaviour of the PNO
truncation error is very similar across basis sets.Liakos  _et al._ (2020);
Sorathia  _et al._ (2023) We therefre use double-zeta quality basis sets in
this work.

Two measures of accuracy are used to assess each model density. The first,
$Q_{d}$, is a direct measure of the proportion of the CCSD density discarded
when truncating the space using PNOs from a model density; the second,
$Q_{e}$, is the error in the correlation energy incurred when the selected
model PNOs are used in a PNO-CCSD(T) calculation.

We define $Q_{d}(\mathcal{T})$ as the average density lost per pair using PNO
truncation threshold $\mathcal{T}$

$\displaystyle Q_{d}(\mathcal{T})$ $\displaystyle=\frac{1}{n_{ij}}\sum_{i\leq
j}(\mathrm{Tr}[\mathbf{D}^{ij}]-\mathrm{Tr}[\tilde{\mathbf{D}}^{ij}(\mathcal{T})])/\mathrm{Tr}[\mathbf{D}^{ij}]$
(13) $\displaystyle\tilde{D}^{ij}_{\tilde{a}\tilde{b}}$
$\displaystyle=\sum_{ab}^{v}U^{ij}_{a\tilde{a}}D^{ij}_{ab}U^{ij}_{b\tilde{b}}$
(14)

Here the CCSD pair correlation density $\mathbf{D}^{ij}$ formally has the
$v\times v$ dimension of the full virtual space, $\tilde{\mathbf{D}}^{ij}(T)$
is the truncated pair density that has the $v_{ij}\times v_{ij}$ dimension of
the retained PNOs at treshold $\mathcal{T}$, and $U^{ij}_{a\tilde{a}}$ is the
transformation matrix from $v$ canonical virtuals to the $v_{ij}$ PNOs for
pair $ij$ obtained from the chosen model density at threshold $\mathcal{T}$.
By plotting $Q_{d}(\mathcal{T})$ against the average dimension
$v_{ij}(\mathcal{T})$ as a function of $\mathcal{T}$ for each model density,
we can compare the relative accuracy of each model.

We define $Q_{e}(\mathcal{T})$ as the error in the CCSD(T) correlation energy
using the PNOs from the model density and truncation threshold $\mathcal{T}$.
It is composed of three terms

$\displaystyle
Q_{e}(\mathcal{T})=\mathcal{E}_{\text{SD}}(\mathcal{T})+\mathcal{E}_{\text{T}}(\mathcal{T})-\Delta(\mathcal{T})$
(15)

$\mathcal{E}_{\text{SD}}$ is the error in the CCSD correlation energy and
$\mathcal{E}_{\text{T}}$ is the error in the (T) energy. The
$\mathcal{E}_{\text{T}}$ contribution depends primarily on the truncation of
the triples space through the selection of triple natural orbitals,Riplinger
_et al._ (2013) which is independent of the chosen pair density model in the
current implementation.Schmitz and Hättig (2016) Variations in
$\mathcal{E}_{\text{T}}$ among model densities are a reflection of the
accuracy of the $T_{2}$ amplitudes used to compute the (T) energy.

The $\Delta$ term corrects for $\mathcal{E}_{\text{SD}}$ and is the energy
estimate for the discarded PNOs based on the model amplitudes

$\displaystyle\Delta=\sum_{ij}\left(\sum_{ab}^{v}\tilde{t}^{ij}_{ab}\bar{g}_{ij}^{ab}-\sum_{\tilde{a}\tilde{b}}^{v_{ij}}\tilde{t}^{ij}_{\tilde{a}\tilde{b}}\bar{g}_{ij}^{\tilde{a}\tilde{b}}\right)$
(16)

Here $\bar{g}_{ij}^{ab}=2g_{ij}^{ab}-g_{ij}^{ba}$. We report error statistics
in units of mEh per valence electron for total energies, and in kcal/mol for
reaction energies and barrier heights.

Figure 1:  The loss of information against average number of PNOs for
different model densities for reactant 29 of the MOBH35 set

##  II Computational details

The structures for the MOBH35 test set were taken from the supporting
information of Ref Semidalas and Martin, 2022, where the transition state
structures for reactions 11 and 12, and all species of reaction 14 are
modified from the original data set, as recommended by Dohm _et al_.Dohm  _et
al._ (2020) We use the def2-SVPWeigend and Ahlrichs (2005) basis and density
fitting basisWeigend (2006) and the Stuttgart relativistic effective core
potentialsDolg and Cao (2012) for molecules containing second- and third-row
transition metal atoms. Reactant 16 converges to the incorrect state if the
default extended Hückel orbital guess is applied and an orbital swap was
required to obtain the correct initial electronic configuration. Canonical
CCSD(T) energies for 19 of the 35 reactions are available from earlier
work,Sorathia  _et al._ (2023) so we use this subset MOBH19 for the current
study. These are reactions 3, 4, 6, 7, 13, 14, 15, 16, 21, 23, 26, 27, 29, 30,
31, 32, 33, 34 and 35.

The structures for the POLYPYR21 data set were taken from the supporting
information of Ref. Sylvetsky  _et al._ , 2020. To enable direct comparison
with the CCSD and CCSD(T) energies of the original work, we use the cc-pVDZ
basis setDunning (1989) with the $d$ functions on the H atoms removed, and the
cc-pVDZ density fitting basis.Weigend  _et al._ (2002)

All calculations were performed using the the `dscf`Häser and Ahlrichs  and
`pnoccsd` modules,Schmitz  _et al._ (2014); Schmitz and Hättig (2016); Tew
(2019) of the Turbomole package. Density fitting was not used in the
Hartree–Fock calculations. Orbitals were localised using the intrinsic bond
orbital method.Knizia (2013) and in all cases the multilevel
approximationMasur  _et al._ (2013); Schütz  _et al._ (2014) were not applied.
The $Q_{d}$ and $Q_{e}$ measures were computed from near-canonical model
densities obtained from the corresponding PNO coupled cluster calculation
using a very tight PNO threshold of $\mathcal{T}$=$10^{-9}$. Where timings are
reported, these were performed on a single Intel(R) Xeon(R) Gold 6248R CPU @
3.00GHz node with 48 cores, 380 Gb RAM and 1.8Tb SSD.

##  III Results

Figs. 1 and 2 display the average information lost $Q_{d}(\mathcal{T})$
against the average number of PNOs $v_{ij}(\mathcal{T})$ retained with
$\mathcal{T}$=$10^{-5}$–$10^{-8}$ for the hierarchy of model densities. The
two chosen molecules r29 and r32 are the extremal cases of the MOBH19 set
where the SCMP2 density performs best and worst relative to CCSD,
respectively. The CCSD model loses the least possible information per PNO
discarded and is the reference against which the other methods are judged.

Figure 2:  The loss of information against average number of PNOs for
different model densities for reactant 32 of the MOBH35 set

The LCCSD and LCCD models are only marginally inferior to CCSD and CCD,
respectively. The influence of the quadratic $T_{2}$ terms in the CC model are
therefore small on average. We find that they are more significant for strong
pairs than weak pairs, but in most cases have less impact than the singles
cluster amplitudes. Particulary for molecules with substantial correlation-
induced orbital relaxation, the CCD and LCCD models lose accuracy relative to
CCSD and CC2 is more accurate than MP2. In almost all cases MP2 is
significantly better than semi-canonical MP2. Neglecting the Fock matrix
elements that couple the localised orbitals is a poor approximation. Overall,
the accuracy of the CC model densities is in line with that expected from a
perturbative analysis of the the amplitude equations, and the magnitude of the
interaction follows $f_{i}^{j}>[\hat{H},T_{2}]>T_{1}>[[\hat{H},T_{2}],T_{2}]$.

Figure 3:  Mean absolute PNO truncation errors $\mathcal{E}_{\text{SD}}$
($\times$) and $\mathcal{E}_{\text{SD}}-\Delta$ ($\square$) against the
average number of PNOs for the MOBH19 set using various model densities

Fig. 3 shows the mean absolute truncation error per valence electron in the
CCSD correlation energy against the mean number of PNOs per pair, averaged
over the dat set, for $\mathcal{T}$=$10^{-5}$–$10^{-8}$. For each model
density, we plot the mean absolute values of $\mathcal{E}_{\text{SD}}$ and
$\mathcal{E}_{\text{SD}}-\Delta$, to ascertain the relative accuracy of the
energy estimate for the discarded PNOs. The mean $\mathcal{E}_{\text{SD}}$
values for each method mirror the trend observed for $Q_{d}$, with
CCSD$\sim$LCCSD$<$CCD$\sim$LCCD$<$CC2$\sim$MP2$<$SCMP2.

Figure 4:  Box plots for PNO-CCSD errors $\mathcal{E}_{\text{SD}}$ (left) and
$\mathcal{E}_{\text{SD}}-\Delta$ (right) in MOBH19 barrier heights using
$\mathcal{T}$=$10^{-6}$.

A greater variation in accuracy is observed for
$\mathcal{E}_{\text{SD}}-\Delta$, but again the same trend is observed.
However, we find that the CC2 model is significantly more accurate than MP2
when the $\Delta$ correction is taken into account. Interestingly, for CCSD,
the $\Delta$ correction does not entirely eliminate the PNO trunctation error,
and a small error remains due to relaxation of the $T_{2}$ amplitudes when
solving the non-linear CCSD equations in the retained subspace.

Table 1: PNO truncation error normal distributions $(\mu,\sigma_{n})$ for barrier heights of the MOBH19 set in kcal/mol |  |  |  | SCMP2 |  | MP2 |  | CC2 |  | LCCD |  | CCD |  | LCCSD |  | CCSD
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Contr. |  | $\mathcal{T}$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$ |  | $\mu$ | $\sigma$
$\mathcal{E}_{\text{SD}}$ |  | $10^{-5}$ |  | 1.23 | 2.06 |  | 0.94 | 1.72 |  | 0.78 | 1.50 |  | 1.12 | 1.99 |  | 1.13 | 2.03 |  | 0.98 | 1.79 |  | 1.03 | 1.78
|  | $10^{-6}$ |  | 0.37 | 0.74 |  | 0.27 | 0.61 |  | 0.20 | 0.46 |  | 0.30 | 0.64 |  | 0.29 | 0.65 |  | 0.24 | 0.50 |  | 0.23 | 0.48
|  | $10^{-7}$ |  | 0.12 | 0.24 |  | 0.08 | 0.19 |  | 0.04 | 0.12 |  | 0.09 | 0.20 |  | 0.09 | 0.20 |  | 0.06 | 0.12 |  | 0.06 | 0.12
|  | $10^{-8}$ |  | 0.03 | 0.09 |  | 0.02 | 0.07 |  | 0.01 | 0.05 |  | 0.03 | 0.07 |  | 0.02 | 0.07 |  | 0.02 | 0.05 |  | 0.01 | 0.04
$\mathcal{E}_{\text{SD}}-\Delta$ |  | $10^{-5}$ |  | 0.20 | 0.73 |  | 0.08 | 0.48 |  |  $-$0.03 | 0.34 |  | 0.16 | 0.55 |  | 0.16 | 0.55 |  | 0.09 | 0.28 |  | 0.04 | 0.19
|  | $10^{-6}$ |  | 0.15 | 0.38 |  | 0.08 | 0.27 |  | 0.02 | 0.14 |  | 0.10 | 0.30 |  | 0.09 | 0.30 |  | 0.05 | 0.14 |  | 0.04 | 0.10
|  | $10^{-7}$ |  | 0.08 | 0.18 |  | 0.03 | 0.12 |  | 0.01 | 0.06 |  | 0.04 | 0.13 |  | 0.04 | 0.13 |  | 0.02 | 0.06 |  | 0.01 | 0.05
|  | $10^{-8}$ |  | 0.02 | 0.08 |  | 0.01 | 0.06 |  | 0.01 | 0.04 |  | 0.02 | 0.06 |  | 0.01 | 0.07 |  | 0.01 | 0.04 |  | 0.01 | 0.04
$\mathcal{E}_{\text{(T)}}$ |  | $10^{-5}$ |  | 0.78 | 1.41 |  | 0.62 | 1.09 |  | 0.61 | 1.05 |  | 0.62 | 1.14 |  | 0.64 | 1.17 |  | 0.64 | 1.09 |  | 0.64 | 1.09
|  | $10^{-6}$ |  | 0.41 | 0.73 |  | 0.28 | 0.55 |  | 0.28 | 0.51 |  | 0.29 | 0.58 |  | 0.30 | 0.59 |  | 0.29 | 0.53 |  | 0.29 | 0.54
|  | $10^{-7}$ |  | 0.16 | 0.31 |  | 0.11 | 0.24 |  | 0.11 | 0.22 |  | 0.12 | 0.25 |  | 0.12 | 0.26 |  | 0.11 | 0.22 |  | 0.11 | 0.22
|  | $10^{-8}$ |  | 0.05 | 0.11 |  | 0.04 | 0.09 |  | 0.03 | 0.08 |  | 0.04 | 0.09 |  | 0.04 | 0.10 |  | 0.04 | 0.08 |  | 0.04 | 0.09
$Q_{e}$ |  | $10^{-5}$ |  | 0.98 | 1.84 |  | 0.70 | 1.34 |  | 0.58 | 1.06 |  | 0.78 | 1.52 |  | 0.80 | 1.56 |  | 0.72 | 1.21 |  | 0.68 | 1.14
|  | $10^{-6}$ |  | 0.56 | 1.06 |  | 0.36 | 0.77 |  | 0.30 | 0.59 |  | 0.38 | 0.82 |  | 0.39 | 0.84 |  | 0.34 | 0.63 |  | 0.33 | 0.60
|  | $10^{-7}$ |  | 0.24 | 0.47 |  | 0.14 | 0.34 |  | 0.12 | 0.26 |  | 0.16 | 0.36 |  | 0.16 | 0.37 |  | 0.13 | 0.26 |  | 0.13 | 0.26
|  | $10^{-8}$ |  | 0.08 | 0.19 |  | 0.05 | 0.14 |  | 0.04 | 0.11 |  | 0.05 | 0.15 |  | 0.05 | 0.15 |  | 0.04 | 0.11 |  | 0.05 | 0.11

In Table 1 we report statistical measures for the distribution of PNO
truncation errors in the barrier heights of the MOBH19 set. From the data for
the CCSD model, we see that for both $\mathcal{E}_{\text{SD}}$ and
$\mathcal{E}_{\text{(T)}}$ the $2\sigma$ 95% confidence interval falls within
1 kcal/mol only for $\mathcal{T}$=$10^{-7}$ and tighter. The
$\mathcal{E}_{\text{SD}}$ contribution is almost entirely removed by the
$\Delta$ correction and the error in the (T) energy becomes the limiting
factor. All of the approximate models follow this same pattern.

Fig. 4 displays box plots for the PNO truncation errors in the PNO-CCSD
barrier heights of the MOBH19 set, with and without the correction for
discarded PNOs, using the various model densities with
$\mathcal{T}$=$10^{-6}$. Fig. 5 displays the corresponding information for
$\mathcal{T}$=$10^{-7}$. To indicate the accuracy relative to the subspace
size, the data for each method is positioned along the $x$ axis at the average
of $\langle v_{ij}\rangle$ for the transition states.

Figure 5:  Box plots for PNO-CCSD errors $\mathcal{E}_{\text{SD}}$ (left) and
$\mathcal{E}_{\text{SD}}-\Delta$ (right) in MOBH19 barrier heights using
$\mathcal{T}$=$10^{-7}$.

From the table and the figures, we see that for a given threshold
$\mathcal{T}$, the CC2 model gives the smallest errors, but this is because it
predicts a larger density and therefore more PNOs are retained. Nevertheless,
the $\Delta$ correction with the CC2 model appears to be particularly
accurate. The ratio of accuracy to number of PNOs for barrier heights using
the model densities follows the same trend as that seen for total energies and
information content. The MP2 model is remarkably accurate overall and the
benefit of using higher-order CC models appears marginal, especially since the
(T) error is the dominant contribution.

The $\mathcal{E}_{\text{(T)}}$ statistics are very similar for all methods.
The primary factor in $\mathcal{E}_{\text{(T)}}$ is the model density used for
the TNO truncation, which is independent of the PNO model density in this
work. The error in the (T) energy due to errors in the $T_{2}$ amplitudes are
smaller, but not insignificant, and give rise to some variation among the
model densities tested. We find that the maximum errors are generally smaller
for CC2 and LCCSD than for LCCD and CCD. Overall, the SCMP2 performs the
worst, and the majority of the available improvement is obtained already at
the MP2 level. Using MP2 or CC2 model densities with $\mathcal{T}$=$10^{-7}$
is sufficient to reduce the maximum total PNO truncation errors to within 1
kcal/mol.

The cost of using a higher-order CC model to obtain an approximate density has
two sources: the overhead of performing the higher-order CC model using a
tight threshold; and a cost or benefit arising from a larger or smaller PNO
set for the subsequent PNO-CCSD(T) calculation. The overhead can be reduced by
loosening the PNO threshold used to obtain the model density.

Table 2: Mean absulte (and maximum) PNO-CCSD errors for MOBH19 barriers in kcal/mol Density | $\mathcal{T}=10^{-5}$ | $\mathcal{T}=10^{-6}$ | $\mathcal{T}=10^{-7}$
---|---|---|---
SCMP2 $(10^{-9})$ | $0.51$ | $(2.61)$ | $0.28$ | $(1.47)$ | $0.13$ | $(0.78)$
MP2 $(0.1\mathcal{T})$ | $0.38$ | $(1.45)$ | $0.22$ | $(1.20)$ | $0.09$ | $(0.54)$
MP2 $(10^{-9})$ | $0.31$ | $(1.80)$ | $0.18$ | $(1.08)$ | $0.08$ | $(0.48)$
CC2 $(0.1\mathcal{T})$ | $0.26$ | $(0.78)$ | $0.14$ | $(0.80)$ | $0.06$ | $(0.34)$
CC2 $(10^{-9})$ | $0.23$ | $(1.11)$ | $0.08$ | $(0.43)$ | $0.04$ | $(0.21)$

In Table 2 we list the mean absolute and maximum errors using the MP2 and CC2
model densities computed using a threshold of $0.1\mathcal{T}$, compared to
the near-canonical densities obtained with $10^{-9}$. The errors obtained
using the SCMP2 density are included for comparison. Evidently, using a
threshold of $0.1\mathcal{T}$ is sufficient to recover more than two-thirds of
the available improvement over SCMP2 on average, but the outliers are
significantly improved if the $\Delta$ correction spans the larger PNO space.

Figure 6:  Box plots for PNO-CCSD(T) errors of the POLYPYR21 set using
different model densities with $\mathcal{T}$=$10^{-6}$, $10^{-7}$ and
$10^{-8}$.

In Fig. 6 we report the error statistics for the POLYPYR21 set using
$\mathcal{T}$=$10^{-6}$, $10^{-7}$ and $10^{-8}$ and a threshold of
$0.1\mathcal{T}$ for the model density. To provide a direct measure of the
cost, we compute a weighted sum of the wall times in seconds for the test set,
where the weight for each molecule is the reciprocol of the size of the AO
basis. The box plot for each model and threshold is positioned along the x
axis at the value of the corresponding cost measure.

The PNO truncation errors for $\mathcal{T}$=$10^{-6}$ are up to 8 kcal/mol.
The PNO approximation with this threshold introduces an error so large that
almost no benefit from the sophisticated coupled-cluster wavefunction model
remains. Even with $\mathcal{T}$=$10^{-7}$, errors are as much as 4 kcal/mol
for the Möbius structures. With the SCMP2 model density and
$\mathcal{T}$=$10^{-8}$, some PNO errors of 2 kcal/mol remain. However, if the
MP2, CC2 or LCCSD model densities are used, errors approach the reliability
target of 1 kcal/mol. Among these, the MP2 model is the cheapest.

The timings used for the figure are only indicative, since the implementation
of each approach has not been fully optimised. Nevertheless, we note some
interesting trends in the relative costs of using the model densities. For
example, LCCD is more expensive than CCD. This is because, although each
iteration is cheaper, the LCCD equations take longer to converge. The CC2
model is much more expensive than MP2, both because the number of retained
PNOs is higher, and because the singles residual contains long-ranged Fock-
like terms that are expensive to evaluate. The cost of the MP2 model can
probably be reduced to close to that of the SCMP2 model by use of the Laplace
transform approach.Häser and Almlöf (1992)

Table 3: PNO truncation errors for CCSD energies of the POLYPYR21 set in
kcal/mol

System | Molpro | Orca | MRCC | SCMP2a | MP2a,b
---|---|---|---|---|---
24H${}_{\text{a}}$ | $0.1$ | $0.0$ | $0.0$ | $0.0$ | $0.0$
24H${}_{\text{b}}$ | $0.0$ | $0.0$ | $0.0$ | $0.0$ | $0.0$
24M | $0.3$ | $0.1$ | $0.1$ | $0.1$ | $0.1$
24TS1 | $-0.1$ | $0.0$ | $0.1$ | $0.0$ | $0.0$
24TS2 | $0.1$ | $0.0$ | $0.1$ | $0.1$ | $0.0$
28F | $-0.2$ | $-0.6$ | $-0.5$ | $-0.2$ | $0.0$
28M${}_{1\text{A}}$ | $2.0$ | $0.5$ | $0.3$ | $0.3$ | $0.3$
28M | $1.8$ | $0.7$ | $0.2$ | $0.4$ | $0.4$
28M${}_{1\text{B}}$ | $2.2$ | $0.5$ | $0.2$ | $0.4$ | $0.4$
28TS3 | $0.0$ | $-0.1$ | $-0.4$ | $0.0$ | $0.0$
28H | $0.0$ | $0.0$ | $0.0$ | $0.0$ | $0.0$
28TS${}_{1\text{A}}$ | $0.2$ | $0.3$ | $0.3$ | $0.2$ | $0.1$
28TS${}_{1\text{B}}$ | $-0.1$ | $-0.1$ | $-0.4$ | $0.0$ | $0.0$
28TS${}_{2\text{A}}$ | $2.4$ | $1.0$ | $0.2$ | $0.6$ | $0.4$
28TS${}_{2\text{B}}$ | $2.0$ | $0.8$ | $0.2$ | $0.6$ | $0.5$
32F | $0.0$ | $0.0$ | $0.0$ | $0.0$ | $0.0$
32H | $-0.7$ | $0.5$ | $0.4$ | $0.3$ | $-0.1$
32TS2 | $-0.7$ | $0.5$ | $0.3$ | $0.2$ | $-0.1$
32M${}_{\text{a}}$ | $1.4$ | $0.8$ | $1.2$ | $0.8$ | $0.4$
32M${}_{\text{b}}$ | $2.2$ | $1.2$ | $1.0$ | $0.9$ | $0.6$
32TS1 | $-0.3$ | $0.1$ | $0.1$ | $0.1$ | $0.0$
RMS | $1.3$ | $0.6$ | $0.5$ | $0.4$ | $0.3$

a Turbomole b This work

In Table 3 we report the PNO truncation errors using the SCMP2 and MP2 model
densities for PNO-CCSD energies of the POLYPYR21 set, compared to those
reported in Ref. Sylvetsky  _et al._ , 2020 using the PNO-CCSD implementations
in ORCAGuo  _et al._ (2018) and MolproMa and Werner (2018b), and using the
LNO-CCSD implementation in MRCC.Nagy  _et al._ (2018) In all cases, a
threshold corresponding to $\mathcal{T}$=$10^{-7}$ is used. The differences
among the Molpro, Orca and Turbomole PNO implementations are most likely due
to the differing domain sizes and multi-level pair approximations. Using the
full MP2 density instead of the semi-canonical density to determine the PNOs
ensures that the root mean squared (RMS) CCSD error is 0.3 kcal/mol with a
maximum error of only 0.6 kcal/mol for the challenging Möbius structures.

The truncation error in the PNO-CCSD(T) energies, in Fig. 6 is dominated by
the (T) energy contribution. The RMS error rises from 0.3 kcal/mol for CCSD to
1.6 kcal/mol for CCSD(T), when using the MP2 model density for the PNO space.
The maximum error rises from 0.6 to 3.0 kcal/mol. The truncation error can be
reduced through two-point extrapolation assuming Eq. 5. Application of two-
point extrapolation of PNO-CCSD(T) energies using $\mathcal{T}$=$10^{-6}$ and
$\mathcal{T}$=$10^{-7}$ data reduces RMS and maximum errors to 0.9 kcal/mol
and 1.7 kcal/mol, respectively. Neverthelss, the primary remaining bottleneck
in accuracy for PNO calculations on systems with moderate static correlation
is treatment of the (T) energy contribution.

##  IV Conclusions

We set out to investigate the extent to which the model density used to
construct the PNOs impacts the PNO truncation errors in local correlation
theory and to answer the title question ”which model density is best?” by
assessing the relative cost-benefit ratios of a hierarchy of model densities,
starting with semi-canonical MP2 theory and ending with full CCSD. We were
motivated to address this question by the observation that large truncation
errors are observed for transition metal containing systems with moderate
static correlation, where canonical coupled-cluster theory is expected to
provide accurate predictions.Tew (2016); Giner  _et al._ (2018); Li Manni  _et
al._ (2019)

We find that even when using the best possible model density, CCSD itself, the
energy converges slowly with PNO truncation threshold for systems with a high
degree of multi-reference character. Provided that the Fock terms $f_{i}^{j}$
are not neglected in the model, there is only a modest improvement of the PNO
space when using CCSD instead of MP2 for these systems. However, the quality
of the model density does greatly impact the accuracy of the energy correction
for discarded PNOs. We find that coupling to the singles amplitudes is more
important than higher-order doubles terms in the CC model.

In the context of practical low-scaling electronic structure approaches, the
two candidates for best cost-to-accuracy ratios are therefore the MP2 and CC2
models. Due to the significant additional cost of computing the singles
residual for CC2, our recommendation is to use the MP2 model. The application
of the semi-canonical approximation leads to a substantial deterioration of
both the PNOs and the energy correction term and should be avoided.

Using the MP2 density, errors in the CCSD energy are well below 1 kcal/mol
with $\mathcal{T}$=$10^{-7}$ and the primary remaining error arises from the
truncation of the triples space when computing the (T) energy. Although not
tested in this work, these conclusions are equally applicable to the
explicitly-correlated variants of local correlation theory.Tew  _et al._
(2011); Schmitz  _et al._ (2014); Tew (2021); Ma  _et al._ (2017); Pavošević
_et al._ (2017)

## References

  * Riplinger  _et al._ (2016) Christoph Riplinger, Peter Pinski, Ute Becker, Edward F. Valeev, and Frank Neese, “Sparse maps—A systematic infrastructure for reduced-scaling electronic structure methods. II. Linear scaling domain based pair natural orbital coupled cluster theory,” J. Chem. Phys. 144, 024109 (2016).
  * Schmitz  _et al._ (2014) Gunnar Schmitz, Christof Hättig, and David P. Tew, “Explicitly correlated pno-mp2 and pno-ccsd and their application to the s66 set and large molecular systems,” Phys. Chem. Chem. Phys. 16, 22167–22178 (2014).
  * Ma and Werner (2018a) Qianli Ma and Hans-Joachim Werner, “Explicitly correlated local coupled-cluster methods using pair natural orbitals,” WIREs Comput. Mol. Sci. 8, e1371 (2018a).
  * Nagy and Kállay (2019) Péter R. Nagy and Mihály Kállay, “Approaching the basis set limit of ccsd(t) energies for large molecules with local natural orbital coupled-cluster methods,” J. Chem. Theory Comput. 15, 5275–5298 (2019).
  * Franzke  _et al._ (0) Yannick J. Franzke, Christof Holzer, Josefine H. Andersen, Tomislav Begušić, Florian Bruder, Sonia Coriani, Fabio Della Sala, Eduardo Fabiano, Daniil A. Fedotov, Susanne Fürst, Sebastian Gillhuber, Robin Grotjahn, Martin Kaupp, Max Kehry, Marjan Krstić, Fabian Mack, Sourav Majumdar, Brian D. Nguyen, Shane M. Parker, Fabian Pauly, Ansgar Pausch, Eva Perlt, Gabriel S. Phun, Ahmadreza Rajabi, Dmitrij Rappoport, Bibek Samal, Tim Schrader, Manas Sharma, Enrico Tapavicza, Robert S. Treß, Vamsee Voora, Artur Wodyński, Jason M. Yu, Benedikt Zerulla, Filipp Furche, Christof Hättig, Marek Sierka, David P. Tew, and Florian Weigend, “Turbomole: Today and tomorrow,” J. Chem. Theory Comput. 0, null (0).
  * Sinanoǧlu (1964) Oktay Sinanoǧlu, “Many-electron theory of atoms, molecules and their interactions,” in  _Adv. Chem. Phys._ (John Wiley & Sons, Ltd, 1964) pp. 315–412.
  * Pulay (1983) Peter Pulay, “Localizability of dynamic electron correlation,” Chem. Phys. Lett. 100, 151 – 154 (1983).
  * Sæbø and Pulay (1985) Svein Sæbø and Peter Pulay, “Local configuration interaction: An efficient approach for larger molecules,” Chem. Phys. Lett. 113, 13–18 (1985).
  * Schütz  _et al._ (1999) Martin Schütz, Georg Hetzer, and Hans Joachim Werner, “Low-order scaling local electron correlation methods. I. Linear scaling local MP2,” J. Chem. Phys. 111, 5691–5705 (1999).
  * Boughton and Pulay (1993) James W. Boughton and Peter Pulay, “Comparison of the boys and pipek–mezey localizations in the local correlation approach and automatic virtual basis selection,” J. Comput. Chem. 14, 736–740 (1993).
  * Neese  _et al._ (2009) Frank Neese, Frank Wennmohs, and Andreas Hansen, “Efficient and accurate local approximations to coupled-electron pair approaches: An attempt to revive the pair natural orbital method,” J. Chem. Phys. 130, 114108 (2009).
  * Löwdin (1955) Per-Olov Löwdin, “Quantum theory of many-particle systems. i. physical interpretations by means of density matrices, natural spin-orbitals, and convergence problems in the method of configurational interaction,” Phys. Rev. 97, 1474–1489 (1955).
  * Löwdin and Shull (1956) Per-Olov Löwdin and Harrison Shull, “Natural orbitals in the quantum theory of two-electron systems,” Phys. Rev. 101, 1730–1739 (1956).
  * Meyer (1973) Wilfried Meyer, “PNO–CI Studies of electron correlation effects. I. Configuration expansion by means of nonorthogonal orbitals, and application to the ground state and ionized states of methane,” J. Chem. Phys. 58, 1017–1035 (1973).
  * Sorathia and Tew (2020) Kesha Sorathia and David P. Tew, “Basis set extrapolation in pair natural orbital theories,” J. Chem. Phys. 153, 174112 (2020).
  * Sorathia  _et al._ (2023) Kesha Sorathia, Damyan Frantzov, and David P. Tew, “Improved cps and cbs extrapolation of pno-ccsd(t) energies: The mobh35 and isol24 data sets,” (2023), arXiv:2309.02639 [physics.chem-ph] .
  * Liakos  _et al._ (2020) Dimitrios G. Liakos, Yang Guo, and Frank Neese, “Comprehensive benchmark results for the domain based local pair natural orbital coupled cluster method (dlpno-ccsd(t)) for closed- and open-shell systems,” J. Phys. Chem. A 124, 90–100 (2020).
  * Goerigk  _et al._ (2017) Lars Goerigk, Andreas Hansen, Christoph Bauer, Stephan Ehrlich, Asim Najibi, and Stefan Grimme, “A look at the density functional theory zoo with the advanced gmtkn55 database for general main group thermochemistry, kinetics and noncovalent interactions,” Phys. Chem. Chem. Phys. 19, 32184–32215 (2017).
  * Altun  _et al._ (2020) Ahmet Altun, Frank Neese, and Giovanni Bistoni, “Extrapolation to the limit of a complete pair natural orbital space in local coupled-cluster calculations,” J. Chem. Theory Comput. 16, 6142–6149 (2020).
  * Semidalas and Martin (2022) Emmanouil Semidalas and Jan M.L. Martin, “The mobh35 metal–organic barrier heights reconsidered: Performance of local-orbital coupled cluster approaches in different static correlation regimes,” J. Chem. Theory Comput. 18, 883–898 (2022).
  * Iron and Janes (2019) Mark A. Iron and Trevor Janes, “Evaluating transition metal barrier heights with the latest density functional theory exchange–correlation functionals: The mobh35 benchmark database,” J. Phys. Chem. A 123, 3761–3781 (2019).
  * Sylvetsky  _et al._ (2020) Nitai Sylvetsky, Ambar Banerjee, Mercedes Alonso, and Jan M. L. Martin, “Performance of localized coupled cluster methods in a moderately strong correlation regime: Hückel–möbius interconversions in expanded porphyrins,” J. Chem. Theory Comput. 16, 3641–3653 (2020).
  * (23) “TURBOMOLE V7.7 2023, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989–2007, TURBOMOLE GmbH, since 2007; available from http://www.turbomole.com,” .
  * Riplinger  _et al._ (2013) Christoph Riplinger, Barbara Sandhoefer, Andreas Hansen, and Frank Neese, “Natural triple excitations in local coupled cluster calculations with pair natural orbitals,” J. Chem. Phys. 139, 134101 (2013).
  * Schmitz and Hättig (2016) Gunnar Schmitz and Christof Hättig, “Perturbative triples correction for local pair natural orbital based explicitly correlated ccsd(f12*) using laplace transformation techniques,” J. Chem. Phys. 145, 234107 (2016).
  * Dohm  _et al._ (2020) Sebastian Dohm, Markus Bursch, Andreas Hansen, and Stefan Grimme, “Semiautomated transition state localization for organometallic complexes with semiempirical quantum chemical methods,” J. Chem. Theory Comput. 16, 2002–2012 (2020).
  * Weigend and Ahlrichs (2005) Florian Weigend and Reinhart Ahlrichs, “Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for h to rn: Design and assessment of accuracy,” Phys. Chem. Chem. Phys. 7, 3297–3305 (2005).
  * Weigend (2006) Florian Weigend, “Accurate coulomb-fitting basis sets for h to rn,” Phys. Chem. Chem. Phys. 8, 1057–1065 (2006).
  * Dolg and Cao (2012) Michael Dolg and Xiaoyan Cao, “Relativistic pseudopotentials: Their development and scope of applications,” Chem. Rev. 112, 403–480 (2012).
  * Dunning (1989) Thom H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen,” J. Chem. Phys. 90, 1007–1023 (1989).
  * Weigend  _et al._ (2002) Florian Weigend, Andreas Köhn, and Christof Hättig, “Efficient use of the correlation consistent basis sets in resolution of the identity mp2 calculations,” J. Chem. Phys. 116, 3175–3183 (2002).
  * (32) Marco Häser and Reinhart Ahlrichs, “Improvements on the direct scf method,” J. Comput. Chem. .
  * Tew (2019) David P. Tew, “Principal domains in local correlation theory,” J. Chem. Theory Comput. 15, 6597–6606 (2019).
  * Knizia (2013) Gerald Knizia, “Intrinsic atomic orbitals: An unbiased bridge between quantum theory and chemical concepts,” J. Chem. Theory Comput. 9, 4834–4843 (2013).
  * Masur  _et al._ (2013) Oliver Masur, Denis Usvyat, and Martin Schütz, “Efficient and accurate treatment of weak pairs in local CCSD(T) calculations,” J. Chem. Phys. 139, 164116 (2013).
  * Schütz  _et al._ (2014) Martin Schütz, Oliver Masur, and Denis Usvyat, “Efficient and accurate treatment of weak pairs in local CCSD(T) calculations. II. Beyond the ring approximation,” J. Chem. Phys. 140, 244107 (2014).
  * Häser and Almlöf (1992) Marco Häser and Jan Almlöf, “Laplace transform techniques in møller–plesset perturbation theory,” J. Chem. Phys. 96, 489–494 (1992).
  * Guo  _et al._ (2018) Yang Guo, Christoph Riplinger, Ute Becker, Dimitrios G. Liakos, Yury Minenkov, Luigi Cavallo, and Frank Neese, “Communication: An improved linear scaling perturbative triples correction for the domain based local pair-natural orbital based singles and doubles coupled cluster method [DLPNO-CCSD(T)],” J. Chem. Phys. 148, 011101 (2018).
  * Ma and Werner (2018b) Qianli Ma and Hans-Joachim Werner, “Scalable electron correlation methods. 5. parallel perturbative triples correction for explicitly correlated local coupled cluster with pair natural orbitals,” J. Chem. Theory Comput. 14, 198–215 (2018b).
  * Nagy  _et al._ (2018) Péter R. Nagy, Gyula Samu, and Mihály Kállay, “Optimization of the linear-scaling local natural orbital ccsd(t) method: Improved algorithm and benchmark applications,” Journal of Chemical Theory and Computation 14, 4193–4215 (2018).
  * Tew (2016) David P. Tew, “Explicitly correlated coupled-cluster theory with Brueckner orbitals,” J. Chem. Phys. 145, 074103 (2016).
  * Giner  _et al._ (2018) Emmanuel Giner, David P. Tew, Yann Garniron, and Ali Alavi, “Interplay between electronic correlation and metal–ligand delocalization in the spectroscopy of transition metal compounds: Case study on a series of planar cu2+ complexes,” J. Chem. Theory Comput. 14, 6240–6252 (2018).
  * Li Manni  _et al._ (2019) Giovanni Li Manni, Daniel Kats, David P. Tew, and Ali Alavi, “Role of valence and semicore electron correlation on spin gaps in fe(ii)-porphyrins,” J. Chem. Theory Comput. 15, 1492–1497 (2019).
  * Tew  _et al._ (2011) David P. Tew, Benjamin Helmich, and Christof Hättig, “Local explicitly correlated second-order møller–plesset perturbation theory with pair natural orbitals,” J. Chem. Phys. 135, 074107 (2011).
  * Tew (2021) David P. Tew, “Chapter four - principal domains in f12 explicitly correlated theory,” in  _New Electron Correlation Methods and their Applications, and Use of Atomic Orbitals with Exponential Asymptotes_, Advances in Quantum Chemistry, Vol. 83, edited by Monika Musial and Philip E. Hoggan (Academic Press, 2021) pp. 83–106.
  * Ma  _et al._ (2017) Qianli Ma, Max Schwilk, Christoph Köppl, and Hans-Joachim Werner, “Scalable electron correlation methods. 4. parallel explicitly correlated local coupled cluster with pair natural orbitals (pno-lccsd-f12),” J. Chem. Theory Comput. 13, 4871–4896 (2017).
  * Pavošević  _et al._ (2017) Fabijan Pavošević, Chong Peng, Peter Pinski, Christoph Riplinger, Frank Neese, and Edward F. Valeev, “SparseMaps—A systematic infrastructure for reduced scaling electronic structure methods. V. Linear scaling explicitly correlated coupled-cluster method with pair natural orbitals,” J. Chem. Phys. 146, 174108 (2017).