#  Efficient implementation of the single-reference algebraic diagrammatic
construction theory for charged excitations: Applications to the TEMPO radical
and DNA base pairs

Samragni Banerjee  Department of Chemistry and Biochemistry, The Ohio State
University, Columbus, Ohio 43210, USA  Alexander Yu. Sokolov
<EMAIL_ADDRESS>Department of Chemistry and Biochemistry, The Ohio State
University, Columbus, Ohio 43210, USA

###### Abstract

We present an efficient implementation of the second- and third-order single-
reference algebraic diagrammatic construction theory for electron attachment
(EA) and ionization (IP) energies and spectra (EA/IP-ADC($n$), $n=2,3$). Our
new EA/IP-ADC program features spin adaptation for closed-shell systems,
density fitting for efficient handling of the two-electron integral tensors,
as well as vectorized and parallel implementation of tensor contractions. We
demonstrate capabilities of our efficient implementation by applying the
EA/IP-ADC($n$) (n = 2, 3) methods to compute the photoelectron spectrum of the
TEMPO radical, as well as the vertical and adiabatic electron affinities of
TEMPO and two DNA base pairs (guanine–cytosine and adenine–thymine). The
spectra and electron affinities computed using large diffuse basis sets with
up to 1028 molecular orbitals are found to be in a good agreement with the
best available results from the experiment and theoretical simulations.

##  I Introduction

Accurate theoretical simulations of charged electronic states are crucial for
reliable predictions of many important properties of molecules and materials,
such as electron affinities (EA), ionization potentials (IP), band gaps, and
photoelectron spectra.Cederbaum:1975p290 ; VonNiessen:1984p57 ; Ortiz:2012p123
; Hirata:2015p1595  Computations of charged excitations face many challenges
associated with the description of charge distribution and open-shell
electronic states that require accurate treatment of orbital relaxation and
electron correlation effects. Reliable simulations of electron attachment and
ionization can be performed using high-level ab initio methods, such as state-
specific or equation-of-motion (EOM) coupled cluster theory
(CC)Nooijen:1992p55 ; Nooijen:1993p15 ; Nooijen:1995p1681 ; Stanton:1993p7029
; Krylov:2008 ; Kowalski:2014p094102 ; BhaskaranNair:2016p144101 ;
Peng:2018p4335  and multireference configuration interaction,Fleig:2008p270 ;
Streit:2011p22  but these methods are usually limited to small systems due to
their high computational cost. On the other hand, computationally efficient
approaches, such as those based on time-dependent many-body perturbation
theory, often combined with density functional theory,Mckechnie:2015p194114 ;
Nguyen:2012p081101 ; Hedin:1965p796 ; Faleev:2004p126406 ;
vanSchilfgaarde:2006p226402 ; Neaton:2006p216405 ; Samsonidze:2011p186404 ;
vanSetten:2013p232 ; Reining:2017pe1344  can be applied to larger systems but
may be unreliable and are difficult to improve systematically. The development
of accurate yet efficient computational methods for charged electronic states
is an ongoing challenge and is an active area of
research.Pavosevic:2017p121101 ; Lange:2018p4224 ; Backhouse:2020p6294 ;
Dempwolff:2019p064108 ; Dempwolff:2020p024113 ; Dempwolff:2020p024125

An attractive alternative to traditional coupled cluster theory for simulating
charged electronic states is algebraic diagrammatic construction (ADC)
theorySchirmer:1982p2395 ; Schirmer:1983p1237 ; Schirmer:1991p4647 ;
Mertins:1996p2140 ; Schirmer:2004p11449 ; Dreuw:2014p82  that offers an
efficient approach for calculating excitation energies and transition
intensities from poles and residues of a many-body propagator approximated in
a time-independent perturbation series. Since its formulation in 1982, single-
reference ADC theory has been applied to a wide range of problems in chemistry
and spectroscopy, including simulating UV/VisSchirmer:1982p2395 ;
Starcke:2009p024104 ; Harbach:2014p064113  and X-ray absorption,Barth:1985p867
; Wenzel:2014p1900  photoelectron,Schirmer:1983p1237 ; Schirmer:1998p4734 ;
Trofimov:2005p144115 ; Angonoa:1987p6789 ; Schirmer:2001p10621 ;
Thiel:2003p2088 ; Dempwolff:2019p064108 ; Dempwolff:2020p024113 ;
Dempwolff:2020p024125  and two-photon absorption
spectra.Knippenberg:2012p064107  The ADC methods for charged excitations
(EA/IP), first proposed in 1983, were based on the Dyson framework that
involves perturbative self-energy expansion of the Dyson equation and couples
the EA and IP components of the one-electron propagator.Schirmer:1983p1237
Later work by Schirmer, Trofimov, and co-workersSchirmer:1998p4734 ;
Trofimov:2005p144115  established the non-Dyson EA/IP-ADC formalism, where the
EA and IP components of the propagator are decoupled. Recently, our group
developed a multireference formulation of the non-Dyson EA/IP-ADC theory for
simulations of charged excitations in strongly correlated chemical
systems.Sokolov:2018p204113 ; Chatterjee:2019p5908 ; Chatterjee:2020p6343

Several implementations of the single-reference non-Dyson EA/IP-ADC methods
were reported in the literature, most of them limited to calculations of only
ionization processes (IP). In 2005, Trofimov and Schirmer presented the first
pilot implementation of the non-Dyson IP-ADC($n$) method up to third order in
perturbation theory. More recently, two efficient implementations of non-Dyson
IP-ADC($n$) ($n\leq 3$) were reported by Schneider et al.Schneider:2015p144103
and Dempwolff et al.Dempwolff:2019p064108 ; Dempwolff:2020p024113 ;
Dempwolff:2020p024125  and the latter is available in the Q-Chem
package.qchem:44  Implementations of the non-Dyson EA-ADC($n$) ($n\leq 3$)
methods were limited to two early benchmark studies,Trofimov:2011p77 ;
Schneider:2015  a more extensive benchmark of EA-ADC($n$) ($n\leq 3$) for
closed- and open-shell systems performed in our group,Banerjee:2019p224112
and a recent efficient implementation of EA-ADC(2) by Liu et
al.Liu:2020p174109  To the best of our knowledge, an efficient implementation
of EA-ADC(3) has not been reported. While applications of non-Dyson IP-ADC(3)
were presented for systems with up to 639 molecular
orbitals,Dempwolff:2019p064108 ; Dempwolff:2020p024113  the largest non-Dyson
EA-ADC(3) calculations reported to date was limited to 298
orbitals.Banerjee:2019p224112

In this work, we present an efficient implementation of the non-Dyson EA- and
IP-ADC($n$) methods ($n\leq 3$), now available as a module in the PySCF
software package.Sun:2020p024109  Our EA/IP-ADC($n$) program features the
spin-restricted (RADC) and unrestricted (UADC) codes for closed- and open-
shell systems, respectively, takes advantage of the vectorized and parallel
implementation of tensor contractions, and employs the density fitting
approximationWhitten:1973p4496 ; Vahtras:1993p514 ; Feyereisen:1993p359 ;
Rendell:1994p400 ; Weigend:2002p4285  for efficient disk and memory
management. We demonstrate the efficiency of our implementation by performing
the EA- and IP-UADC(3) calculations for the open-shell
(2,2,6,6-tetramethylpiperidin-1-yl)oxyl (TEMPO) radical with up to 758
orbitals and EA-RADC(3) calculations for the closed-shell DNA base pairs with
up to 1028 basis functions, which are the largest non-Dyson EA/IP-ADC(3)
calculations reported to date.

We begin our discussion with a brief overview of ADC (Section II), followed by
the overview of our implementation (Section III), where we discuss the general
algorithm (Section III.1) along with the details of spin adaptation (Section
III.2) and the density fitting approximation (Section III.3). We next present
computational details in Section IV and benchmark the accuracy of the density
fitting approximation along with the efficiency of our implementation in
Section V.1. Finally, we present our results for the TEMPO radical in Section
V.2 and two DNA base pairs (guanine–cytosine, adenine–thymine) in Section V.3.
We present our conclusions in Section VI.

##  II Theory

###  II.1 Non-Dyson ADC theory for the one-particle propagator

We first give a brief overview of the non-Dyson algebraic diagrammatic
construction (ADC) theorySchirmer:1982p2395 ; Schirmer:1983p1237 ;
Schirmer:1991p4647 ; Mertins:1996p2140 ; Schirmer:2004p11449 ; Dreuw:2014p82 ;
Banerjee:2019p224112  for electron attachment and removal, where the central
object of interest is the retarded single-particle Green’s function

$\displaystyle G_{pq}(\omega)$
$\displaystyle=G_{pq}^{+}(\omega)+G_{pq}^{-}(\omega)$
$\displaystyle=\langle\Psi_{0}^{N}|a_{p}(\omega-H+E_{0}^{N})^{-1}a_{q}^{\dagger}|\Psi_{0}^{N}\rangle$
$\displaystyle+\langle\Psi_{0}^{N}|a_{q}^{\dagger}(\omega+H-E_{0}^{N})^{-1}a_{p}|\Psi_{0}^{N}\rangle$
(1)

Here, the forward ($G_{pq}^{+}(\omega)$) and backward ($G_{pq}^{-}(\omega)$)
components of the propagator describe electron attachment (EA) and ionization
(IP) processes, respectively, $|\Psi_{0}^{N}\big{\rangle}$ and E${}_{0}^{N}$
are the eigenfunction and eigenvalue of the electronic Hamiltonian H of
$N$-electron system, and a frequency $\omega\equiv\omega^{\prime}+i\eta$ is
expressed in terms of its real ($\omega^{\prime}$) and infinitesimal imaginary
($i\eta$) components. The operators $a_{p}^{{\dagger}}$ and $a_{p}$ are the
usual fermionic creation and annihilation operators, respectively.

The forward and backward propagators in Eq. 1 can be expressed in the matrix
form

$\displaystyle\textbf{G}_{\pm}(\omega)=\mathbf{\tilde{X}}_{\pm}(\omega\textbf{1}-\boldsymbol{\tilde{\Omega}}_{\pm})^{-1}\mathbf{\tilde{X}}_{\pm}^{\dagger}$
(2)

where $\boldsymbol{\tilde{\Omega}}_{\pm}$ are the diagonal matrices of the
exact vertical attachment ($\omega_{n}=E_{n}^{N+1}-E_{0}^{N}$) and ionization
($\omega_{n}=E_{0}^{N}-E_{n}^{N-1}$) energies and $\mathbf{\tilde{X}}_{\pm}$
are the matrices of the spectroscopic amplitudes with elements
$\tilde{X}_{+pn}=\langle\Psi_{0}^{N}|a_{p}|\Psi_{n}^{N+1}\rangle$ and
$\tilde{X}_{-qn}=\langle\Psi_{0}^{N}|a_{q}^{\dagger}|\Psi_{n}^{N-1}\rangle$,
where $\ket{\Psi_{n}^{N+1}}$ and $\ket{\Psi_{n}^{N-1}}$ are the exact
eigenstates and $E_{n}^{N+1}$ and $E_{n}^{N-1}$ are the exact eigenvalues of
the Hamiltonian $H$ for the $(N+1)$\- and $(N-1)$-electron system,
respectively. Eq. 2 is known as the spectral (or Lehmann) representation of
the propagator.Kobe:1962p448

The propagator in Eqs. 1 and 2 expressed in the basis of the exact eigenstates
$\ket{\Psi_{n}^{N+1}}$ and $\ket{\Psi_{n}^{N-1}}$ is computationally
expensive. In non-Dyson ADC, the propagator is rewritten in an approximate (in
general, non-orthogonal) basis of the ($N\pm 1$)-electron states

$\displaystyle\textbf{G}_{\pm}(\omega)=\textbf{T}_{\pm}(\omega\textbf{S}_{\pm}-\textbf{M}_{\pm})^{-1}\textbf{T}_{\pm}^{\dagger}$
(3)

