MIT]Department of Chemistry, Massachusetts Institute of Technology, Cambridge,
MA XPS, XAS, XANES, NEXAFS, RIXS, XES, CBS, CEBE, MP2, CC, DFT, TDDFT, EOM-CC,
ADC, X2C, HF, MOM, DIIS, RHF, SCF, RI-MP2, CBS, MAE, RMSE

# MP2-based composite extrapolation schemes can predict core-ionization
energies for first-row elements with coupled-cluster level accuracy

Anton Morgunov  Henry K. Tran  Oinam Romesh Meitei  Yu-Che Chien  Troy Van
Voorhis<EMAIL_ADDRESS>[

###### Abstract

X-ray photoelectron spectroscopy (XPS) measures core-electron binding energies
(CEBEs) to reveal element-specific insights into chemical environment and
bonding. Accurate theoretical CEBE prediction aids XPS interpretation but
requires proper modeling of orbital relaxation and electron correlation upon
core-ionization. This work systematically investigates basis set selection for
extrapolation to the complete basis set (CBS) limit of CEBEs from $\Delta$MP2
and $\Delta$CC energies across 94 K-edges in diverse organic molecules. We
demonstrate that an alternative composite scheme using $\Delta$MP2 in a large
basis corrected by $\Delta$CC-$\Delta$MP2 difference in a small basis can
quantitatively recover optimally extrapolated $\Delta$CC CEBEs within 0.02 eV.
Unlike $\Delta$CC, MP2 calculations do not suffer from convergence issues and
are computationally cheaper, and, thus, the composite $\Delta$MP2/$\Delta$CC
scheme balances accuracy and cost, overcoming limitations of solely using
either method. We conclude by providing a comprehensive analysis of the choice
of small and large basis sets for the composite schemes and provide practical
recommendations for highly accurate (within 0.10-0.15 eV MAE) ab initio
prediction of XPS spectra.

###### keywords:

X-ray photoelectron spectroscopy, core-electron binding energy, core-
ionization potential, complete basis set, coupled-cluster, maximum overlap
method, second-order Møller–Plesset perturbation theory

##  1 Introduction

Developments in the generation of ultra-short X-ray pulses1, 2, 3 have sparked
a surge of interest in the X-ray Absorption Spectroscopy (XAS), in which
details of atomistic structure are revealed by the ionization or excitation of
core electrons upon x-ray radiation.4 By measuring the kinetic energy of
ionized electrons in an X-ray Photoelectron Spectroscopy (XPS), one can
determine the core-electron binding energy (CEBE). Remarkably, these CEBEs are
not only element-specific but also are affected by the ionized atom’s chemical
(and hence electronic) environment, earning XPS the alternative name of
electron spectroscopy for chemical analysis (ESCA).5 Such specificity enables
XPS to perform surface composition analysis,6, 7, 8, 9 analyze chemical
distribution in quantum dots,10 and infer oxidation states and coordination
numbers of active centers in catalysis.11

The energies at which core electrons are ionized are called edges, with a
preceding K, L, or M specifying the principal quantum number corresponding to
those core orbitals. Excitations of electrons from core to unoccupied orbitals
tend to have similar energies and can be studied through X-ray Absorption
Near-Edge Structure (XANES) or Near-Edge X-ray Absorption Fine Structure
(NEXAFS). NEXAFS focuses on the very near-edge region (within 10 eV) and
reveals the structure and orientation of molecules adsorbed on surfaces.12
While NEXAFS has been traditionally more associated with surface chemistry, it
has also been recently applied to study ultrafast chemical dynamics.13, 14 In
contrast, XANES studies regions farther from the edge (within 10-50 eV) and
allows to infer oxidation states, coordination number, and types of ligands
surrounding the absorbing atom.15 XANES has beeen also applied to chemical
dynamics on attosecond timescales16 and photoinduced proton-coupled electron
transfer.17

X-ray spectroscopy is not limited to XAS, however. An electron in an occupied
orbital of a core-ionized or core-excited particle may transition into the
vacant core orbital, emitting X-ray radiation. Energies of these emitted
photons can be studied through X-ray Emission Spectroscopy (XES), which has
been recently used to enhance understanding of bonding in transition metals.18
A related Resonant Inelastic X-ray Scattering (RIXS) technique measures the
energy and momentum of photons emitted after the core-ionized or core-excited
particle undergoes other low-energy transitions (because the emitted photon
has a different energy, it has been effectively scattered).19 Such scattering
can provide novel insight into mechanisms of photochemical reactions.20

This paper focuses on the ab initio prediction of CEBE, which is often
required for a reliable interpretation of XPS experiments.5 Computational
methods are benchmarked against experimental CEBEs (also known as core-
ionization potentials) aiming to predict experimental values within 0.2 eV.21
Notably, the experimental resolution of XPS instruments determines the
threshold above,22 implying that even greater accuracy will be expected from
theoretical methods as more precise sources of X-ray radiation are developed.

Ejection of a core electron exposes surrounding atoms and electrons to a
positive charge, which results in a significant redistribution of electronic
density resulting in an increased electronic density in the vicinity of an
ionized atom. In order to accurately describe the energetics of this process,
a theoretical model must be able to model the changes in molecular orbitals
(usually referred to as the relaxation effect) and the effects on electron
correlation that this relaxation imposes. Inclusion of relativistic effects
and spin-orbit coupling allows to improve the accuracy of predictions
especially for heavier elements.23 There are two conceptually different
approaches to calculating core-electron binding energy5. In the first family
of methods, response operators modeling the effects of ionization or
excitation are applied to ground-state wave functions. These methods allowed
successful prediction of energies and physical properties of valence-ionized
states,24 without explicit calculation of the wave function of said ionized
states, avoiding the issue of variational collapse to the ground state.
Unsurprisingly, significant effort has been made to investigate the
applications of response methods to the description of core-ionized states.
Time-dependent density functional theory (TDDFT25), a widely used
representative of response methods, with commonly used functionals
systematically underestimates CEBEs by 10 or more eV, an error typically
attributed to the self-interaction and delocalization error characteristic of
DFT methods.26, 27, 28 These errors can be brought down to 1-3 eV if exchange
and correlation functionals are optimized explicitly for core excitations;29
such calibrations, however, are heavily empiric. Similar accuracy without
experimental calibration can be achieved with equation-of-motion coupled-
cluster (EOM-CC) methods30 within a core-valence separation (CVS)
approximation.31 The inclusion of CVS approximation allows one to skip
valence-to-virtual transitions and directly calculate CEBEs, significantly
reducing the computational cost of EOM-CC calculations. CVS-EOM-CCSD can
describe CEBEs with a mean absolute error (MAE) around 1.75 eV, and CVS-EOM-
CCSDT in a quadruple zeta basis set brings down the MAE to 0.15 eV32. A
quantitative agreement with the experiment, with a MAE of 0.07 eV, can be
achieved with CVS-EOM-CCSDTQ32. Sub-electronvolt accuracy can also be achieved
at a lower cost with second and third-order algebraic diagrammatic
construction (ADC) methods within an intermediate state representation
approach.33, 34, 35, 36

Large errors in response theory methods are usually attributed to the
insufficient relaxation of molecular orbitals in the presence of a core hole.
The second family of delta methods aims to properly account for orbital
relaxation by explicitly optimizing the core-ionized state. In such
calculations, the difference between the energy of the core-ionized and ground
state is reported as the CEBE. Remarkably, even at the Hartree-Fock
($\Delta$HF) level of theory, the CEBEs can be calculated within 1 eV, an
error comparable to the accuracy of CVS-EOM-CCSD.37 inclusion of correlation
effects at the level of second-order Møller–Plesset perturbation theory, known
as $\Delta$MP2 energies, can bring the mean absolute error to 0.5 eV.38
Further improvements, which bring the MAE below 0.2 eV, include the use of a
specially calibrated basis set;39 application of the spin-component-scaled
technique;38 or employment of restricted open-shell MP2 theory.40, 41 Recently
developed square gradient minimization algorithm42 is capable of achieving
similar precision by solving restricted open-shell Kohn-Sham equations.26

So far, the highest accuracy method for computing CEBEs is the coupled-cluster
theory,21, 43 often referred to as the golden standard of computational
chemistry. Until recently, however, the application of the $\Delta$CC method
has been sparse primarily because of convergence issues: the creation of a
hole in a core orbital opens room for core-to-virtual and valence-to-core
double transitions that have very similar energies with different signs
(sometimes called a near-degeneracy issue), resulting in the divergence of
coupled-cluster equations. To circumvent these issues, inspired by the CVS
scheme from response theory, Zheng and Cheng in 2019 have shown that exclusion
of such troublesome transitions can solve convergence issues without
significantly affecting the accuracy of results.21 Additional schemes solving
the near-degeneracy issue have been developed and tested by Arias-Martinez et
al. in 2022, who benchmarked the accuracy of $\Delta$CC calculations
extrapolated to the CBS limit by using energies in the aug-cc-pCVnZ
(heavy)/aug-cc-pVDZ (hydrogen) with $n=T,Q$ basis set for 18 organic
molecules. Notably, these extrapolations are only approximations because
calculation of the true CBS limit requires using $\Delta$CC results in
pentuple, and sometimes even hexuple zeta basis sets. Such $\Delta$CC
calculations, however, may be prohibitively expensive even in quadruple-zeta
basis sets for large molecules because of high computational cost of CCSD44
and CCSD(T)45, which scale as $O(N^{6})$ and $O(N^{7})$ respectively, where
$N$ is the number of electrons in the system. A protocol that avoids these
high costs without sacrificing accuracy is thus highly desirable.

In this work, we study the composite wave function based schemes that allow
one to quantitatively (within 0.02 eV) recover the accuracy of $\Delta$CC
calculations by correcting large basis $\Delta$MP2 calculations with a
$\Delta$MP2-$\Delta$CC difference in a small basis. Because MP2 calculations
do not suffer from convergence issues and are computationally less expensive,
such composite schemes make predictions of XPS spectra for large molecules
more feasible. We begin by benchmarking the accuracy (measured as mean
absolute error, MAE) of $\Delta$HF, $\Delta$MP2, $\Delta$CCSD, and
$\Delta$CCSD(T) at calculating CEBEs in basis sets of different sizes. We
continue by analyzing the dependence of the performance of different
extrapolations to the CBS limit of $\Delta$CCSD and $\Delta$CCSD(T)
calculations on the sizes of the basis sets included in the extrapolation and
uncover element-specific trends. We then investigate the impact of the choice
of the large and small basis for the MP2-based composite extrapolation scheme.
We conclude with practical recommendations for highly accurate (within 0.15
eV) ab initio prediction of K-edge CEBEs.

##  2 Theory and computational details

K-edge CEBEs were calculated for carbon (26 molecules, hereafter referred to
as C-series), nitrogen (30 molecules, N-series), oxygen (25 molecules,
O-series), and fluorine (13 molecules, F-series) for a total of 94 data
points. For some molecules, CEBEs were calculated for the ionization of
different atoms within the molecule. Molecules (see Table S1 for a list) were
selected from the table of experimental K-edge CEBEs compiled by Jolly et
al.46 with an attempt to include molecules of varying sizes and in which the
ionized orbitals are in different chemical environments. Whenever multiple
values were reported for one compound, an arithmetic average was used.

Experimental geometries from Computational Chemistry Comparison and Benchmark
DataBase47 were used whenever available. If experimental geometries were
unavailable, the results of full MP2 geometry optimizations with aug-cc-pVQZ
or aug-cc-pVTZ basis sets were taken from the same database. If such results
were not available, an RI-MP248 geometry optimization was performed in cc-
pVQZ49, 50, 51 basis set (Table S1 specifies which geometries were used for
each molecule). The geometry of a neutral species was used to calculate both
ground and core-ionized state energy.

The calculations of core-ionization energies were performed in Dunning basis
sets cc-pVnZ/cc-pCVnZ, with $n=D,T,Q,5$. In this mixed basis set, elements for
which CEBEs are calculated are all treated in core-enhanced basis sets cc-
pCVnZ,52, 53 while the other elements were treated in cc-pVnZ.49, 50 For
example, in the calculation of a CEBE of one of the oxygen atoms in formic
acid, all oxygen atoms were treated in cc-pCVnZ basis set, while carbon and
hydrogen atoms were treated in a cc-pVnZ basis.

The core-electron binding energy (CEBE) in the delta-methods family is
calculated as a difference between the single-point energy of the core-ionized
and ground state (eq. 1).

$E_{\textsc{CEBE}}=E_{\textsc{ion}}-E_{\textsc{ground}}$ (1)

Energies and wave functions of the ground- and ionized states are calculated
using Hartree-Fock (HF). The algorithm for the optimization of the ionized
state is modified to proceed through the maximum overlap method (MOM)54 to
avoid the variational collapse to the ground state. SCF equations for the
ionized state are converged with the direct inversion in the iterative
subspace55 (see Table S2 for specification of DIIS parameters). Resulting HF
wave functions are used as a starting point for calculating correlation
energies with MP2, CCSD, and CCSD(T). The empty core orbital was excluded
(implemented as freezing in PySCF56, 57) from amplitude calculations within
MP2 and CC to improve their convergence. Finally, a 1-electron spin-free X2C
approximation58, 59 has been applied to all calculations to account for scalar
relativistic effects.

If the molecule of interest contains symmetrically equivalent atoms, the core
orbitals in the RHF solution are delocalized over those equivalent atoms, and
a vacation of such orbital leads to the delocalization of the core hole as
well. Practically, this results in an inaccurate orbital relaxation and an
overestimation of core-ionization energies by more than 10 eV. The problem is
resolved by applying a localization scheme to all atoms, e.g., the Boys
localization,60 or by explicitly localizing the core orbitals of symmetrically
identical atoms.

As mentioned previously, precise extrapolation to the CBS limit of coupled-
cluster energies is impractical due to the high computational cost of coupled-
cluster calculations in large basis sets. Hence, lower-cost approximation
schemes must be applied. One commonly used scheme is the two-point
extrapolation from $\Delta$CC energies calculated in triple and quadruple-zeta
basis sets using a two-parameter equation $E=a+bn^{-3}$ described by
Helgaker.61 Energies from this scheme are denoted as X-Y-CCSD and X-Y-CCSD(T),
in which $X$ and $Y$ refer to the size of the basis set based on which the
extrapolation was made.

An alternative approach to approximate CBS values is to use a composite method
(eq. 2-3), which corrects MP2 energies in a large basis set with a difference
between CC and MP2 energies in a small basis set.62 This technique utilizes
the observed linear relationship between MP2 and CC energies, offset with a
constant factor $\delta_{\text{CC-MP2}}$. This scheme has proven to be a good
approximation for energies of non-covalent interactions in large systems.63,
64

$\displaystyle E_{\text{CC}}^{\text{large basis}}$ $\displaystyle\approx
E_{\text{MP2}}^{\text{large basis}}+\delta_{\text{CC-MP2}}^{\text{small
basis}}$ (2) $\displaystyle\delta_{\text{CC-MP2}}^{\text{small basis}}$
$\displaystyle=E_{\text{CC}}^{\text{small basis}}-E_{\text{MP2}}^{\text{small
basis}}$ (3)

We will refer to this scheme as MP2[X Y]+DifZ where X, Y denote the size of
the basis set based on which the large basis MP2 value was obtained, Z refers
to the size of the small basis set, and DifZ or DifZ(T) denotes whether the
difference in the small basis was taken between CCSD and MP2 or CCSD(T) and
MP2 respectively. Calculations of CEBEs were performed in cc-pCVnZ/cc-pVnZ
(described above, $n=D,T,Q,5$), STO-3G,65, 66 STO-6G,65, 66 3,-21G,67, 68
4-31G,69, 70, 71 and 6-31G.68, 69, 72, 73, 74

The latest (2.1.1) version of PySCF56, 57 was used for all CEBE calculations.
Geometry optimizations were performed in the latest (5.0.3) version of
ORCA.75, 76, 77, 78

##  3 Results and Discussion

###  3.1 Convergence of basis for HF, MP2, and CC CEBEs.

The mean absolute errors (MAE) for the CEBEs rely heavily on the basis set’s
size and the method. Remarkably, $\Delta$HF calculations in a double-zeta
basis set produce CEBEs with a MAE of 0.41 eV, significantly lower than the
MAE of $\Delta$MP2 (1.81 eV) or $\Delta$CCSD (1.64 eV) in the same basis (see
Fig. 1). Such a small error could be explained by fortunate cancellation of
systematic errors during the calculation of energy difference between the
ionized and ground states. Surprisingly, the MAE for $\Delta$HF increases with
the basis set size: for triple, quadruple, and pentuple zeta, it is 0.84,
0.91, and 0.92 eV, respectively. In contrast, $\Delta$MP2 and $\Delta$CC
calculations improve as the basis size increases (Table S3). As expected,
$\Delta$CC calculations provide the most accurate CEBEs with a MAE of 0.20 eV
in a quadruple basis set. Inclusion of a perturbative triples correction to
$\Delta$CCSD result increases the MAE, a result in agreement with the work of
Arias-Martinez et al..43

When a similar analysis is performed for CEBEs grouped based on the element on
which the ionized orbitals are located, it is observed that the trends differ
between element series (Fig. 2 and Tables S4-S7). For example, the MAE for
$\Delta$HF increases with the size of the basis set for N-, O-, and F-series
but decreases for C-series. Another peculiarity of the carbon series is that
$\Delta$HF CEBEs have smaller MAE than $\Delta$MP2 CEBEs in triple, quadruple,
and pentuple bases. Finally, $\Delta$CCSD(T) energies (MAE 0.24 eV) are 20%
more accurate than $\Delta$CCSD energies (MAE 0.30 eV) for C-series, but
20-30% less accurate for N- and O-series. The effect of perturbative triples
for the F-series is ambiguous and dependent on the basis size. In effect, the
lack of improvement upon introduction of perturbative triples, described by
Arias-Martinez et al.,43 seems to be an artifact of the selection of molecules
for benchmark studies, which tend to have an abundance of oxygen-based
ionizations.

Figure 1: Accuracy (mean absolute error in eV relative to experimental values)
of $\Delta$HF and $\Delta$post-HF methods in cc-pVXZ/cc-pCVXZ basis sets (X=D,
T, Q, 5) at the calculation of 94 K-edge CEBEs. Error bars show standard
deviations of absolute errors. Figure 2: Accuracy (mean absolute error in eV
relative to experimental values) of $\Delta$HF and $\Delta$post-HF methods in
cc-pVXZ/cc-pCVXZ basis sets (X=D, T, Q, 5) at the calculation of core (K-edge)
ionization energies of electrons localized on carbon (25 molecules), nitrogen
(31 molecules), oxygen (25 molecules), and fluorine (13 molecules). Error bars
show standard deviations of absolute errors.

###  3.2 Choice of the basis sets for $\Delta$CC extrapolations.

Averaged over all molecules, two-point extrapolations T-Q-CCSD and T-Q-CCSD(T)
result in a MAE of 0.18 and 0.17 eV, respectively (Table S8). The inclusion of
double-zeta results worsens the accuracy: D-T-Q-CCSD and D-T-Q-CCSD(T) have a
MAE of 0.26 and 0.20 eV, respectively. However, as seen in Fig. 3 and Tables
S9-S12, these global averages hide the element-specific trends. For example,
the inclusion of double zeta results for carbon series reduces the MAE by a
factor of two from 0.27 to 0.12 eV for T-Q-CCSD vs. D-T-Q-CCSD and from 0.21
to 0.09 eV for T-Q-CCSD(T) vs. D-T-Q-CCSD(T). In contrast, for nitrogen
series, D-T-Q-CCSD is 50% less accurate than T-Q-CCSD (0.21 vs 0.14 eV
respectively). The systematic nature of the increase in accuracy upon
inclusion (for C-series) or exclusion (for all other series) of double zeta
results is confirmed by correlation plots (Fig. S5-6). It should be noted that
such element-specific impact of double zeta basis may be an artifact of the
differences in the parametrization of double zeta basis for different elements
in the Dunning family. Nonetheless, because Dunning basis sets are a
predominant (if not exclusive) choice for $\Delta$CC calculations, such effect
bears practical importance. For all series, the effect of triples correction
is marginal and is more significant when suboptimal (for example, including
double zeta results when it is better not to) extrapolation schemes are used.
Notably, the effect of triples is systematic, here defined as having the same
sign and similar in relative magnitude for all molecules in a series, for
molecules within N-, O-, and F-series (Fig. S9-10), but is irregular for
carbon series: while CEBEs for all but 3 molecules is improved upon inclusion
of perturbative triples, the magnitude of the change is not uniform.

Figure 3: Mean absolute errors (in eV) of D-T-Q-CCSD, D-T-Q-CCSD(T), T-Q-CCSD,
and T-Q-CCSD(T) extrapolations at the calculation of core (K-edge) ionization
energies of electrons localized on carbon (25 molecules), nitrogen (31
molecules), oxygen (25 molecules), and fluorine (13 molecules). Error bars
show standard deviations of absolute errors.

###  3.3 MP2-based composite extrapolation scheme

Extrapolation of $\Delta$MP2 CEBEs to the CBS limit does not result in a
significant increase of accuracy: the MAE for $\Delta$MP2 values extrapolated
from results in $2,3,4,5$-zeta basis sets is 0.26 eV (Tables S13-S16). These
results are improved significantly for carbon and nitrogen series if a
$\delta_{\text{CC-MP2}}^{\text{D}}$ correction is introduced in the double-
zeta polarized basis (D refers to cc-pCVDZ/cc-pVDZ basis): the MAE lowers to
0.17 eV (Table S18). Greater accuracy can be achieved if extrapolation schemes
are chosen differently for each series: the errors can be then reduced to 0.10
eV for C-series (Table S19), 0.14 eV for N-series (Table S20), 0.13 eV for
O-series (Table S21), and 0.15 eV for F-series (Table S22).

Figure 4: Mean absolute errors (in eV) of MP2-based extrapolation schemes in
comparison to CCSD-based extrapolation. Calculations in 2-$\zeta$ are included
only for C-based molecules. MP2 extrapolations do not include results in
pentuple basis. Perturbative triples are included for C- and F-based
molecules. Error bars show standard deviations of absolute errors. Figure 5:
Comparisons of CEBEs for all 94 molecules in the study calculated with
different extrapolations to the CBS limit. Energies are reported as shifts (in
eV) from the lowest value in the element-specific series (set to equal 0.5
eV). The values from MP2-based extrapolation are reported on the x-axis, and
CCSD-based extrapolation is reported on the y-axis. The dotted line is the
curve $y=x$. The 2-$\zeta$ calculations are included only for carbon-based
molecules. MP2 extrapolations do not include results in pentuple basis.

In the previous subsection, we have demonstrated (based on $\Delta$CC CEBEs)
that it is optimal to include double zeta results only in extrapolations for
carbon-based molecules and include perturbative triples correction in carbon-
and fluorine-based molecules. Following these insights, hereafter, we will use
MP2[D T Q]+DifZ(T) scheme for C-series, MP2[T Q]+DifZ for N- and O-series, and
MP2[T Q]+DifZ(T) for F-series (the MP2-based versions of Fig. 3 are Fig.
S2-3). As seen in Fig. 4, such an approach results in a MAE for an MP2-based
composite extrapolation scheme of 0.128 eV, practically identical to the MAE
from CC-based extrapolations of 0.125 eV, when the $\delta_{\text{CC-MP2}}$
correction is calculated in cc-pCVDZ/cc-pVDZ. Notably, if a significantly
smaller regular double zeta basis (3-21G) is used to calculate the
$\delta_{\text{CC-MP2}}$ correction, only a marginally larger MAE of 0.142 eV
is observed.

The CC- and MP2-based extrapolation schemes are similar not only on average
but on a per-molecule basis, as seen in Fig. 5. In fact, MP2[$\infty$]+DifD
recovers $\Delta$CC-based CBS values within 0.04 eV (as measured by RMSE) and
0.03 eV (as measured by MAE), while MP2[$\infty$]+DifD(T) recovers $\Delta$CC
CBS values within 0.03 eV (RMSE) & 0.02 (MAE). Such small differences
establish a quantitative equivalency between the MP2- and CC-based
extrapolation schemes. If the correction is calculated in a regular double-
zeta basis, i.e., 3-21G, an MP2-based scheme recovers CC-based values within
0.10 eV (RMSE) & 0.07 eV (MAE) and 0.09 eV (RMSE) & 0.06 eV (MAE) without and
with perturbative triples respectively. These results are remarkable as
$\Delta$MP2 calculations are both asymptotically and practically faster given
that MP2 scales as $O(N^{5})$ and it only requires a single computation of
perturbative correction, unlike the inherently iterative coupled-cluster
calculations.

####  3.3.1 Effect of the size of the large basis on the accuracy of composite
schemes

Equations 2-3 require the calculation of a $\Delta$MP2 CEBE in a large basis
set. To establish the quantitative equivalence of the MP2-based extrapolation
scheme, we have used MP2 results in double (for C-series only), triple, and
quadruple zeta basis sets. Two questions warrant further investigation: 1. Can
the accuracy be improved even further if MP2 results in the pentuple basis set
are included? 2. Can we use fewer basis sets in MP2 extrapolation? Somewhat
surprisingly, as seen in Fig. 6 and Tables S19-S22, the results in pentuple
zeta basis either do not lower the MAE at all (as for oxygen and fluorine
series) or even slightly increase it. Perhaps even more surprisingly, the
MP2[T Q]+DifD, MP2[T Q 5]+DifD, MP2[Q 5]+DifD, and MP2[5]+DifD all result in
practically identical MAEs for oxygen and fluorine series. The MAEs for
nitrogen series increase slightly as the number of basis sets included in MP2
extrapolation decreases: from 0.14 eV for MP2[T Q]+DifD to 0.18 eV for
MP2[Q]+DifD. In other words, accurate CEBE predictions can be made with less
computationally expensive calculations if extrapolation schemes are chosen
separately for each ionized element.

Figure 6: Mean absolute errors (in eV) of MP2-based extrapolations using
different extrapolations of $E_{\text{MP2}}^{\text{large basis}}$ value for
the calculation of core (K-edge) ionization energies of electrons localized on
carbon (25 molecules), nitrogen (31 molecules), oxygen (25 molecules), and
fluorine (13 molecules). The $\delta_{\text{CC-MP2}}^{\text{D}}$ correction is
calculated with cc-pVDZ/cc-pCVDZ basis set. Error bars show standard
deviations of absolute errors.

####  3.3.2 Effect of the size of the small basis on the accuracy of composite
schemes

We now investigate the question of how small should be the basis set in the
$\delta_{\text{CC-MP2}}^{\text{small basis}}$ correction in the equation 2. As
seen in Fig. 7, the minimal atomic orbital basis sets (such as STO-3G or
STO-6G) are insufficiently large for the use in MP2-based extrapolation
schemes as the MAEs are roughly twice as large as when cc-pCVDZ/cc-pVDZ is
used. In contrast, 3-21G, 4-31G and 6-31G all result in comparable accuracy.

Given the results in Fig. 7 and Tables S18-S22, K-edge CEBEs can be predicted
ab initio within 0.19 eV by simply using the MP2[Q]+DifD scheme, involving an
MP2 calculation in a double zeta and quadruple zeta basis sets and a single
CCSD calculation in the double zeta basis. These results are effectively
identical to the significantly more expensive T-Q-CCSD calculations. To
improve the accuracy of ab initio predictions, one must employ element-
specific extrapolation schemes: MP2[D T Q]+DifD scheme for C-series and MP2[T
Q]+DifD for N-, O-, and F-series. If one is willing to sacrifice up to 10%
accuracy in favor of a lower computational cost, the small basis correction
can be calculated in 3-21G instead.

Figure 7: Mean absolute errors (in eV) of MP2-based extrapolations using
different $\delta_{\text{CC-MP2}}$ corrections for the calculation of (K-edge)
ionization energies of electrons localized on carbon (25 molecules), nitrogen
(31 molecules), oxygen (25 molecules), and fluorine (13 molecules). The
$E_{\text{MP2}}^{\text{large}}$ value is calculated with $2,3,4$-$\zeta$ for
C-series and $3,4$-$\zeta$ for C-, O-, and F-series. Error bars show standard
deviations of absolute errors.

##  4 Conclusions

This paper reports a systematic element-specific study of the impact of basis
set selection for extrapolations to the CBS limit both for $\Delta$MP2 and
$\Delta$CC calculations. Firstly, we reproduce the conclusions of Zheng and
Cheng21 and Arias-Martinez et al.43 that T-Q-CCSD is the best extrapolation
scheme on average, but we note that this result hides element-specific trends:
we found that carbon-based CEBEs are predicted significantly more accurately
if energies from double zeta basis are included in the extrapolation. This
result is important both for experimentalists, especially those who are mostly
concerned with carbon-based XPS, and theoreticians, who select molecules for
benchmark studies, as the performance of a given method on average could
depend not only on the quality of the method itself but also on the relative
proportions of elements that are ionized in the dataset. We also find that
$\Delta$MP2 extrapolated to the CBS limit corrected by $\delta_{\text{CC-
MP2}}^{\text{small basis}}$ quantitatively reproduces optimally extrapolated
$\Delta$CC CEBEs at a fraction of their computational cost. For example, MP2[T
Q]+DifD can be used to estimate CEBEs of molecules for which T-Q-CCSD
calculations are prohibitively expensive. Just as with $\Delta$CC
extrapolations, the choice of the basis for the MP2-based composite method is
element-specific, so the inclusion of a double-zeta basis in the extrapolation
of $\Delta$MP2 results will be useful for carbon-based ionizations. Finally,
this work shows that neither perturbative triples nor MP2 calculations in
pentuple basis set systematically and significantly improve the quality of
predictions. Our results suggest that highly accurate ab initio prediction of
the XPS spectra of large molecules is feasible with currently available
methods.

Future work may be done in three directions. First, a more efficient
implementation of MP2 or CCSD can be used to quantitatively assess the
practical computational costs of proposed composite MP2-based extrapolation
schemes relative to CCSD calculations in triple and quadruple basis sets. An
investigation of the change in accuracy if MP2 results are replaced with a
resolution of identity (RI) approximation of MP2 (RI-MP2) can also be
performed. It would also be beneficial to search for ways to eliminate
dependence on iterative $\Delta$CC results altogether. Second, an
investigation of the accuracy of MP2-based extrapolation schemes for third-row
elements and transition metals is of great practical importance. Third and
finally, all calculations reported in this paper relied upon manual freezing
of core orbitals to avoid the near-degeneracy issue of coupled-cluster
calculations. A black box implementation is especially desirable for L-edge
CEBEs and beyond as the density of orbitals increases with the principal
quantum number.

We thank National Science Foundation for funding this project (CHE-2154938).
Anton Morgunov and Yu-Che Chien are grateful for the additional financial
support by the Undergraduate Research Opportunities Program (UROP) at
Massachusetts Institute of Technology.

The following files are available free of charge.

  * •

Supporting Information: addditional tables and figures for the statistical
metrics for different methods and extrapolation schemes

  * •

https://github.com/anmorgunov/cebe_prediction: code used to run calculations,
analyze the results, and generate the tables and figures in this paper. Also
includes the CEBE values for all molecules in the study.

## References

  * Bergmann et al. 2017 Bergmann, U., Yachandra, V., Yano, J., Eds. _X-Ray Free Electron Lasers_ ; Energy and Environment Series; The Royal Society of Chemistry, 2017; pp P001–463
  * Pellegrini et al. 2016 Pellegrini, C.; Marinelli, A.; Reiche, S. The physics of x-ray free-electron lasers. _Rev. Mod. Phys._ 2016, _88_ , 015006
  * Young et al. 2018 Young, L.; Ueda, K.; Gühr, M.; Bucksbaum, P. H.; Simon, M.; Mukamel, S.; Rohringer, N.; Prince, K. C.; Masciovecchio, C.; Meyer, M. et al. Roadmap of ultrafast X-ray atomic and molecular physics. _Journal of Physics B: Atomic, Molecular and Optical Physics_ 2018, _51_ , 032003
  * Chergui and Collet 2017 Chergui, M.; Collet, E. Photoinduced Structural Dynamics of Molecular Systems Mapped by Time-Resolved X-ray Methods. _Chemical Reviews_ 2017, _117_ , 11025–11065, PMID: 28692268
  * Norman and Dreuw 2018 Norman, P.; Dreuw, A. Simulating X-ray spectroscopies and calculating core-excited states of molecules. _Chemical Reviews_ 2018, _118_ , 7208–7248
  * Fadley 2010 Fadley, C. X-ray photoelectron spectroscopy: Progress and perspectives. _Journal of Electron Spectroscopy and Related Phenomena_ 2010, _178-179_ , 2–32, Trends in X-ray Photoelectron Spectroscopy of solids (theory, techniques and applications)
  * Oswald 2013 Oswald, S. X-ray photoelectron spectroscopy in analysis of surfaces. _Encyclopedia of Analytical Chemistry_ 2013,
  * Matthew 2004 Matthew, J. Surface analysis by Auger and X-ray photoelectron spectroscopy. D. Briggs and J. T. Grant (EDS). Impublications, Chichester, UK and Surfacespectra, Manchester, UK, 2003. 900 pp., ISBN 1-901019-04-7, 900 PP. _Surface and Interface Analysis_ 2004, _36_ , 1647–1647
  * Watts and Wolstenholme 2019 Watts, J. F.; Wolstenholme, J. An introduction to surface analysis by XPS and AES. 2019,
  * Zorn et al. 2011 Zorn, G.; Dave, S. R.; Gao, X.; Castner, D. G. Method for determining the elemental composition and distribution in semiconductor core-shell quantum dots. _Analytical Chemistry_ 2011, _83_ , 866–873
  * Nguyen et al. 2019 Nguyen, L.; Tao, F. F.; Tang, Y.; Dou, J.; Bao, X.-J. Understanding catalyst surfaces during catalysis through near ambient pressure X-ray photoelectron spectroscopy. _Chemical Reviews_ 2019, _119_ , 6822–6905
  * Hähner 2006 Hähner, G. Near edge X-ray absorption fine structure spectroscopy as a tool to probe electronic and structural properties of thin organic films and liquids. _Chem. Soc. Rev._ 2006, _35_ , 1244–1255
  * Segatta et al. 2020 Segatta, F.; Nenov, A.; Orlandi, S.; Arcioni, A.; Mukamel, S.; Garavelli, M. Exploring the capabilities of optical pump X-ray probe NEXAFS spectroscopy to track photo-induced dynamics mediated by conical intersections. _Faraday Discuss._ 2020, _221_ , 245–264
  * Attar et al. 2017 Attar, A. R.; Bhattacherjee, A.; Pemmaraju, C. D.; Schnorr, K.; Closser, K. D.; Prendergast, D.; Leone, S. R. Femtosecond x-ray spectroscopy of an electrocyclic ring-opening reaction. _Science_ 2017, _356_ , 54–59
  * Iglesias-Juez et al. 2022 Iglesias-Juez, A.; Chiarello, G. L.; Patience, G. S.; Guerrero-Pérez, M. O. Experimental methods in chemical engineering: X-ray absorption spectroscopy—XAS, XANES, EXAFS. _The Canadian Journal of Chemical Engineering_ 2022, _100_ , 3–22
  * Kraus et al. 2018 Kraus, P. M.; Zürch, M.; Cushing, S. K.; Neumark, D. M.; Leone, S. R. The ultrafast X-ray spectroscopic revolution in chemical dynamics. _Nature Reviews Chemistry_ 2018, _2_ , 82–94
  * Soley et al. 2022 Soley, M. B.; Videla, P. E.; Nibbering, E. T. J.; Batista, V. S. Ultrafast Charge Relocation Dynamics in Enol–Keto Tautomerization Monitored with a Local Soft-X-ray Probe. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 8254–8263, PMID: 36018775
  * Britz et al. 2019 Britz, A.; Gawelda, W.; Assefa, T. A.; Jamula, L. L.; Yarranton, J. T.; Galler, A.; Khakhulin, D.; Diez, M.; Harder, M.; Doumy, G. et al. Using Ultrafast X-ray Spectroscopy To Address Questions in Ligand-Field Theory: The Excited State Spin and Structure of [Fe(dcpp)2]2+. _Inorganic Chemistry_ 2019, _58_ , 9341–9350, PMID: 31241335
  * Ament et al. 2011 Ament, L. J. P.; van Veenendaal, M.; Devereaux, T. P.; Hill, J. P.; van den Brink, J. Resonant inelastic x-ray scattering studies of elementary excitations. _Rev. Mod. Phys._ 2011, _83_ , 705–767
  * Wernet et al. 2015 Wernet, P.; Kunnus, K.; Josefsson, I.; Rajkovic, I.; Quevedo, W.; Beye, M.; Schreck, S.; Grübel, S.; Scholz, M.; Nordlund, D. et al. Orbital-specific mapping of the ligand exchange dynamics of fe(co)5 in solution. _Nature_ 2015, _520_ , 78–81
  * Zheng and Cheng 2019 Zheng, X.; Cheng, L. Performance of delta-coupled-cluster methods for calculations of core-ionization energies of first-row elements. _Journal of Chemical Theory and Computation_ 2019, _15_ , 4945–4955
  * Greczynski and Hultman 2020 Greczynski, G.; Hultman, L. X-ray photoelectron spectroscopy: Towards reliable binding energy referencing. _Progress in Materials Science_ 2020, _107_ , 100591
  * Kotsis and Staemmler 2006 Kotsis, K.; Staemmler, V. Ab initio calculations of the O1s XPS spectra of ZnO and Zn oxo compounds. _Phys. Chem. Chem. Phys._ 2006, _8_ , 1490–1498
  * Dreuw and Head-Gordon 2005 Dreuw, A.; Head-Gordon, M. Single-reference ab initio methods for the calculation of excited states of large molecules. _Chemical Reviews_ 2005, _105_ , 4009–4037
  * Marques and Gross 2004 Marques, M.; Gross, E. Time-dependent density functional theory. _Annual Review of Physical Chemistry_ 2004, _55_ , 427–455
  * Hait and Head-Gordon 2020 Hait, D.; Head-Gordon, M. Highly accurate prediction of core spectra of molecules at density functional theory cost: Attaining sub-electronvolt error from a restricted open-shell kohn–sham approach. _The Journal of Physical Chemistry Letters_ 2020, _11_ , 775–786
  * Hait and Head-Gordon 2018 Hait, D.; Head-Gordon, M. Delocalization errors in density functional theory are essentially quadratic in fractional occupation number. _The Journal of Physical Chemistry Letters_ 2018, _9_ , 6280–6288
  * Mori-Sánchez et al. 2006 Mori-Sánchez, P.; Cohen, A. J.; Yang, W. Many-electron self-interaction error in approximate density functionals. _The Journal of Chemical Physics_ 2006, _125_ , 201102
  * Lestrange et al. 2015 Lestrange, P. J.; Nguyen, P. D.; Li, X. Calibration of energy-specific TDDFT for modeling K-edge xas spectra of light elements. _Journal of Chemical Theory and Computation_ 2015, _11_ , 2994–2999
  * Nooijen and Bartlett 1995 Nooijen, M.; Bartlett, R. J. Description of core-excitation spectra by the open-shell electron-attachment equation-of-motion coupled cluster method. _The Journal of Chemical Physics_ 1995, _102_ , 6735–6756
  * Cederbaum et al. 1980 Cederbaum, L. S.; Domcke, W.; Schirmer, J. Many-body theory of core holes. _Physical Review A_ 1980, _22_ , 206–222
  * Liu et al. 2019 Liu, J.; Matthews, D.; Coriani, S.; Cheng, L. Benchmark calculations of K-edge ionization energies for first-row elements using scalar-relativistic core–valence-separated equation-of-motion coupled-cluster methods. _Journal of Chemical Theory and Computation_ 2019, _15_ , 1642–1651
  * Schirmer 1982 Schirmer, J. Beyond the random-phase approximation: A new approximation scheme for the polarization propagator. _Physical Review A_ 1982, _26_ , 2395–2416
  * Schirmer and Trofimov 2004 Schirmer, J.; Trofimov, A. B. Intermediate state representation approach to physical properties of electronically excited molecules. _The Journal of Chemical Physics_ 2004, _120_ , 11449–11464
  * Wenzel and Dreuw 2016 Wenzel, J.; Dreuw, A. Physical properties, exciton analysis, and visualization of core-excited states: An intermediate state representation approach. _Journal of Chemical Theory and Computation_ 2016, _12_ , 1314–1330
  * Wenzel et al. 2014 Wenzel, J.; Wormit, M.; Dreuw, A. Calculating core-level excitations and x-ray absorption spectra of medium-sized closed-shell molecules with the algebraic-diagrammatic construction scheme for the polarization propagator. _Journal of Computational Chemistry_ 2014, _35_ , 1900–1915
  * Besley et al. 2009 Besley, N. A.; Gilbert, A. T.; Gill, P. M. Self-consistent-field calculations of core excited states. _The Journal of Chemical Physics_ 2009, _130_ , 124308
  * Smiga and Grabowski 2018 Smiga, S.; Grabowski, I. Spin-component-scaled delta MP2 parametrization: Toward a simple and reliable method for ionization energies. _Journal of Chemical Theory and Computation_ 2018, _14_ , 4780–4790
  * Shim et al. 2011 Shim, J.; Klobukowski, M.; Barysz, M.; Leszczynski, J. Calibration and applications of the $\Delta$MP2 method for calculating core electron binding energies. _Phys. Chem. Chem. Phys._ 2011, _13_ , 5703–5711
  * Ye and Van Voorhis 2020 Ye, H.-Z.; Van Voorhis, T. Self-consistent Møller-Plesset Perturbation Theory For Excited States. 2020; https://arxiv.org/abs/2008.10777
  * Garner and Neuscamman 2020 Garner, S. M.; Neuscamman, E. Core excitations with excited state mean field and perturbation theory. _The Journal of Chemical Physics_ 2020, _153_ , 154102
  * Hait and Head-Gordon 2020 Hait, D.; Head-Gordon, M. Excited state orbital optimization via minimizing the square of the gradient: General approach and application to singly and doubly excited states via density functional theory. _Journal of Chemical Theory and Computation_ 2020, _16_ , 1699–1710
  * Arias-Martinez et al. 2022 Arias-Martinez, J. E.; Cunha, L. A.; Oosterbaan, K. J.; Lee, J.; Head-Gordon, M. Accurate core excitation and ionization energies from a state-specific coupled-cluster singles and doubles approach. _Physical Chemistry Chemical Physics_ 2022, _24_ , 20728–20741
  * Purvis and Bartlett 1982 Purvis, G. D.; Bartlett, R. J. A full coupled-cluster singles and doubles model: The inclusion of disconnected triples. _The Journal of Chemical Physics_ 1982, _76_ , 1910–1918
  * Raghavachari et al. 1989 Raghavachari, K.; Trucks, G. W.; Pople, J. A.; Head-Gordon, M. A fifth-order perturbation comparison of electron correlation theories. _Chemical Physics Letters_ 1989, _157_ , 479–483
  * Jolly et al. 1984 Jolly, W.; Bomben, K.; Eyermann, C. Core-electron binding energies for gaseous atoms and molecules. _Atomic Data and Nuclear Data Tables_ 1984, _31_ , 433–493
  * ccc 2022 NIST Computational Chemistry Comparison and Benchmark Database. 2022; https://cccbdb.nist.gov/
  * Stoychev et al. 2018 Stoychev, G. L.; Auer, A. A.; Neese, F. Efficient and accurate prediction of nuclear magnetic resonance shielding tensors with double-hybrid density functional theory. _Journal of Chemical Theory and Computation_ 2018, _14_ , 4756–4771
  * Dunning 1988 Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. i. the atoms boron through neon and hydrogen. _The Journal of Chemical Physics_ 1988, _90_ , 1007–1023
  * Woon and Dunning 1993 Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. III. the atoms aluminum through argon. _The Journal of Chemical Physics_ 1993, _98_ , 1358–1371
  * Weigend et al. 2002 Weigend, F.; Köhn, A.; Hättig, C. Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations. _The Journal of Chemical Physics_ 2002, _116_ , 3175–3183
  * Woon and Dunning 1995 Woon, D. E.; Dunning, T. H. Gaussian basis sets for use in correlated molecular calculations. V. Core-valence basis sets for boron through neon. _The Journal of Chemical Physics_ 1995, _103_ , 4572–4585
  * Peterson and Dunning 2002 Peterson, K. A.; Dunning, T. H. Accurate correlation consistent basis sets for molecular core–valence correlation effects: The second row atoms al–ar, and the first row atoms B–ne revisited. _The Journal of Chemical Physics_ 2002, _117_ , 10548–10560
  * Gilbert et al. 2008 Gilbert, A. T.; Besley, N. A.; Gill, P. M. Self-consistent field calculations of excited states using the maximum overlap method (MOM). _The Journal of Physical Chemistry A_ 2008, _112_ , 13164–13171
  * Pulay 1980 Pulay, P. Convergence acceleration of iterative sequences. the case of scf iteration. _Chemical Physics Letters_ 1980, _73_ , 393–398
  * 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. et al. PySCF: The python-based simulations of Chemistry Framework. _WIREs Computational Molecular Science_ 2017, _8_
  * Sun et al. 2020 Sun, Q.; Zhang, X.; Banerjee, S.; Bao, P.; Barbry, M.; Blunt, N. S.; Bogdanov, N. A.; Booth, G. H.; Chen, J.; Cui, Z.-H. et al. Recent developments in the PySCF program package. _The Journal of Chemical Physics_ 2020, _153_ , 024109
  * Dyall 2001 Dyall, K. G. Interfacing relativistic and nonrelativistic methods. iv. one- and two-electron scalar approximations. _The Journal of Chemical Physics_ 2001, _115_ , 9136–9143
  * Liu and Peng 2009 Liu, W.; Peng, D. Exact two-component Hamiltonians revisited. _The Journal of Chemical Physics_ 2009, _131_ , 031104
  * Boys 1960 Boys, S. F. Construction of some molecular orbitals to be approximately invariant for changes from one molecule to another. _Reviews of Modern Physics_ 1960, _32_ , 296–299
  * Helgaker et al. 1997 Helgaker, T.; Klopper, W.; Koch, H.; Noga, J. Basis-set convergence of correlated calculations on water. _The Journal of Chemical Physics_ 1997, _106_ , 9639–9646
  * Marshall et al. 2011 Marshall, M. S.; Burns, L. A.; Sherrill, C. D. Basis set convergence of the coupled-cluster correction, mp2ccsd(t): Best practices for benchmarking non-covalent interactions and the attendant revision of the S22, NBC10, HBC6, and HSG databases. _The Journal of Chemical Physics_ 2011, _135_ , 194102
  * Sinnokrot and Sherrill 2006 Sinnokrot, M. O.; Sherrill, C. D. High-accuracy quantum mechanical studies of $\pi$-$\pi$ interactions in benzene dimers. _The Journal of Physical Chemistry A_ 2006, _110_ , 10656–10668
  * Pitoňák et al. 2009 Pitoňák, M.; Janowski, T.; Neogrády, P.; Pulay, P.; Hobza, P. Convergence of the ccsd(t) correction term for the stacked complex methyl adenine-methyl thymine: comparison with lower-cost alternatives. _Journal of Chemical Theory and Computation_ 2009, _5_ , 1761–1766
  * Hehre et al. 1969 Hehre, W. J.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. I. Use of Gaussian Expansions of Slater-Type Atomic Orbitals. _J. Chem. Phys._ 1969, _51_ , 2657–2664
  * Hehre et al. 1970 Hehre, W. J.; Ditchfield, R.; Stewart, R. F.; Pople, J. A. Self-Consistent Molecular Orbital Methods. IV. Use of Gaussian Expansions of Slater-Type Orbitals. Extension to Second-Row Molecules. _J. Chem. Phys._ 1970, _52_ , 2769–2773
  * Binkley et al. 1980 Binkley, J. S.; Pople, J. A.; Hehre, W. J. Self-consistent molecular orbital methods. 21. Small split-valence basis sets for first-row elements. _J. Am. Chem. Soc._ 1980, _102_ , 939–947
  * Gordon et al. 1982 Gordon, M. S.; Binkley, J. S.; Pople, J. A.; Pietro, W. J.; Hehre, W. J. Self-consistent molecular-orbital methods. 22. Small split-valence basis sets for second-row elements. _J. Am. Chem. Soc._ 1982, _104_ , 2797–2803
  * Ditchfield et al. 1971 Ditchfield, R.; Hehre, W. J.; Pople, J. A. Self-Consistent Molecular-Orbital Methods. IX. An Extended Gaussian-Type Basis for Molecular-Orbital Studies of Organic Molecules. _J. Chem. Phys._ 1971, _54_ , 724–728
  * Hehre and Lathan 1972 Hehre, W. J.; Lathan, W. A. Self-Consistent Molecular Orbital Methods. XIV. An Extended Gaussian-Type Basis for Molecular Orbital Studies of Organic Molecules. Inclusion of Second Row Elements. _J. Chem. Phys._ 1972, _56_ , 5255–5257
  * Hehre and Pople 1972 Hehre, W. J.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XIII. An Extended Gaussian-Type Basis for Boron. _J. Chem. Phys._ 1972, _56_ , 4233–4234
  * Dill and Pople 1975 Dill, J. D.; Pople, J. A. Self-consistent molecular orbital methods. XV. Extended Gaussian-type basis sets for lithium, beryllium, and boron. _J. Chem. Phys._ 1975, _62_ , 2921–2923
  * Francl et al. 1982 Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; DeFrees, D. J.; Pople, J. A. Self-consistent molecular orbital methods. XXIII. A polarization-type basis set for second-row elements. _J. Chem. Phys._ 1982, _77_ , 3654–3665
  * Hehre et al. 1972 Hehre, W. J.; Ditchfield, R.; Pople, J. A. Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules. _J. Chem. Phys._ 1972, _56_ , 2257–2261
  * Neese 2011 Neese, F. The Orca Program System. _WIREs Computational Molecular Science_ 2011, _2_ , 73–78
  * Neese et al. 2020 Neese, F.; Wennmohs, F.; Becker, U.; Riplinger, C. The Orca Quantum Chemistry Program Package. _The Journal of Chemical Physics_ 2020, _152_ , 224108
  * Valeev 2022 Valeev, E. F. Libint: A library for the evaluation of molecular integrals of many-body operators over Gaussian functions. http://libint.valeyev.net/, 2022; version 2.8.0
  * Lehtola et al. 2018 Lehtola, S.; Steigemann, C.; Oliveira, M. J.; Marques, M. A. Recent developments in libxc — a comprehensive library of functionals for density functional theory. _SoftwareX_ 2018, _7_ , 1–5