where the $\textbf{M}_{\pm}$ matrices are no longer diagonal, in contrast to
$\boldsymbol{\tilde{\Omega}}_{\pm}$ in Eq. 2. In Eq. 3, $\textbf{M}_{\pm}$ and
$\textbf{T}_{\pm}$, known as the effective Hamiltonian and transition moment
matrices, contain information about charged excitation energies and
probabilities of electronic transitions, respectively. Since the approximate
basis states are in general non-orthogonal, information about their overlap is
contained in the $\textbf{S}_{\pm}$ matrices. For all single-reference ADC
methods discussed in this work, the approximate basis is orthonormal and,
thus, we will assume $\textbf{S}_{\pm}=\textbf{1}$ henceforth.

Starting with Eq. 3, the approximate propagator $\textbf{G}_{\pm}(\omega)$ is
expanded in a perturbative series

G${}_{\pm}(\omega)$
$\displaystyle=\textbf{G${}_{\pm}^{(0)}(\omega)$}+\textbf{G${}_{\pm}^{(1)}(\omega)$}+...+\textbf{G${}_{\pm}^{(n)}(\omega)$}+...$
(4)

Truncating this expansion at the $n$-th order defines the propagator of the
$n$-th order non-Dyson ADC approximation for EA or IP (EA/IP-ADC($n$)).
Working equations for G${}_{\pm}(\omega)$ are derived by evaluating
$\textbf{M}_{\pm}$ and $\textbf{T}_{\pm}$ up to the $n$-th order in
perturbation theory using the intermediate state representation
approachSchirmer:1991p4647 ; Mertins:1996p2140 ; Schirmer:2004p11449  or via
the formalism of the effective Liouvillean theory.Mukherjee:1989p257 ;
Sokolov:2018p204113 ; Banerjee:2019p224112

In the following, we briefly overview the derivation of the non-Dyson EA/IP-
ADC($n$) equations using effective Liouvillean theory. In this approach, the
exact ground-state wavefunction is expressed via a unitary transformation

$\displaystyle|\Psi_{0}^{N}\rangle=e^{T-T^{\dagger}}|\Phi\rangle$ (5)

where $|\Phi\rangle$ is a reference Slater determinant and the operator $T$
generates all possible excitations from the occupied to virtual (external)
orbitals labeled using the $i,j,k,l,\ldots$ and $a,b,c,d,\ldots$ indices,
respectively

$\displaystyle T=\sum_{m}^{N}T_{m},\quad
T_{m}=\frac{1}{(m!)^{2}}\sum_{ijab\ldots}t_{ij\ldots}^{ab\ldots}a_{a}^{\dagger}a_{b}^{\dagger}\ldots
a_{j}a_{i}$ (6)

To define perturbative series (4), the electronic Hamiltonian $H$ is
partitioned into the zeroth-order part $H^{(0)}$ and a perturbation
$V=H-H^{(0)}$. This allows to determine contributions to $\textbf{M}_{\pm}$
and $\textbf{T}_{\pm}$ at each order in perturbation theory.

The $n$th-order contributions to the EA-ADC matrices have the form:

$\displaystyle M_{+\mu\nu}^{(n)}$
$\displaystyle=\sum_{klm}^{k+l+m=n}\langle\Phi|[h_{+\mu}^{(k)},[\tilde{H}^{(l)},h_{+\nu}^{(m)\dagger}]]_{+}|\Phi\rangle$
(7) $\displaystyle T_{+p\nu}^{(n)}$
$\displaystyle=\sum_{kl}^{k+l=n}\langle\Phi|[\tilde{a}_{p}^{(k)},h_{+\nu}^{(l)\dagger}]_{+}|\Phi\rangle$
(8)

where operators $\tilde{H}^{(k)}$ and $\tilde{a}_{p}^{(k)}$ are the $k$-th-
order contributions to the effective Hamiltonian
$\tilde{{H}}=e^{-(T-T^{\dagger})}He^{(T-T^{\dagger})}$ and observable
$\tilde{a}_{p}=e^{-(T-T^{\dagger})}a_{p}e^{(T-T^{\dagger})}$ operators, while
$[\ldots]$ and $[\ldots]_{+}$ denote commutator and anticommutator,
respectively. The operators $h_{+\mu}^{(k)\dagger}$ in Eqs. 7 and 8 are used
to construct a set of basis states
$\ket{\Psi^{(k)}_{+\mu}}=h_{+\mu}^{(k)\dagger}\ket{\Phi}$ necessary for
expanding the eigenstates of the $(N+1)$-electron system in the EA-ADC
equations. Similarly, the $n$th-order contributions to the IP-ADC matrices are
given as:

$\displaystyle M_{-\mu\nu}^{(n)}$
$\displaystyle=\sum_{klm}^{k+l+m=n}\langle\Phi|[h_{-\mu}^{(k)\dagger},[\tilde{H}^{(l)},h_{-\nu}^{(m)}]]_{+}|\Phi\rangle$
(9) $\displaystyle T_{-p\nu}^{(n)}$
$\displaystyle=\sum_{kl}^{k+l=n}\langle\Phi|[\tilde{a}_{p}^{(k)},h_{-\nu}^{(l)}]_{+}|\Phi\rangle$
(10)

where the ionization operators $h_{-\mu}^{(k)\dagger}$ define the
$(N-1)$-electron basis states
$\ket{\Psi^{(k)}_{-\mu}}=h_{-\mu}^{(k)\dagger}\ket{\Phi}$. For the low-order
EA/IP-ADC($n$) approximations ($n$ $\leq$ 3), only the zeroth- and first-order
components of $h_{\pm\mu}^{(k)\dagger}$ are required and are given as follows:
$h_{+\mu}^{(0)\dagger}=a_{a}^{\dagger}$,
${h}_{+\mu}^{(1)\dagger}=a_{b}^{\dagger}a_{a}^{\dagger}a_{i}$,
$h_{-\mu}^{(0)\dagger}=a_{i}$, and
${h}_{-\mu}^{(1)\dagger}=a_{a}^{\dagger}a_{j}a_{i}$. The $\tilde{H}^{(k)}$ and
$\tilde{a}_{p}^{(k)}$ operators in Eqs. 7, 8, 9 and 10 are obtained by
expanding $\tilde{H}$ and $\tilde{a}_{p}$ using the Baker–Campbell–Hausdorff
(BCH) formula

$\displaystyle\tilde{H}$
$\displaystyle=H^{(0)}+V+[H^{(0)},T^{(1)}-T^{(1){\dagger}}]+[H^{(0)},T^{(2)}-T^{(2){\dagger}}]$
$\displaystyle+\frac{1}{2!}[V+(V+[H^{(0)},T^{(1)}-T^{(1){\dagger}}]),T^{(1)}-T^{(1){\dagger}}]+...$
(11) $\displaystyle\tilde{a}_{p}$
$\displaystyle=a_{p}+[a_{p},T^{(1)}-T^{(1){\dagger}}]+[a_{p},T^{(2)}-T^{(2){\dagger}}]$
$\displaystyle+\frac{1}{2!}[[a_{p},T^{(1)}-T^{(1){\dagger}}],T^{(1)}-T^{(1){\dagger}}]+...$
(12)

and collecting terms at the $k$-th order. These equations depend on the
amplitudes of the excitation operators $T^{(k)}$ (Eq. 6). The low-order EA/IP-
ADC($n$) ($n$ $\leq$ 3) approximations require calculating up to the $n$-th
order single-excitation amplitudes ($t_{i}^{a(k)}$) and up to the ($n-1$)-th
order double-excitation amplitudes ($t_{ij}^{ab(k-1)}$) . At each order $k$,
these amplitudes are computed by solving a system of projected amplitude
equations:

$\displaystyle\langle\Phi|a_{i}^{\dagger}a_{a}\tilde{H}^{(k)}|\Phi\rangle=0$
(13)
$\displaystyle\langle\Phi|a_{i}^{\dagger}a_{j}^{\dagger}a_{b}a_{a}\tilde{H}^{(k)}|\Phi\rangle=0$
(14)

The resulting amplitudes are equivalent to those that define the $k$-th-order
wavefunction in the single-reference Møller–Plesset perturbation theory.

Once the equations for $\textbf{M}_{\pm}$ and $\textbf{T}_{\pm}$ are obtained,
the EA/IP-ADC($n$) transition energies are computed by solving the eigenvalue
problem

$\displaystyle\textbf{M}_{\pm}\textbf{Y}_{\pm}=\textbf{Y}_{\pm}\boldsymbol{\Omega}_{\pm}$
(15)

where $\boldsymbol{\Omega}_{\pm}$ is a diagonal matrix of eigenvalues and Y±
is a matrix of eigenvectors that are used to compute the spectroscopic
amplitudes

$\displaystyle\textbf{X}_{\pm}=\textbf{T}_{\pm}\textbf{Y}_{\pm}$ (16)

which provide information about transition intensities. The elements of
$\textbf{X}_{\pm}$ can be used to calculate the spectroscopic factors

$\displaystyle P_{\pm\mu}=\sum_{p}|X_{\pm p\mu}|^{2}$ (17)

as well as the ADC propagator and density of states

$\displaystyle\mathbf{G}_{\pm}(\omega)$
$\displaystyle=\mathbf{X}_{\pm}\left(\omega\mathbf{1}-\boldsymbol{\Omega}_{\pm}\right)^{-1}\mathbf{X}_{\pm}^{\dagger}$
(18) $\displaystyle A(\omega)$
$\displaystyle=-\frac{1}{\pi}\mathrm{Im}\left[\mathrm{Tr}\,\mathbf{G}_{\pm}(\omega)\right]$
(19)

Working equations for the non-Dyson EA/IP-ADC($n$) ($n$ = 2, 3) methods in the
spin-orbital basis have been presented in Ref. 53.

##  III Implementation

###  III.1 Overview of the EA/IP-ADC implementation

We have implemented the EA/IP-ADC($n$) (n = 2, 3) methods as a module in the
PySCF program.Sun:2020p024109  Our efficient implementation consists of the
spin-restricted (EA/IP-RADC) and unrestricted (EA/IP-UADC) codes, for closed-
and open-shell systems, respectively. The RADC($n$)/UADC($n$) computation is
preceded by the corresponding restricted/unrestricted (RHF/UHF) Hartree–Fock
calculation and follows the algorithm described below:

Table 1: Comparison of formal computational scaling of the EA/IP-ADC($n$) methods ($n\leq 3$) for each step of the algorithm outlined in Section III.1. $O$ and $V$ denote the number of occupied and virtual molecular orbitals, respectively. Also shown are the corresponding computational prefactors relative to those of RADC(3). Component | EA | IP | Relative prefactor
---|---|---|---
|  ADC(2) |  ADC(3) |  ADC(2) |  ADC(3) |  RADC(3) |  UADC(3)
Amplitudes and MP$n$ energy |  $O^{2}V^{3}$ |  $O^{2}V^{4}$ |  $O^{2}V^{3}$ |  $O^{2}V^{4}$ |  1 |  3
Mab(1p-1p) or Mij(1h-1h) |  $O^{2}V^{3}$ |  $O^{2}V^{4}$ |  $O^{3}V^{2}$ |  $O^{2}V^{4}$ |  1 |  3
$\boldsymbol{\sigma}_{\pm}^{(i+1)}=\mathbf{M_{\pm}Y_{\pm}}^{(i)}$ |  $OV^{3}$ |  $O^{2}V^{3}$ |  $O^{3}V$ |  $O^{3}V^{2}$ |  1 |  4
Spectroscopic factors ($P_{\pm\mu}$) |  $O^{2}V^{3}$ |  $O^{2}V^{3}$ |  $O^{3}V^{2}$ |  $O^{3}V^{2}$ |  1 |  3
Overall |  $O^{2}V^{3}$ |  $O^{2}V^{4}$ |  $O^{2}V^{3}$ |  $O^{2}V^{4}$ |  1 |  3

  1. 1.

Transform two-electron repulsion integrals $(pq|rs)$ from the atomic to the
RHF/UHF molecular orbital basis.

  2. 2.

Compute the amplitudes of the effective Hamiltonian ($t_{i}^{a(k)}$ and
$t_{ij}^{ab(k-1)}$, $k\leq n$) by solving Eqs. 13 and 14 and evaluate the
$n$-th-order Møller-Plesset (MP$n$) correlation energy, where $n$ is the order
of the ADC($n$) approximation.

  3. 3.

Compute small blocks of the effective Hamiltonian matrix $\textbf{M}_{\pm}$
with respect to the zeroth-order operators ($h_{\pm\mu}^{(0)\dagger}$, Section
II.1). These are the 1p-1p ($M_{+ab}$, Eq. 7) and 1h-1h ($M_{-ij}$, Eq. 9)
blocks of $\textbf{M}_{+}$ and $\textbf{M}_{-}$ for EA- and IP-ADC($n$),
respectively.

  4. 4.

Define a function for calculating a matrix-vector product
$\boldsymbol{\sigma}_{\pm}=\mathbf{M}_{\pm}\mathbf{Y}_{\pm}$ using the
approximate eigenvectors $\mathbf{Y}_{\pm}$ and the small blocks of the
effective Hamiltonian matrix precomputed in step 3. Solve the eigenvalue
problem (15) by optimizing the eigenvectors $\mathbf{Y}_{\pm}$ until
convergence using an iterative diagonalization algorithm to obtain low-energy
EA/IP’s.

  5. 5.

Using the converged eigenvectors $\mathbf{Y}_{\pm}$, compute spectroscopic
amplitudes $\mathbf{X}_{\pm}$ (Eq. 16) and spectroscopic factors $P_{\pm\mu}$
(Eq. 17). If desired, compute spectral density $A(\omega)$ (Eq. 18), Dyson
orbitals, and other molecular properties.

Our efficient EA/IP-ADC($n$) program written entirely in Python takes
advantage of the highly optimized Basic Linear Algebra Subroutines (BLAS) and
Open Multi-Processing (OpenMP) parallelization for the most expensive tensor
contractions at each step of the algorithm described above. Implementation of
the inexpensive but tedious tensor contractions was performed using the
capabilities of the NumPy package.Harris:2020p357  Our code features three
algorithms for the efficient management of the memory and disk space: (i) in-
core algorithm where all tensors such as two-electron integrals and amplitudes
are stored in memory, (ii) out-of-core algorithm that stores all large tensors
on disk, and (iii) density-fitted algorithm that significantly reduces disk
usage and enables memory savings by approximating two-electron integrals as
discussed in Section III.3. We use the multiroot Davidson
algorithmDavidson:1975p87 ; Liu:1978p49  available in PySCF for iterative
solution of the ADC eigenvalue problem.

Table 1 shows the computational scaling of each step of the ADC algorithm with
the number of occupied ($O$) and virtual ($V)$ molecular orbitals. For $n$ = 2
and 3, the overall computational scaling of EA/IP-ADC($n$) is
$\mathcal{O}(O^{2}V^{3})$ and $\mathcal{O}(O^{2}V^{4})$, respectively.
Importantly, precomputing the 1p-1p and 1h-1h blocks of the effective
Hamiltonian matrix $\mathbf{M}_{\pm}$ at step 3 and reusing them in the
Davidson diagonalization at step 4 lowers the cost of computing the
$\boldsymbol{\sigma}_{\pm}=\mathbf{M}_{\pm}\mathbf{Y}_{\pm}$ matrix-vector
products from $\mathcal{O}(O^{2}V^{4})$ to $\mathcal{O}(O^{2}V^{3})$ and
$\mathcal{O}(O^{3}V^{2})$ for EA- and IP-ADC(3), respectively. As a result,
EA/IP-RADC(3) and EA/IP-UADC(3) require a total of two and six
$\mathcal{O}(O^{2}V^{4})$ contractions, respectively, independent of the
number of iterations in the Davidson algorithm. A similar approach that allows
to reduce the scaling of the matrix-vector products has been employed by
Dempwolff and co-workers in the efficient implementation of IP-
ADC(3).Dempwolff:2019p064108

We note that the computational scaling of solving the EA/IP-ADC($n$)
eigenvalue problem can in principle be lowered if, instead of precomputing the
1p-1p and 1h-1h blocks of the effective Hamiltonian matrix $\mathbf{M}_{\pm}$
at step 3, the contributions of these blocks to the matrix-vector products
(i.e., $\sum_{b}M_{+ab}Y_{b}$ and $\sum_{j}M_{-ij}Y_{j}$) are recomputed at
every iteration of step 4 by forming suitable efficient intermediates. This
avoids step 3 and changes the scaling of step 4 to $\mathcal{O}(O^{2}V^{2})$
for EA/IP-ADC(2) and to $\mathcal{O}(O^{2}V^{3})$/$\mathcal{O}(OV^{4})$ for
EA/IP-ADC(3), respectively, but does not change the overall formal scaling of
the EA/IP-ADC($n$) methods due to the higher computational scaling of step 2.

In the following, we discuss the derivation of the spin-adapted equations used
in our EA/IP-RADC(3) implementation (Section III.2) and give an overview of
the density fitting approximation used for the two-electron integrals (Section
III.3). The efficiency of our implementation and the accuracy of the density
fitting approximation are analyzed in Section V.1.

###  III.2 Spin-restricted EA/IP-ADC for closed-shell systems (EA/IP-RADC)

In our earlier work,Banerjee:2019p224112  we presented a general spin-orbital
derivation and implementation of EA/IP-ADC($n$) ($n\leq 3$) applicable to
closed- and open-shell systems. Here, we discuss spin adaptation of the EA/IP-
ADC equations used in our efficient spin-restricted EA/IP-RADC implementation.
While spin adaptation of the EA-ADC equations has not been discussed
previously in the literature, it has been presented for other ADC
formulations, such as ADC for ionizationThiel:2003p2088  and
neutralBarth:1985p867 ; Trofimov:2002p6402  excitation energies.

In our derivation of the EA/IP-RADC equations we follow the approach described
in Ref. 67. Starting with spin-orbital EA/IP-ADC equations, we obtain the
spin-adapted equations for the effective Hamiltonian ($\mathbf{M}_{\pm}$) and
effective transition moments ($\mathbf{T}_{\pm}$) matrix elements in Eqs. 7,
8, 9 and 10 by integrating over the spin variables. We denote each spin-
orbital with an index $p_{\sigma}$, where $\sigma=\\{\alpha,\beta\\}$ is a
spin label and $p$ is an index of the spatial molecular orbital (MO) from the
restricted Hartree–Fock (RHF) calculation. We use the $i,j,k,l,\ldots$ and
$a,b,c,d,\ldots$ indices to label the occupied and virtual MO’s, respectively,
and reserve $p,q,r,s,\ldots$ to denote a general (occupied or virtual) MO.

Integrating over the spin variables allows to express the spin-orbital
$\mathbf{M}_{\pm}$ and $\mathbf{T}_{\pm}$ matrix elements in terms of their
counterparts that depend solely on spatial orbitals. As an example, the matrix
elements of $\mathbf{M}_{-}$ and $\mathbf{T}_{-}$ in IP-ADC can be expressed
as:

$\displaystyle M_{-i_{\alpha},j_{\alpha}}$
$\displaystyle=M_{-i_{\beta},j_{\beta}}=M_{-i,j}$ (20) $\displaystyle
M_{-i_{\alpha},j_{\beta}k_{\alpha}a_{\beta}}$
$\displaystyle=M_{-i_{\beta},j_{\alpha}k_{\beta}a_{\alpha}}=M_{-i,jka}$ (21)
$\displaystyle M_{-i_{\alpha},j_{\alpha}k_{\alpha}a_{\alpha}}$
$\displaystyle=M_{-i_{\beta},j_{\beta}k_{\beta}a_{\beta}}$
$\displaystyle=M_{-i,jka}-M_{-i,kja}$ (22) $\displaystyle
M_{-i_{\beta}l_{\alpha}b_{\beta},j_{\beta}k_{\alpha}a_{\beta}}$
$\displaystyle=M_{-i_{\alpha}l_{\beta}b_{\alpha},j_{\alpha}k_{\beta}a_{\alpha}}=M_{-ilb,jka}$
(23) $\displaystyle
M_{-i_{\alpha}l_{\alpha}b_{\alpha},j_{\alpha}k_{\alpha}a_{\alpha}}$
$\displaystyle=M_{-i_{\beta}l_{\beta}b_{\beta},j_{\beta}k_{\beta}a_{\beta}}$
$\displaystyle=M_{-ilb,jka}-M_{-lib,jka}$ $\displaystyle-
M_{-ilb,kja}+M_{-lib,kja}$ (24) $\displaystyle T_{-p_{\alpha},i_{\alpha}}$
$\displaystyle=T_{-p_{\beta},i_{\beta}}=T_{-p,i}$ (25) $\displaystyle
T_{-p_{\alpha},j_{\beta}k_{\alpha}a_{\beta}}$
$\displaystyle=T_{-p_{\beta},j_{\alpha}k_{\beta}a_{\alpha}}=T_{-p,jka}$ (26)
$\displaystyle T_{-p_{\alpha},j_{\alpha}k_{\alpha}a_{\alpha}}$
$\displaystyle=T_{-p_{\beta},j_{\beta}k_{\beta}a_{\beta}}$
$\displaystyle=T_{-p,jka}-T_{-p,kja}$ (27)

where our notation for the $\mathbf{M}_{-}$ and $\mathbf{T}_{-}$ matrix
elements matches that of Ref. 53.

To obtain working equations for $\mathbf{M}_{-}$ and $\mathbf{T}_{-}$ in IP-
RADC, we write their spin-adapted matrix elements in terms of the spin-orbital
Fock matrix eigenvalues ($\varepsilon_{p_{\sigma}}$), antisymmetrized two-
electron integrals
($\braket{p_{\sigma_{1}}q_{\sigma_{2}}}{r_{\sigma_{3}}s_{\sigma_{4}}}$), and
$k$-th-order amplitudes of the effective Hamiltonian
($t_{i_{\sigma_{1}}}^{a_{\sigma_{2}}(k)}$ and
$t_{i_{\sigma_{1}}j_{\sigma_{2}}}^{a_{\sigma_{3}}b_{\sigma_{4}}(k)}$, Eq. 6),
and perform spin-integration as follows:

$\displaystyle\varepsilon_{p_{\alpha}}$
$\displaystyle=\varepsilon_{p_{\beta}}=\varepsilon_{p}$ (28)
$\displaystyle\langle p_{\alpha}q_{\beta}||r_{\alpha}s_{\beta}\rangle$
$\displaystyle=\langle
p_{\beta}q_{\alpha}||r_{\beta}s_{\alpha}\rangle=(pr|qs)$ (29)
$\displaystyle\langle p_{\alpha}q_{\beta}||r_{\beta}s_{\alpha}\rangle$
$\displaystyle=\langle
p_{\beta}q_{\alpha}||r_{\alpha}s_{\beta}\rangle=-(ps|qr)$ (30)
$\displaystyle\langle p_{\alpha}q_{\alpha}||r_{\alpha}s_{\alpha}\rangle$
$\displaystyle=\langle
p_{\beta}q_{\beta}||r_{\beta}s_{\beta}\rangle=(pr|qs)-(ps|qr)$ (31)
$\displaystyle t_{i_{\alpha}}^{a_{\alpha}(k)}$
$\displaystyle=t_{i_{\beta}}^{a_{\beta}(k)}=t_{i}^{a(k)}$ (32) $\displaystyle
t_{i_{\alpha}j_{\beta}}^{a_{\alpha}b_{\beta}(k)}$
$\displaystyle=t_{i_{\beta}j_{\alpha}}^{a_{\beta}b_{\alpha}(k)}=t_{ij}^{ab(k)}$
(33) $\displaystyle t_{i_{\alpha}j_{\beta}}^{a_{\beta}b_{\alpha}(k)}$
$\displaystyle=t_{i_{\beta}j_{\alpha}}^{a_{\alpha}b_{\beta}(k)}=-t_{ij}^{ba(k)}$
(34) $\displaystyle t_{i_{\alpha}j_{\alpha}}^{a_{\alpha}b_{\alpha}(k)}$
$\displaystyle=t_{i_{\beta}j_{\beta}}^{a_{\beta}b_{\beta}(k)}=t_{ij}^{ab(k)}-t_{ij}^{ba(k)}$
(35)

where $(pq|rs)$ denotes a Chemists’ notation two-electron integral in the RHF
MO basis. The resulting equations are used to solve the IP-RADC eigenvalue
problem by computing the matrix-vector products

$\displaystyle\sigma_{-i}$
$\displaystyle=\sum_{j}M_{-i,j}Y_{-j}+2\sum_{jka}M_{-i,jka}Y_{-jka}$
$\displaystyle-\sum_{jka}M_{-i,jka}Y_{-kja}$ (36) $\displaystyle\sigma_{-jka}$
$\displaystyle=\sum_{i}M_{-i,jka}Y_{-i}+\sum_{ilb}M_{-jka,ilb}Y_{-ilb}$ (37)

and iteratively optimizing the spin-adapted eigenvectors with elements
$\mathbf{Y_{-}}=\\{Y_{-i},Y_{-ija}\\}$. The converged eigenvectors are then
used to evaluate the spectroscopic factors as described in the Appendix.
Expressions for $\mathbf{M}_{+}$, $\mathbf{T}_{+}$, and
$\boldsymbol{\sigma}_{+}$ in EA-ADC can be obtained by replacing the 1h ($i$)
and 2h-1p ($ija$) excitation labels with those of the 1p ($a$) and 2p-1h
($abi$) excitations.

Spin adaptation of the EA/IP-ADC equations provides significant computational
savings for all steps of the ADC algorithm outlined in Section III.1, with a
three-fold reduction in the formal computational prefactor, as shown in Table
1. Most significant savings are achieved for the solution of the EA/IP-ADC
eigenvalue problem, where the matrix-vector products
$\boldsymbol{\sigma}_{\pm}$ are expressed in terms of a much smaller number of
the spin-adapted $\mathbf{M}_{\pm}$ matrix elements (e.g., $M_{-i,j}$,
$M_{-i,jka}$, and $M_{-ilb,jka}$ in Eqs. 20, 21, 22, 23 and 24) compared to
that of the spin-orbital implementation. We will analyze the computational
efficiency of our EA/IP-RADC implementation in Section V.1.

###  III.3 Density fitting in EA/IP-ADC

In addition to using the conventional two-electron integrals ($(pq|rs)$ or
$\braket{p_{\sigma_{1}}q_{\sigma_{2}}}{r_{\sigma_{3}}s_{\sigma_{4}}}$), our
EA/IP-RADC and UADC implementations feature support of the density fitting
(DF) approximationWhitten:1973p4496 ; Dunlap:1979p3396 ; Vahtras:1993p514 ;
Feyereisen:1993p359  that expresses the four-index integrals as a product of
the two- and three-index tensors. For example, the two-electron integrals in
EA/IP-RADC can be approximated as:

$\displaystyle(pq|rs)\approx\sum_{Q}^{N_{aux}}b_{pq}^{Q}b_{rs}^{Q}$ (38)

where the $b_{pq}^{Q}$ coefficients are defined in terms of the MO’s
($\phi_{p}$) and the auxiliary basis functions ($\chi_{P}$) as

$\displaystyle b_{pq}^{Q}$
$\displaystyle=\sum_{P}^{N_{aux}}(pq|P)(J^{-\frac{1}{2}})_{PQ}$ (39)
$\displaystyle(pq|P)$
$\displaystyle=\iint\phi_{p}(1)\phi_{q}(1)\frac{1}{r_{12}}\chi_{P}(2)dr_{1}dr_{2}$
(40) $\displaystyle J_{PQ}$
$\displaystyle=\iint\chi_{P}(1)\frac{1}{r_{12}}\chi_{Q}(2)dr_{1}dr_{2}$ (41)

The use of the DF approximation has virtually no effect on the results of the
EA/IP-ADC calculations (see Section V.1 for details) while offering two
significant advantages: i) it greatly reduces the disk usage and the
input/output operation count by avoiding the storage of large four-index
tensors on disk, and ii) it dramatically lowers the computational scaling of
the two-electron integral transformation from $\mathcal{O}(N^{5})$ to
$\mathcal{O}(N_{aux}N^{3})$ (where $N$ and $N_{aux}$ are the numbers of
$\phi_{p}$ and $\chi_{P}$, respectively), providing significant computational
savings for large molecules and basis sets. In our implementation, we store
the $b_{pq}^{Q}$ tensors either in memory or on disk and recompute sectors of
$(pq|rs)$ on demand. Similar to other methods that employ the DF
approximation,Werner:2003p8149 ; Bozkaya:2014er ; Epifanovsky:2013gd ;
Wang:2016p4833  our density-fitted EA/IP-ADC methods support using the JKFIT
and RI auxiliary basis sets. The RI basis functions are used to approximate
two-electron integrals in the EA/IP-ADC calculations, while the JKFIT
auxiliary basis set optimized for fitting the Coulomb and exchange integrals
in the mean-field (SCF) computations is used in the reference Hartree–Fock
calculation. We note that the DF approximation has been previously used in
other efficient single-reference ADC implementations.Helmich:2014p35 ;
Mester:2018p094111 ; Herbst:2020pe1462 ; Liu:2020p174109

##  IV Computational details

All EA/IP-ADC($n$) calculations ($n\leq 3$) were performed using the ADC
module in PySCFSun:2020p024109 ; Banerjee:2019p224112  and used the singly- or
doubly-augmented Dunning’s correlation consistent basis sets denoted as aug-
cc-pV$X$Z or d-aug-cc-pV$X$Z, respectively.Dunning:1989p1007 ;
Kendall:1992p6796 ; Woon:1994p2975  Specifically, in Section V.1, we used aug-
cc-pVQZ to benchmark the accuracy of the density fitting approximation and
aug-cc-pVDZ to test computational efficiency of our implementation. All ADC(2)
and MP2 computations of the TEMPO radical (Section V.2) and DNA base pairs
(Section V.3) employed the aug-cc-pVTZ basis sets. For the ADC(3) and MP3
calculations we used a modified aug-cc-pVTZ basis set where the diffuse
functions were removed from the hydrogen atoms. We denote this basis set as
aug(nH)-cc-pVTZ. Additionally, for the DNA base pairs, we employed the doubly-
augmented d-aug-cc-pVDZ basis set.

Geometries of molecules in Section V.1 were optimized using coupled cluster
theory with single, double, and perturbative triple excitations
(CCSD(T)).Raghavachari:1989p479 ; Bartlett:1990p513  The EA/IP-ADC($n$)
calculations ($n\leq 3$) of the TEMPO radical and anion used structures
optimized using the $n$-th-order Møller-Plesset perturbation
theoryMoller:1934kp618  (MP$n$) and the same basis set. The EA-ADC($n$)
calculations for the base pairs and their corresponding anions were performed
with the MP2/aug-cc-pVTZ optimized structures. All geometry optimizations were
performed using the PSI4 programSmith:2020p184108  and the optimized
structures are reported in the Supporting Information. All EA/IP-ADC($n$) and
MP$n$ calculations employed the density fitting approximation. The correlation
MP$n$ energies and ADC($n$) excitation energies were computed using the
resolution-of-identity (aug-cc-pV$X$Z-RI) auxiliary basis sets, while the SCF
energies and molecular orbitals were obtained using the aug-cc-pV$X$Z-JKFIT
auxiliary basis sets,Weigend:2002p4285 ; Weigend:2002p3175 ; Hattig:2005p59
where $X$ corresponded to that of the main one-electron basis set.

Throughout the manuscript, positive electron affinity implies exothermic
electron attachment (i.e., EA = $E_{N}$ $-$ $E_{N+1}$), whereas a positive
ionization energy corresponds to an endothermic process (IP = $E_{N-1}$ $-$
$E_{N}$). While the vertical electron affinities (VEA’s) were obtained
directly from our EA-ADC($n$) calculations, the corresponding adiabatic
electron attachment (AEA’s) energies were determined as follows:

$\displaystyle AEA=VEA+\Delta E_{MPn}$ (42)

where $\Delta E_{MPn}=E^{N+1}_{MPn}$(neutral) $-$ $E^{N+1}_{MPn}$(anion) is
the difference of the MP$n$ energies of the anionic system computed at the
neutral and anion optimized geometries. The TEMPO radical AEA’s have been
corrected to account for the zero-point vibrational energy (ZPVE)

$\displaystyle AEA(ZPVE)$ $\displaystyle=AEA$
$\displaystyle+ZPVE(N)-ZPVE(N+1)$ (43)

where $ZPVE(N)$ and $ZPVE(N+1)$ are the ZPVE’s of the neutral and anionic
TEMPO, respectively, computed at the MP2/aug-cc-pVDZ level of theory.

The nature of the electron-attached states was analyzed by visualizing the
Dyson orbitals,Melania:2007p234106  which are computed from the overlap of the
ground state wavefunction and the electron-attached/ionized states:

$\displaystyle|\phi^{Dyson}_{\pm\mu}\rangle=\sum_{p}X_{\pm
p\mu}|\phi_{p}\rangle$ (44)

where $|\phi_{p}\rangle$ are the reference Hartree–Fock orbitals and $X_{\pm
p\mu}$ are the EA/IP-ADC spectroscopic amplitudes defined in Eq. 16. All Dyson
orbitals were computed using the aug-cc-pVDZ basis set.

Figure 1: Mean absolute errors (eV) of the density fitting approximation in
the vertical electron attachment or ionization energies for EA-RADC($n$) (a),
IP-RADC($n$) (b), EA-UADC($n$) (c), and IP-UADC($n$) (d) with $n=2,3$ (aug-cc-
pVQZ basis set). For each method, two approximations are considered, with (DF-
ADC($n$)/DF-SCF) and without (DF-ADC($n$)/SCF) density fitting in the
reference SCF calculation. The RADC($n$) and UADC($n$) results were obtained
for sets of five closed-shell and five open-shell molecules, respectively. See
Tables S1 and S2 in the Supporting Information for data on individual
molecules.

The TEMPO radical photoelectron spectra were simulated by plotting the IP-ADC
density of states

$\displaystyle A(\omega)$
$\displaystyle=-\frac{1}{\pi}\mathrm{Im}\left[\sum_{\mu}\frac{P_{-\mu}}{\omega-\omega_{\mu}+i\eta}\right]$
(45)

for a range of frequencies $\omega$, where $\omega_{\mu}$ denotes the
ionization energies obtained by solving Eq. 15, $\eta$ is a small broadening,
and $P_{-\mu}$ are the spectroscopic factors computed using Eq. 17. The
simulated spectra were compared to the experimental photoelectron spectrum
that was digitized using the WebPlotDigitizer program.Rohatgi2020

##  V Results

###  V.1 Accuracy of the density fitting approximation and analysis of the
computational efficiency

In this section, we benchmark the accuracy of the density fitting (DF)
approximation for EA/IP-RADC($n$) ($n$ = 2, 3) and analyze the efficiency of
our restricted (RADC) and unrestricted (UADC) implementations. Figure 1 shows
the errors in the vertical electron attachment (EA) and ionization (IP)
energies computed using the density-fitted EA- or IP-ADC($n$) methods,
respectively, relative to those of the conventional EA/IP-ADC($n$) ($n$ = 2,
3) methods without the density fitting approximation. For each EA/IP-DF-
ADC($n$) method, results are presented for sets of five closed-shell (H2O, HF,
N2, SiH2, H2CO) and five open-shell (OH, NH, O2, NH2, CH3) molecules computed
using the RADC and UADC implementations, respectively (with data for
individual molecules available in Tables S1 and S2 of the Supporting
Information). Furthermore, we consider two types of the DF approximations,
with and without density fitting in the reference SCF calculation, denoted as
DF-ADC($n$)/DF-SCF and DF-ADC($n$)/SCF, respectively.

For all methods considered, the errors of the DF approximation are very small
and do not exceed 0.0021 eV. The largest mean absolute errors (MAE’s) of
$\sim$ 0.001 eV are observed for IP-DF-ADC(3) of closed- and open-shell
molecules and for EA-DF-ADC(3) of open-shell molecules. The EA/IP-DF-ADC(2)
methods consistently produce smaller errors with MAE’s less than 0.0005 eV.
Incorporating the DF approximation in the SCF reference calculations increases
the computed MAE’s by no more than 0.0003 eV. Overall, we find that the DF
approximation in EA/IP-DF-ADC($n$) ($n$ = 2, 3) is very accurate, which is
consistent with findings for other quantum chemical methods when computing
relative energies.Epifanovsky:2013gd ; Bozkaya:2014p2371 ; Wang:2016p4833 ;
Bozkaya:2016p1179 ; Hattig:2000p5154

Figure 2: Comparison of the computational wall time (in sec) for the
restricted (R) and unrestricted (U) implementations of the conventional (out-
of-core) and density-fitted (DF) ADC algorithms for EA-ADC(3) (a) and IP-
ADC(3) (b). Data is shown for the cytosine molecule with the aug-cc-pVDZ basis
(229 basis functions, 3 lowest EA or IP). Wall time for each step of the ADC
algorithm outlined in Section III.1 is shown with a different color. All
calculations were performed on a single computer with a 12-core Intel Xeon
Gold 6136 3.00GHz CPU.

We now analyze the efficiency of our EA/IP-RADC(3) implementations by
comparing the computational wall time for calculating the three lowest-energy
EA or IP of the cytosine molecule with the aug-cc-pVDZ basis set (229 basis
functions). We employ four algorithms for each method: i) out-of-core UADC(3),
ii) density-fitted UADC(3), iii) out-of-core RADC(3), and iv) density-fitted
RADC(3). The computed wall time is broken down into contributions from each
step of the ADC algorithm outlined in Section III.1.

Taking advantage of the spin adaptation in the RADC algorithm reduces the
total time of the EA- and IP-ADC(3) calculations by approximately a factor of
three. Significant savings are observed for all steps of the ADC computation,
with $\sim$ 2.5-fold speedup for calculating the amplitudes of the effective
Hamiltonian and the 1p-1p/1h-1h blocks of the $\mathbf{M}_{\pm}$ matrix and
more than four-fold reduction in time for calculating the
$\boldsymbol{\sigma}_{\pm}$ vector and spectroscopic factors, in a good
agreement with theoretical prefactors reported in Table 1. Introducing the DF
approximation drastically lowers the computational cost of the two-electron
integral transformation, which has a substantial effect on the wall time of
the EA/IP-UADC(3) calculations resulting in a 1.5- to 2-fold speedup. Timings
for other steps of the ADC algorithm are not affected by the density fitting.
Combining the spin adaptation and density fitting allows to achieve nearly a
four-fold speedup for both EA- and IP-ADC(3), pushing the limits of their
applicability to larger chemical systems and basis sets. In Sections V.2 and
V.3, we demonstrate this by applying EA/IP-RADC($n$) ($n$ = 2, 3) to the TEMPO
radical and the DNA base pairs.

###  V.2 TEMPO radical

Figure 3: Molecular structure of the TEMPO radical.

We first study the (2,2,6,6-tetramethylpiperidin-1-yl)oxyl radical, commonly
known as TEMPO (Figure 3), that is widely used as a catalyst in organic
synthesis, as a mediator in polymerization reactions, and as a structural
probe for performing electron spin resonance spectroscopy on biological
systems.Barriga:2001p0563 ; Galli:2008p705 ; Likhtenshtein:2008nitroxides
TEMPO and its derivatives are also used in solar cells, lithium batteries, and
radical sensors.Nakahara:2011p227 ; Janoschka:2015p78 ; Nakahara:2006p921 ;
Blinco:2008p6763 ; Kato:2010p464 ; Buhrmester:2006pa1800

Figure 4: Density of states of the TEMPO radical computed using IP-
UADC(2)/aug-cc-pVTZ (a) and IP-UADC(3)/aug(nH)-cc-pVTZ (b) with a broadening
of 0.15 eV, compared to the experimental photoelectron spectrum from Ref. 103.
Densities of states were shifted by 1.02 ($n$ = 2) and $-$0.3 eV ($n$ = 3) to
reproduce the position of the first peak in the experimental spectrum.

Figure 5: Dyson orbitals for to the three lowest-energy ionization transitions
of the neutral TEMPO radical corresponding to the $1^{1}A^{\prime}$ (a),
$1^{3}A^{\prime\prime}$ (b) and $1^{1}A^{\prime\prime}$ (c) states of the
TEMPO cation computed using the IP-UADC(2) method with the aug-cc-pVTZ basis
set.

The diverse applications of the TEMPO radical have motivated several
experimental and theoretical investigations of its electronic structure.
Photoionization of TEMPO has been thoroughly studied, experimentally and
theoretically, by Ljubić, Novak, and co-workers using UV/Vis and X-ray
photoelectron spectroscopy.Novak:2004p7628 ; Kovavc:2014p10734 ;
Ljubic:2014p2333 ; Ljubic:2016p10207  These studies demonstrated a good
agreement between the experimental photoelectron spectra and theoretical
spectra simulated using density functional theory (DFT). However, electron
attachment to TEMPO has not been investigated as extensively as its
photoionization. Experimentally, the TEMPO anion has been studied by Kubala et
al. using dissociative electron attachment spectroscopy, electron energy-loss
spectroscopy, and vibrational excitation cross section
measurements.Kubala:2013p2033  Their results indicated that the energy of the
electron-attached state is either slightly below or slightly above that of the
neutral TEMPO radical, but did not allow to distinguish between these two
possibilities. These findings were supported by the DFT calculations
(B3LYP/6-311++G(2d,p)) that predicted the endothermic vertical electron
affinity (VEA) of $-$0.3 eV and exothermic zero-point-energy corrected
adiabatic electron affinity (AEA) of 0.14 eV. In another theoretical study,
Sohn et al. used a combination of the B3LYP and MP2 methods with the
6-311++G(3df,2pd) basis set and computed AEA to be 0.64 eV.Sohn:2015p49  A
more recent study by Gunasekara et al. using coupled cluster theory with
single, double, and perturbative triple excitations (CCSD(T)) and the aug-cc-
pVTZ basis set reports AEA of 0.44 eV.Gunasekara:2019p1882

Table 2: Vertical (VEA) and adiabatic (AEA) electron attachment energies (in
eV) of the TEMPO radical computed using MP$n$ and EA-UADC($n$) methods. Also
shown are AEA calculated with the zero-point vibrational energy (ZPVE)
correction and best available results from other theoretical methods.

Method | VEA | AEA | AEA (ZPVE-corrected)
---|---|---|---
MP2/aug-cc-pVTZ |  $-$0.21 |  0.57 |  0.70
MP3/aug(nH)-cc-pVTZ |  $-$0.58 |  0.23 |  0.35
EA-UADC(2)/aug-cc-pVTZ |  $-$0.30 |  0.49 |  0.61
EA-UADC(3)/aug(nH)-cc-pVTZ |  $-$0.72 |  0.08 |  0.21
Reference |  $-$0.3a |  0.64b, 0.44c |  0.14a

  * a

B3LYP/6-311++G(2d,p) from Ref. 103.

  * b

MP2/6-311++G(3df,2pd) from Ref. 104.

  * c

CCSD(T)/aug-cc-pVTZ from Ref. 105.

To check the accuracy of our ADC methods for the TEMPO radical, we first
compute its ionization spectrum using IP-UADC($n$) ($n$ = 2, 3) and compare it
to the experimental photoelectron spectrum reported by Kubala et
al.Kubala:2013p2033  Figure 4 shows the density of occupied states (DOS) of
TEMPO computed using IP-UADC($n$) ($n$ = 2, 3) with a broadening of 0.15 eV.
For each method, DOS was shifted to align the first IP with the maximum of the
first peak in the experimental spectrum. The IP-UADC(3) DOS required a smaller
shift ($-0.3$ eV) compared to that of IP-UADC(2) (1.02 eV), indicating a
smaller error of the third-order approximation in the computed IP’s, due to
its higher-level description of electron correlation effects. The shifted DOS
computed using both IP-UADC(2) and IP-UADC(3) agree well with the experimental
photoelectron spectrum, accurately reproducing peak spacings and relative peak
intensities, with IP-UADC(3) showing a better agreement than IP-ADC(2) for
ionization energies greater than 10 eV. Upon shifting, both methods predict a
signal at $\sim$ 7.2 eV, corresponding to the $X^{2}A^{\prime}\rightarrow
1^{1}A^{\prime}$ photoelectron transition, and two peaks at $\sim$ 9.0 and 9.3
eV that can be assigned as the $X^{2}A^{\prime}\rightarrow
1^{3}A^{\prime\prime}$ and $X^{2}A^{\prime}\rightarrow 1^{1}A^{\prime\prime}$
transitions, respectively. Figure 5 shows the Dyson orbitals for each of these
photoionization transitions computed using IP-UADC(2). All three Dyson
orbitals are localized on the NO group of the TEMPO radical (Figure 3), either
in or out of the plane of its six-membered ring. The out-of-plane orbital
shows a greater localization on the N atom, while the in-plane orbitals are
more localized on the O atom. While IP-UADC(2) predicts both
$X^{2}A^{\prime}\rightarrow 1^{3}A^{\prime\prime}$ and
$X^{2}A^{\prime}\rightarrow 1^{1}A^{\prime\prime}$ transitions to have similar
spectral density, IP-UADC(3) yields a smaller spectral density for the higher-
energy $X^{2}A^{\prime}\rightarrow 1^{1}A^{\prime\prime}$ transition, in a
good agreement with the experiment. The IP-UADC(3) DOS also shows a weak
satellite transition at $\sim$ 10 eV, which is not visible in the experimental
photoelectron spectrum.

We now turn our attention to VEA and AEA of the TEMPO radical computed using
the MP$n$ and EA-UADC($n$) methods ($n$ = 2, 3) and presented in Table 2. In
agreement with a previous DFT study,Kubala:2013p2033  all MP$n$ and EA-
UADC($n$) methods predict a negative VEA ranging between $-$0.72 eV (EA-
UADC(3)) to $-$0.21 eV (MP2), indicating an endothermic vertical electron
attachment. The VEA computed using MP$n$ and EA-UADC($n$) are close to each
other at each order in perturbation theory (within 0.15 eV) and become more
negative as the order increases from $n$ = 2 to 3.

Upon structural relaxation, the electron-attached state of the TEMPO radical
lowers its energy below that of the neutral radical and the electron
attachment becomes energetically favorable, as evidenced by the positive AEA
computed by all four MP$n$ and EA-UADC($n$) approximations. Out of all
methods, the AEA of EA-UADC(2) (0.49 eV) shows the best agreement with the AEA
from the CCSD(T) study by Gunasekara and co-workers (0.44
eV),Gunasekara:2019p1882  which used the same aug-cc-pVTZ basis set, but a
different (DFT-optimized) geometry. The EA-UADC(3) method predicts a small
positive AEA of 0.08 eV, in a good agreement with experimental and theoretical
findings of Kubala et al.Kubala:2013p2033  Table 2 also reports AEA corrected
by the zero-point vibrational energies (ZPVE) of the neutral and anionic
states (see Section IV for details). Including the zero-point vibrational
effects increases the EA-UADC(3) AEA to 0.21 eV, which agrees well with the
experimental results.

Figure 6: Dyson orbitals corresponding to the lowest-energy electron
attachment to the TEMPO radical at the neutral (a) and anion (b) equilibrium
geometries computed using EA-UADC(3) with the aug-cc-pVDZ basis set.

Figures 6 and 6 show the EA-UADC(3) Dyson orbitals of the TEMPO radical
computed at the equilibrium geometries of the neutral and electron-attached
states. At the neutral equilibrium geometry, the Dyson orbital is localized
outside of the molecule, indicating a dipole-bound nature of the electron-
attached state (Figure 6). Once the TEMPO geometry is relaxed to that of its
anion, the lowest-energy electron-attached state acquires the valence-bound
character (Figure 6) where most of the orbital is localized on the NO group of
the radical.

###  V.3 DNA base pairs

Figure 7: Molecular structures of the guanine–cytosine (a) and adenine–thymine
(b) base pairs.

Finally, we apply our efficient density-fitted EA/IP-ADC implementation to
compute VEA and AEA of two DNA nucleotide base pairs: guanine–cytosine (GC)
and adenine–thymine (AT), shown in Figure 7. Investigation of electron
attachment to nucleobases is important for understanding light-induced
biochemical processes, such as electron or hole transfer along the DNA/RNA
strands that can eventually lead to genetic mutation.Teoule:1987p573 ;
Barrios:2002p7991 ; Boudaiffa:2000p1658  In our previous
work,Banerjee:2019p224112  we reported VEA’s for the five isolated nucleobases
[adenine (A), cytosine (C), guanine (G), thymine (T), and uracil (U)] computed
using the EA-ADC($n$) ($n$ = 2, 3), CCSD(T), and EOM-CCSD methods. The
calculated VEA’s of the nucleobases decreased in the order G $>$ U $\gtrsim$ C
$\gtrsim$ T $>$ A, in a good agreement with the G $>$ U $\gtrsim$ T $\gtrsim$
C $>$ A trend from the experiment.Aflatooni:1998p6205  However, the
nucleobases primarily appear in pairs, linked via a network of hydrogen bonds
that can alter the electron affinities of the nucleobases, making it difficult
to predict the site of electron localization in the base
pairs.Colson:1992p9787 ; Richardson:2002p10163

Although no experimental electron affinities have been reported for the
isolated base pairs, several theoretical studies have been
performed.Colson:1992p9787 ; Richardson:2002p10163 ; Richardson:2003p848 ;
Reynisson:2002p5353 ; Al:2000p2994 ; Smets:2001p342 ; Peng:2012p120 ;
Tripathi:2019p10131  Early computational investigations using Hartree–Fock and
density functional theories (DFT)Colson:1992p9787 ; Richardson:2002p10163 ;
Richardson:2003p848  reported that the electron attachment to the GC base pair
preferentially localizes the electron on the C nucleobase, violating the G $>$
C trend observed for the isolated nucleobases. For the AT base pair, the
electron was found to localize on T, in agreement with the T $>$ A trend in
the experimental EA’s. Another DFT study reported A to be the primary site for
the electron localization in the AT anion.Reynisson:2002p5353  Electron
attachment to GC and AT was studied by the AdamowiczAl:2000p2994 ;
Smets:2001p342  and SchaeferPeng:2012p120  groups using the MP2 method with
the double-zeta basis sets incorporating additional diffuse functions. They
reported AEA’s of the GC/AT pairs to be $-$0.06/$-$0.4 eV and 0.01/$-$0.42 eV,
respectively. Both of these studies found the evidence of electron
localization on A in the AT anion. Very recently, Tripathi et al. performed a
study of electron attachment to the GC and AT base pairs using the EOM-CCSD
method.Tripathi:2019p10131  Using the aug-cc-pVTZ basis set with additional
5s5p4d diffuse functions, they reported the computed VEA’s for the GC/AT pair
to be 0.099/0.003 eV and the AEA’s of 0.237/$-$0.014 eV, respectively. The
authors of this study also demonstrated that at the equilibrium geometry of
the neutral base pairs the electron attachment occurs in the dipole-bound
state, while it occurs in the valence-bound state at the equilibrium
geometries of the AT and GC anions.

Figure 8: Dyson orbitals for the lowest-energy electron attachment to the
guanine–cytosine base pair at the neutral (a) and anion (b) equilibrium
geometries computed using the EA-ADC(3) method with the aug-cc-pVDZ basis set.

Figure 9: Dyson orbitals for the lowest-energy electron attachment to the
adenine–thymine base pair at the neutral (a) and anion (b) equilibrium
geometries computed using the EA-ADC(3) method with the aug-cc-pVDZ basis set.

Figures 8 and 9 show the EA-ADC(3) Dyson orbitals for the lowest-energy
electron attachment to the GC and AT base pairs computed at the neutral (a)
and anion (b) equilibrium geometries. In agreement with the results of
Tripathi et al.,Tripathi:2019p10131  the Dyson orbitals computed at the
neutral equilibrium geometries are localized outside of each molecule, showing
the formation of the dipole-bound states. Relaxation to the anion equilibrium
geometries leads to a formation of the valence-bound states with the Dyson
orbitals localized on C of GC (Figure 8) and T of AT (Figure 9), which agrees
well with the results of other post-Hartree–Fock calculations reported in the
literature.Al:2000p2994 ; Smets:2001p342 ; Peng:2012p120 ; Tripathi:2019p10131

Table 3: Vertical (VEA) and adiabatic (AEA) electron attachment energies (in
eV) of the guanine–cytosine and adenine–thymine base pairs computed using
MP$n$ and EA-ADC($n$) ($n$ = 2, 3) methods. Also shown are the EOM-CCSD/aug-
cc-pV$X$Z+5s5p4d ($X$ = D, T) reference results from Ref. 117.

Method | Guanine–Cytosine | Adenine–Thymine
---|---|---
|  VEA |  AEA |  VEA |  AEA
MP2/aug-cc-pVTZ |  $-$0.43 |  0.07 |  $-$0.69 |  $-$0.29
MP3/aug(nH)-cc-pVTZ |  $-$0.40 |  0.09 |  |
EA-ADC(2)/aug-cc-pVTZ |  0.19 |  0.69 |  $-$0.06 |  0.34
EA-ADC(3)/aug(nH)-cc-pVTZ |  $-$0.34 |  0.15 |  $-$0.54 |  $-$0.14a
MP2/d-aug-cc-pVDZ |  0.06 |  0.05 |  $-$0.07 |  $-$0.33
MP3/d-aug-cc-pVDZ |  0.06 |  0.10 |  $-$0.07 |  $-$0.41
EA-ADC(2)/d-aug-cc-pVDZ |  0.14 |  0.14 |  0.02 |  $-$0.24
EA-ADC(3)/d-aug-cc-pVDZ |  0.09 |  0.13 |  $-$0.03 |  $-$0.38
EOM-CCSD/aug-cc-pVDZ+5s5p4db |  0.10 |  0.18 |  0.00 |  $-$0.10
EOM-CCSD/aug-cc-pVTZ+5s5p4db |  0.10 |  0.24 |  0.00 |  $-$0.01

  * a

The AEA was obtained using Eq. 42 where $\Delta E_{MPn}$ was computed using
MP2/aug-cc-pVTZ.

  * b

Ref. 117.

Table 3 presents the vertical (VEA) and adiabatic (AEA) electron affinities of
GC and AT computed using the second- and third-order MP$n$ and ADC($n$)
methods. Due to the dipole-bound nature of the vertically-attached electronic
states, the computed VEA’s are expected to exhibit a strong basis set
dependence. For this reason, in our calculations we employ two types of basis
sets: singly-augmented aug-cc-pVTZ and aug(nH)-cc-pVTZ and doubly-augmented
d-aug-cc-pVDZ. Due to the limitations of the available unrestricted MP3
implementations, we do not report VEA and AEA of AT computed at the
MP3/aug(nH)-cc-pVTZ level of theory. We compare our results to VEA’s and AEA’s
calculated using the EOM-CCSD method with the aug-cc-pV$X$Z+5s5p4d ($X$ = D,
T) basis sets from Ref. 117.

As expected, the computed VEA’s show a strong dependence on the basis set. For
both GC and AT base pairs, the MP$n$ and ADC($n$) VEA’s computed using the
d-aug-cc-pVDZ basis set are in a very good agreement with the reference
results, while the results obtained using the singly-augmented triple-zeta
basis sets show much larger deviations, indicating that the additional diffuse
orbitals are critical for accurate predictions of VEA’s of these molecules.
The best results are demonstrated by the EA-ADC(3)/d-aug-cc-pVDZ method, which
predicts VEA’s of 0.09/$-$0.03 eV for the GC/AT pair, in a close agreement
with the reference VEA’s of 0.10/0.00 eV. The EA-ADC(2) method combined with
the d-aug-cc-pVDZ basis set shows somewhat larger deviations from the
reference results with VEA’s of 0.14/0.02 eV. The MP$n$/d-aug-cc-pVDZ ($n$ =
2, 3) methods exhibit larger errors with VEA’s of $\sim$ 0.06/$-$0.07 eV.

The computed AEA’s exhibit much weaker basis set dependence. Except for EA-
ADC(2)/aug-cc-pVTZ, all levels of theory with either of the two basis sets
predict a positive AEA for GC and a negative AEA for AT, in agreement with the
reference results. For GC, the EA-ADC(3) method shows the best agreement with
the reference AEA (0.18 to 0.24 eV) predicting AEA’s of 0.15 and 0.13 eV
computed using the aug(nH)-cc-pVTZ and d-aug-cc-pVDZ basis sets, respectively.
When combined with the aug(nH)-cc-pVTZ basis set, the EA-ADC(3) method yields
an accurate AEA of AT ($-$0.14 eV), which agrees well with the reference AEA’s
in the $-0.10$ to $-0.01$ range. At the EA-ADC(3)/d-aug-cc-pVDZ level of
theory, the computed AEA has a lower value of $-$0.38 eV. When using the
d-aug-cc-pVDZ basis set, the MP$n$ methods predict AEA’s that are close to
those of EA-ADC($n$), but further away from the reference. The fact that the
VEA’s and AEA’s computed using EA-ADC($n$) and MP$n$ with the d-aug-cc-pVDZ
basis set agree well with the reference results for GC, but show larger
deviations for AT suggests that the AEA of the latter base pair requires a
higher-level description of electron correlation effects, beyond those
included in the finite-order perturbation theories.

##  VI Conclusions

In this work, we presented an efficient implementation of the single-reference
algebraic diagrammatic construction (ADC) theory for simulating electron
attachment (EA) and ionization (IP) energies and spectra of molecules (EA/IP-
ADC($n$), $n=2,3$). Our new implementation, available in the PySCF program,
consists of the spin-restricted (RADC) and unrestricted (UADC) codes for
closed- and open-shell systems, respectively. The RADC and UADC programs are
capable of efficient memory and disk management by employing the density
fitting approximation for the two-electron integrals and take advantage of the
highly optimized subroutines for tensor contractions and parallelization.
Combining the RADC implementation with density fitting results in a four-fold
speedup of the EA/IP-ADC(3) calculations for closed-shell molecules, relative
to the UADC implementation with unapproximated two-electron integrals. Our
numerical tests show that the density fitting approximation has a negligible
effect on the accuracy of the ADC methods but dramatically lowers the
computational cost of the integral transformation and the input/output
operation count.

We demonstrated the capabilities of our efficient implementation by applying
the EA/IP-ADC($n$) ($n=2,3$) methods to the TEMPO radical and the EA-ADC($n$)
methods to the guanine–cytosine (GC) and adenine–thymine (AT) DNA base pairs.
Our EA- and IP-ADC(3) computations of closed-shell molecules with up to 1028
basis functions and open-shell molecules with up to 758 molecular orbitals
are, to the best of our knowledge, the largest calculations using these
methods reported to date. For the TEMPO radical, the photoelectron spectra
simulated using IP-ADC($n$) ($n=2,3$) are in a good agreement with an
experimental spectrum. The EA-ADC($n$) ($n=2,3$) methods predict that the
TEMPO radical has a negative vertical electron affinity (VEA), suggesting that
the electron attachment is energetically unfavorable at the equilibrium
geometry of the neutral radical. Upon structural relaxation, the electron
attachment becomes favorable, as indicated by the small positive values (0.08
to 0.49 eV) of the adiabatic electron affinities (AEA) computed using EA-
ADC($n$) ($n=2,3$). For the DNA base pairs, our EA-ADC results confirm the
dipole-bound nature of the vertically-attached electronic states reported in
the earlier EOM-CCSD study.Tripathi:2019p10131  The VEA’s computed using the
EA-ADC methods show a strong dependence on the number of diffuse basis
functions employed in the calculation. Using the doubly-augmented d-aug-cc-
pVDZ basis set, the EA-ADC($n$) ($n=2,3$) VEA’s are in a good agreement with
the available reference data. The computed AEA’s show a weaker basis set
dependence. For GC, both EA-ADC(2) and EA-ADC(3) predict accurate AEA’s
($\sim$ 0.14 eV), which agree well with AEA from EOM-CCSD (0.18 to 0.24 eV).
The computed AEA of AT show larger errors ($\sim$ 0.15 to 0.30 eV),
highlighting the importance of infinite-order electron correlation effects for
the electron affinity of this system.

In summary, our work demonstrates that EA/IP-ADC($n$) ($n=2,3$) are efficient
and accurate theoretical methods for simulations of charged excitations and
spectra that can be routinely applied to closed- and open-shell molecules with
more than 1000 basis functions. Our efficient implementation enables new
applications that were previously out of reach and can be combined with
periodic boundary conditions for simulations of crystalline chemical systems.
Work along this direction is ongoing in our group.

##  VII Supplementary Material

See supplementary material for the optimized Cartesian geometries of the TEMPO
radical, DNA base pairs, tables showing the density fitting approximation
errors for EA/IP-DF-ADC($n$) ($n$ = 2, 3), and tabulated photoelectron spectra
of the TEMPO radical computed using IP-UADC($n$) ($n$ = 2, 3).

##  VIII Acknowledgements

This work was supported by the start-up funds provided by the Ohio State
University. Additionally, S.B. was supported by a fellowship from the
Molecular Sciences Software Institute under NSF Grant ACI-1547580.
Computations were performed at the Ohio Supercomputer Center under projects
PAS1583.OhioSupercomputerCenter1987

##  IX Data availability

The data that supports the findings of this study are available within the
article and its supplementary material. Additional data can be made available
upon reasonable request.

##  X Appendix: Calculation of spectroscopic factors in the RADC
implementation

Here, we demonstrate how to calculate the spectroscopic factors in the RADC
implementation. As discussed in Section III.2, solution of the RADC eigenvalue
problem (15) delivers the spin-adapted eigenvectors $\mathbf{Y_{\pm}}$. Since
the RADC eigenvectors contain fewer elements than those in the UADC
implementation, they need to be renormalized to ensure that the computed
spectroscopic factors are consistent with those calculated using UADC. For IP-
RADC, the elements of the renormalized eigenvectors $\mathbf{\tilde{Y}}_{-}$
can be obtained from $\mathbf{Y_{-}}$ as follows:

$\displaystyle\tilde{Y}_{-i,\mu}$
$\displaystyle=\frac{Y_{-i,\mu}}{\sqrt{Y^{norm}_{-\mu}}}$ (46)
$\displaystyle\tilde{Y}_{-ija,\mu}$
$\displaystyle=\frac{Y_{-ija,\mu}}{\sqrt{Y^{norm}_{-\mu}}}$ (47)

where $Y_{-\mu}^{norm}$ is a norm of root $\mu$, defined as

$\displaystyle Y_{-\mu}^{norm}$ $\displaystyle=\sum_{i}Y_{-i,\mu}Y_{-i,\mu}$
$\displaystyle+\sum_{ija}(2Y_{-ija,\mu}Y_{-ija,\mu}-Y_{-ija,\mu}Y_{-jia,\mu})$
(48)

Once computed, the renormalized eigenvectors $\mathbf{\tilde{Y}}_{-}$ are
contracted with the elements of the spin-adapted transition moments matrix
$\mathbf{{T}}_{-}$ to obtain the spin-adapted spectroscopic amplitudes
($\mathbf{{X}}_{-}$)

$\displaystyle X_{-p\mu}$ $\displaystyle=\sum_{i}T_{-p,i}\tilde{Y}_{-i,\mu}$
$\displaystyle+\sum_{ija}(2T_{-p,ija}\tilde{Y}_{-ija,\mu}-T_{-p,jia}\tilde{Y}_{-ija,\mu})$
(49)

For EA-RADC, expressions for the elements of the renormalized eigenvectors
$\mathbf{\tilde{Y}}_{+}$ and spectroscopic amplitudes $\mathbf{{X}}_{+}$ can
be obtained from Eqs. 46, 47, 48 and 49 by replacing the 1h ($i$) and 2h-1p
($ija$) excitation labels with those of the 1p ($a$) and 2p-1h ($abi$)
excitations. Using the elements of $\mathbf{{X}}_{\pm}$, the spectroscopic
factors can be obtained as follows

$\displaystyle P_{\pm\mu}$ $\displaystyle=2\sum_{p}|X_{\pm p\mu}|^{2}$ (50)

## References

  * (1) L. S. Cederbaum, J. Phys. B: At. Mol. Phys. 8, 290 (1974).
  * (2) W. Von Niessen, J. Schirmer, and L. S. Cederbaum, Comput. Phys. Rep. 1, 57 (1984).
  * (3) J. V. Ortiz, WIREs Comput. Mol. Sci. 3, 123 (2012).
  * (4) S. Hirata, M. R. Hermes, J. Simons, and J. V. Ortiz, J. Chem. Theory Comput. 11, 1595 (2015).
  * (5) M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. 44, 55 (1992).
  * (6) M. Nooijen and J. G. Snijders, Int. J. Quantum Chem. 48, 15 (1993).
  * (7) M. Nooijen and J. G. Snijders, J. Chem. Phys. 102, 1681 (1995).
  * (8) J. F. Stanton and R. J. Bartlett, J. Chem. Phys. 98, 7029 (1993).
  * (9) A. I. Krylov, Annu. Rev. Phys. Chem. 59 (2008).
  * (10) K. Kowalski, K. Bhaskaran-Nair, and W. A. Shelton, J. Chem. Phys. 141, 094102 (2014).
  * (11) K. Bhaskaran-Nair, K. Kowalski, and W. A. Shelton, J. Chem. Phys. 144, 144101 (2016).
  * (12) B. Peng and K. Kowalski, J. Chem. Theory Comput. 14, 4335 (2018).
  * (13) T. Fleig, D. Edvardsson, S. T. Banks, and J. H. Eland, Chem. Phys. 343, 270 (2008).
  * (14) L. Streit, F. B. Machado, and R. Custodio, Chem. Phys. Lett. 506, 22 (2011).
  * (15) S. McKechnie, G. H. Booth, A. J. Cohen, and J. M. Cole, J. Chem. Phys. 142, 194114 (2015).
  * (16) H.-V. Nguyen, T. A. Pham, D. Rocca, and G. Galli, Phys. Rev. B 85, 081101 (2012).
  * (17) L. Hedin, Phys. Rev. A 139, 796 (1965).
  * (18) S. V. Faleev, M. van Schilfgaarde, and T. Kotani, Phys. Rev. Lett. 93, 126406 (2004).
  * (19) M. van Schilfgaarde, T. Kotani, and S. Faleev, Phys. Rev. Lett. 96, 226402 (2006).
  * (20) J. B. Neaton, M. S. Hybertsen, and S. G. Louie, Phys. Rev. Lett. 97, 216405 (2006).
  * (21) G. Samsonidze, M. Jain, J. Deslippe, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 107, 186404 (2011).
  * (22) M. J. van Setten, F. Weigend, and F. Evers, J. Chem. Theory Comput. 9, 232 (2013).
  * (23) L. Reining, WIREs Comput. Mol. Sci. 8, e1344 (2017).
  * (24) F. Pavošević, C. Peng, J. V. Ortiz, and E. F. Valeev, J. Chem. Phys. 147, 121101 (2017).
  * (25) M. F. Lange and T. C. Berkelbach, J. Chem. Theory Comput. 14, 4224 (2018).
  * (26) O. J. Backhouse and G. H. Booth, J. Chem. Theory Comput. 16, 6294 (2020).
  * (27) A. L. Dempwolff, M. Schneider, M. Hodecker, and A. Dreuw, J. Chem. Phys. 150, 064108 (2019).
  * (28) A. L. Dempwolff, A. C. Paul, A. M. Belogolova, A. B. Trofimov, and A. Dreuw, J. Chem. Phys. 152, 024113 (2020).
  * (29) A. L. Dempwolff, A. C. Paul, A. M. Belogolova, A. B. Trofimov, and A. Dreuw, J. Chem. Phys. 152, 024125 (2020).
  * (30) J. Schirmer, Phys. Rev. A 26, 2395 (1982).
  * (31) J. Schirmer, L. S. Cederbaum, and O. Walter, Phys. Rev. A 28, 1237 (1983).
  * (32) J. Schirmer, Phys. Rev. A 43, 4647 (1991).
  * (33) F. Mertins and J. Schirmer, Phys. Rev. A 53, 2140 (1996).
  * (34) J. Schirmer and A. B. Trofimov, J. Chem. Phys. 120, 11449 (2004).
  * (35) A. Dreuw and M. Wormit, WIREs Comput. Mol. Sci. 5, 82 (2014).
  * (36) J. H. Starcke, M. Wormit, and A. Dreuw, J. Chem. Phys. 130, 024104 (2009).
  * (37) P. H. P. Harbach, M. Wormit, and A. Dreuw, J. Chem. Phys. 141, 064113 (2014).
  * (38) A. Barth and J. Schirmer, J. Phys. B: At. Mol. Phys. 18, 867 (1985).
  * (39) J. Wenzel, M. Wormit, and A. Dreuw, J. Comput. Chem. 35, 1900 (2014).
  * (40) J. Schirmer, A. B. Trofimov, and G. Stelter, J. Chem. Phys. 109, 4734 (1998).
  * (41) A. B. Trofimov and J. Schirmer, J. Chem. Phys. 123, 144115 (2005).
  * (42) G. Angonoa, O. Walter, and J. Schirmer, J. Chem. Phys. 87, 6789 (1987).
  * (43) J. Schirmer and A. Thiel, J. Chem. Phys. 115, 10621 (2001).
  * (44) A. Thiel, J. Schirmer, and H. Köppel, J. Chem. Phys. 119, 2088 (2003).
  * (45) S. Knippenberg, D. R. Rehn, M. Wormit, J. H. Starcke, I. L. Rusakova, A. B. Trofimov, and A. Dreuw, J. Chem. Phys. 136, 064107 (2012).
  * (46) A. Y. Sokolov, J. Chem. Phys. 149, 204113 (2018).
  * (47) K. Chatterjee and A. Y. Sokolov, J. Chem. Theory Comput. 15, 5908 (2019).
  * (48) K. Chatterjee and A. Y. Sokolov, J. Chem. Theory Comput. 16, 6343 (2020).
  * (49) M. Schneider, D. Y. Soshnikov, D. M. P. Holland, I. Powis, E. Antonsson, M. Patanen, C. Nicolas, C. Miron, M. Wormit, A. Dreuw, and A. B. Trofimov, J. Chem. Phys. 143, 144103 (2015).
  * (50) Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng, X. Feng, D. Ghosh, M. Goldey, P. R. Horn, L. D. Jacobson, I. Kaliman, R. Z. Khaliullin, T. Kuś, A. Landau, J. Liu, E. I. Proynov, Y. M. Rhee, R. M. Richard, M. A. Rohrdanz, R. P. Steele, E. J. Sundstrom, H. L. Woodcock III, P. M. Zimmerman, D. Zuev, B. Albrecht, E. Alguire, B. Austin, G. J. O. Beran, Y. A. Bernard, E. Berquist, K. Brandhorst, K. B. Bravaya, S. T. Brown, D. Casanova, C.-M. Chang, Y. Chen, S. H. Chien, K. D. Closser, D. L. Crittenden, M. Diedenhofen, R. A. Distasio JR., H. Do, A. D. Dutoi, R. G. Edgar, S. Fatehi, L. Fusti-Molnar, A. Ghysels, A. Golubeva-Zadorozhnaya, J. Gomes, M. W. D. Hanson-Heine, P. H. P. Harbach, A. W. Hauser, E. G. Hohenstein, Z. C. Holden, T.-C. Jagau, H. Ji, B. Kaduk, K. Khistyaev, J. Kim, J. Kim, R. A. King, P. Klunzinger, D. Kosenkov, T. Kowalczyk, C. M. Krauter, K. U. Lao, A. D. Laurent, K. V. Lawler, S. V. Levchenko, C. Y. Lin, F. Liu, E. Livshits, R. C. Lochan, A. Luenser, P. Manohar, S. F. Manzer, S.-P. Mao, N. Mardirossian, A. V. Marenich, S. A. Maurer, N. J. Mayhall, E. Neuscamman, C. M. Oana, R. Olivares-Amaya, D. P. O’Neill, J. A. Parkhill, T. M. Perrine, R. Peverati, A. Prociuk, D. R. Rehn, E. Rosta, N. J. Russ, S. M. Sharada, S. Sharma, D. W. Small, A. Sodt, T. Stein, D. Stück, Y.-C. Su, A. J. W. Thom, T. Tsuchimochi, V. Vanovschi, L. Vogt, O. Vydrov, T. Wang, M. A. Watson, J. Wenzel, A. F. White, C. F. Williams, J. Yang, S. Yeganeh, S. R. Yost, Z.-Q. You, I. Y. Zhang, X. Zhang, Y. Zhao, B. R. Brooks, G. K.-L. Chan, D. M. Chipman, C. J. Cramer, W. A. Goddard III, M. S. Gordon, W. J. Hehre, A. Klamt, H. F. Schaefer, M. W. Schmidt, C. D. Sherrill, D. G. Truhlar, A. Warshel, X. Xu, A. Aspuru-Guzik, R. Baer, A. T. Bell, N. A. Besley, J.-D. Chai, A. Dreuw, B. D. Dunietz, T. R. Furlani, S. R. Gwaltney, C.-P. Hsu, Y. Jung, J. Kong, D. S. Lambrecht, W. Liang, C. Ochsenfeld, V. A. Rassolov, L. V. Slipchenko, J. E. Subotnik, T. Van Voorhis, J. M. Herbert, A. I. Krylov, P. M. W. Gill, and M. Head-Gordon, Mol. Phys. 113, 184 (2014).
  * (51) A. B. Trofimov and J. Schirmer, in Proceedings of 14th European Symposium on Gas Phase Electron Diffraction (2011), p. 77.
  * (52) M. Schneider, Weiterentwicklung und Implementierung quantenchemischer Methoden zur direkten Berechnung von Ionisationspotentialen und Elektroaffinitäten (Heidelberg, 2015), zs.fassung in engl. und dt. Sprache.
  * (53) S. Banerjee and A. Y. Sokolov, J. Chem. Phys. 151, 224112 (2019).
  * (54) J. Liu, C. Hättig, and S. Höfener, J. Chem. Phys. 152, 174109 (2020).
  * (55) Q. Sun, X. Zhang, S. Banerjee, P. Bao, M. Barbry, N. S. Blunt, N. A. Bogdanov, G. H. Booth, J. Chen, Z.-H. Cui, J. J. Eriksen, Y. Gao, S. Guo, J. Hermann, M. R. Hermes, K. Koh, P. Koval, S. Lehtola, Z. Li, J. Liu, N. Mardirossian, J. D. McClain, M. Motta, B. Mussard, H. Q. Pham, A. Pulkin, W. Purwanto, P. J. Robinson, E. Ronca, E. R. Sayfutyarova, M. Scheurer, H. F. Schurkus, J. E. T. Smith, C. Sun, S.-N. Sun, S. Upadhyay, L. K. Wagner, X. Wang, A. F. White, J. D. Whitfield, M. J. Williamson, S. Wouters, J. Yang, J. M. Yu, T. Zhu, T. C. Berkelbach, S. Sharma, A. Y. Sokolov, and G. K.-L. Chan, J. Chem. Phys. 153, 024109 (2020).
  * (56) J. L. Whitten, J. Chem. Phys. 58, 4496 (1973).
  * (57) O. Vahtras, J. Almlöf, and M. W. Feyereisen, Chem. Phys. Lett. 213, 514 (1993).
  * (58) M. Feyereisen, G. Fitzgerald, and A. Komornicki, Chem. Phys. Lett. 208, 359 (1993).
  * (59) A. P. Rendell and T. J. Lee, J. Chem. Phys. 101, 400 (1994).
  * (60) F. Weigend, Phys. Chem. Chem. Phys. 4, 4285 (2002).
  * (61) D. H. Kobe, Ann. Phys. 19, 448 (1962).
  * (62) D. Mukherjee and W. Kutzelnigg, in Many-Body Methods in Quantum Chemistry (Springer, Berlin, Heidelberg, Berlin, Heidelberg, 1989), pp. 257–274.
  * (63) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R’ıo, M. Wiebe, P. Peterson, P. G’erard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357 (2020).
  * (64) E. R. Davidson, J. Comput. Phys. 17, 87 (1975).
  * (65) B. Liu, The Simultaneous Expansion Method for the Iterative Solution of Several of the Lowest-Lying Eigenvalues and Corresponding Eigenvectors of Large Real-Symmetric Matrices, Technical report, Berkeley, CA, USA (1978).
  * (66) A. B. Trofimov, G. Stelter, and J. Schirmer, J. Chem. Phys. 117, 6402 (2002).
  * (67) M. Nooijen and R. J. Bartlett, J. Chem. Phys. 102, 3629 (1995).
  * (68) B. I. Dunlap, J. W. D. Connolly, and J. R. Sabin, J. Chem. Phys. 71, 3396 (1979).
  * (69) H.-J. Werner, F. R. Manby, and P. J. Knowles, J. Chem. Phys. 118, 8149 (2003).
  * (70) U. Bozkaya, J. Chem. Phys. 141, 124108 (2014).
  * (71) E. Epifanovsky, D. Zuev, X. Feng, K. Khistyaev, Y. Shao, and A. I. Krylov, J. Chem. Phys. 139, 134105 (2013).
  * (72) X. Wang, A. Y. Sokolov, J. M. Turney, and H. F. Schaefer, J. Chem. Theory Comput. 12, 4833 (2016).
  * (73) B. Helmich and C. Hättig, Comput. Theor. Chem 1040-1041, 35 (2014).
  * (74) D. Mester, P. R. Nagy, and M. Kállay, J. Chem. Phys. 148, 094111 (2018).
  * (75) M. F. Herbst, M. Scheurer, T. Fransson, D. R. Rehn, and A. Dreuw, WIREs Comput. Mol. Sci. p. e1462 (2020).
  * (76) T. H. Dunning Jr, J. Chem. Phys. 90, 1007 (1989).
  * (77) R. A. Kendall, T. H. Dunning Jr, and R. J. Harrison, J. Chem. Phys. 96, 6796 (1992).
  * (78) D. E. Woon and T. H. Dunning Jr, J. Chem. Phys. 100, 2975 (1994).
  * (79) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 (1989).
  * (80) R. J. Bartlett, J. Watts, S. Kucharski, and J. Noga, Chem. Phys. Lett. 165, 513 (1990).
  * (81) C. Møller and M. S. Plesset, Phys. Rev. 46, 618 (1934).
  * (82) D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrinceIII, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford, and C. D. Sherrill, J. Chem. Phys. 152, 184108 (2020).
  * (83) F. Weigend, A. Köhn, and C. Hättig, J. Chem. Phys. 116, 3175 (2002).
  * (84) C. Hättig, Phys. Chem. Chem. Phys. 7, 59 (2005).
  * (85) C. Melania Oana and A. I. Krylov, J. Chem. Phys. 127, 234106 (2007).
  * (86) A. Rohatgi, Webplotdigitizer: Version 4.3 (2020).
  * (87) U. Bozkaya, J. Chem. Theory Comput. 10, 2371 (2014).
  * (88) U. Bozkaya, J. Chem. Theory Comput. 12, 1179 (2016).
  * (89) C. Hättig and F. Weigend, J. Chem. Phys. 113, 5154 (2000).
  * (90) S. Barriga, Synlett 2001, 0563 (2001).
  * (91) C. Galli, Z. Rappoport, and J. F. Liebman, Nitroxyl radicals in The Chemistry of Hydroxylamines, Oximes and Hydroxamic Acids, Volume 1, vol. 175 (John Wiley & Sons, 2008).
  * (92) G. I. Likhtenshtein, J. Yamauchi, S. Nakatsuji, A. Smirnov, and R. Tamura, Nitroxides in Physicochemistry (Wiley Online Library, 2008), chapter 8, pp. 269–301.
  * (93) K. Nakahara, K. Oyaizu, and H. Nishide, Chem. Lett. 40, 222 (2011).
  * (94) T. Janoschka, N. Martin, U. Martin, C. Friebe, S. Morgenstern, H. Hiller, M. D. Hager, and U. S. Schubert, Nature 527, 78 (2015).
  * (95) K. Nakahara, S. Iwasa, J. Iriyama, Y. Morioka, M. Suguro, M. Satoh, and E. J. Cairns, Electrochim. Acta 52, 921 (2006).
  * (96) J. P. Blinco, J. L. Hodgson, B. J. Morrow, J. R. Walker, G. D. Will, M. L. Coote, and S. E. Bottle, J. Org. Chem. 73, 6763 (2008).
  * (97) F. Kato, N. Hayashi, T. Murakami, C. Okumura, K. Oyaizu, and H. Nishide, Chem. Lett. 39, 464 (2010).
  * (98) C. Buhrmester, L. Moshurchak, R. Wang, and J. Dahn, J. Electrochem. Soc. 153, A1800 (2006).
  * (99) I. Novak, L. J. Harrison, B. Kovač, and L. M. Pratt, J. Org. Chem. 69, 7628 (2004).
  * (100) B. Kovač, I. Ljubić, A. Kivimäki, M. Coreno, and I. Novak, Phys. Chem. Chem. Phys. 16, 10734 (2014).
  * (101) I. Ljubić, J. Chem. Theory Comput. 10, 2333 (2014).
  * (102) I. Ljubić, A. Kivimäki, and M. Coreno, Phys. Chem. Chem. Phys. 18, 10207 (2016).
  * (103) D. Kubala, K. Regeta, R. Janečková, J. Fedor, S. Grimme, A. Hansen, P. Nesvadba, and M. Allan, Mol. Phys. 111, 2033 (2013).
  * (104) C. H. Sohn, S. Yin, I. Peng, J. A. Loo, and J. Beauchamp, Int. J. Mass Spectrom. 390, 49 (2015).
  * (105) T. Gunasekara, G. P. Abramo, A. Hansen, H. Neugebauer, M. Bursch, S. Grimme, and J. R. Norton, J. Am. Chem. Soc. 141, 1882 (2019).
  * (106) R. Teoule, Int. J. Radiat. Biol. Relat. Stud. Phys. Chem. Med. 51, 573 (1987).
  * (107) R. Barrios, P. Skurski, and J. Simons, J. Phys. Chem. B 106, 7991 (2002).
  * (108) B. Boudaıffa, P. Cloutier, D. Hunting, M. A. Huels, and L. Sanche, Science 287, 1658 (2000).
  * (109) K. Aflatooni, G. Gallup, and P. Burrow, J. Phys. Chem. A 102, 6205 (1998).
  * (110) A. O. Colson, B. Besler, and M. D. Sevilla, J. Phys. Chem. 96, 9787 (1992).
  * (111) N. A. Richardson, S. S. Wesolowski, and H. F. Schaefer, J. Am. Chem. Soc. 124, 10163 (2002).
  * (112) N. A. Richardson, S. S. Wesolowski, and H. F. Schaefer, J. Phys. Chem. B 107, 848 (2003).
  * (113) J. Reynisson and S. Steenken, Phys. Chem. Chem. Phys. 4, 5353 (2002).
  * (114) I. Al-Jihad, J. Smets, and L. Adamowicz, J. Phys. Chem. A 104, 2994 (2000).
  * (115) J. Smets, A. F. Jalbout, and L. Adamowicz, Chem. Phys. Lett. 342, 342 (2001).
  * (116) B. Peng, S. R. McNew, Q.-s. Li, Y. Xie, and H. F. Schaefer, Chem. Phys. Lett. 523, 120 (2012).
  * (117) D. Tripathi and A. K. Dutta, J. Phys. Chem. A 123, 10131 (2019).
  * (118) Ohio Supercomputer Center, Technical report (1987).