We present an efficient implementation of the second- and third-order single-reference algebraic diagrammatic construction (ADC) theory for electron attachment and ionization 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, and 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 (2,2,6,6-tetramethylpiperidin-1-yl)oxyl (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 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 (EAs), ionization potentials (IPs), bandgaps, and photoelectron spectra.^{1–4} 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 (CC)^{5–12} theory and multireference configuration interaction,^{13,14} 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,^{15–23} 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.^{24–29}

An attractive alternative to traditional coupled cluster theory for simulating charged electronic states is algebraic diagrammatic construction (ADC) theory^{30–35} 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/Vis^{30,36,37} and x-ray absorption,^{38,39} photoelectron,^{27–29,31,40–44} and two-photon absorption spectra.^{45} The ADC methods for charged excitations (EA/IP), first proposed in 1983, were based on the Dyson framework that involves the perturbative self-energy expansion of the Dyson equation and couples the EA and IP components of the one-electron propagator.^{31} Later work by Schirmer, Trofimov, and co-workers^{40,41} 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.^{46–48}

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 the ionization potentials (IPs). 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* ≤ 3) were reported by Schneider *et al.*^{49} and Dempwolff *et al.,*^{27–29} and the latter is available in the Q-CHEM package.^{50} Implementations of the non-Dyson EA-ADC(*n*) (*n* ≤ 3) methods were limited to two early benchmark studies,^{51,52} a more extensive benchmark of EA-ADC(*n*) (*n* ≤ 3) for closed- and open-shell systems performed in our group,^{53} and a recent efficient implementation of EA-ADC(2) by Liu *et al.*^{54} 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,^{27,28} the largest non-Dyson EA-ADC(3) calculations reported to date were limited to 298 orbitals.^{53}

In this work, we present an efficient implementation of the non-Dyson EA- and IP-ADC(*n*) methods (*n* ≤ 3), now available as a module in the PySCF software package.^{55} 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 approximation^{56–60} 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 (Sec. II), followed by the overview of our implementation (Sec. III), where we discuss the general algorithm (Sec. III A) along with the details of spin adaptation (Sec. III B) and the density fitting approximation (Sec. III C). We next present computational details in Sec. IV and benchmark the accuracy of the density fitting approximation along with the efficiency of our implementation in Sec. V A. Finally, we present our results for the TEMPO radical in Sec. V B and two DNA base pairs (guanine–cytosine and adenine–thymine) in Sec. V C. We present our conclusions in Sec. VI.

## II. THEORY

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

We first give a brief overview of the non-Dyson algebraic diagrammatic construction (ADC) theory^{30–35,53} for electron attachment and removal, where the central object of interest is the retarded single-particle Green’s function,

Here, the forward ($Gpq+(\omega )$) and backward ($Gpq\u2212(\omega )$) components of the propagator describe electron attachment (EA) and ionization potential (IP) processes, respectively, $|\Psi 0N\u2009$ and $E0N$ are the eigenfunction and eigenvalue of the electronic Hamiltonian *H* of the *N*-electron system, and a frequency *ω* ≡ *ω*′ + *iη* is expressed in terms of its real (*ω*′) and infinitesimal imaginary (*iη*) components. The operators $ap\u2020$ 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,

where $\Omega \u0303\xb1$ are the diagonal matrices of the exact vertical attachment ($\omega n=EnN+1\u2212E0N$) and ionization ($\omega n=E0N\u2212EnN\u22121$) energies and $X\u0303\xb1$ are the matrices of the spectroscopic amplitudes with elements $X\u0303+pn=\u2009\Psi 0N|ap|\Psi nN+1\u2009$ and $X\u0303\u2212qn=\u2009\Psi 0N|aq\u2020|\Psi nN\u22121\u2009$, where $|\Psi nN+1\u2009$ and $|\Psi nN\u22121\u2009$ are the exact eigenstates and $EnN+1$ and $EnN\u22121$ are the exact eigenvalues of the Hamiltonian *H* for the (*N* + 1)- and (*N* − 1)-electron system, respectively. Equation (2) is known as the spectral (or Lehmann) representation of the propagator.^{61}

The propagator in Eqs. (1) and (2) expressed in the basis of the exact eigenstates $|\Psi nN+1\u2009$ and $|\Psi nN\u22121\u2009$ is computationally expensive. In non-Dyson ADC, the propagator is rewritten in an approximate (in general, non-orthogonal) basis of the (*N* ± 1)-electron states,

where the **M**_{±} matrices are no longer diagonal, in contrast to $\Omega \u0303\xb1$ in Eq. (2). In Eq. (3), **M**_{±} and **T**_{±}, 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 **S**_{±} matrices. For all single-reference ADC methods discussed in this work, the approximate basis is orthonormal, and thus, we will assume **S**_{±} = **1** henceforth.

Starting with Eq. (3), the approximate propagator **G**_{±}(*ω*) is expanded in a perturbative series,

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_{±}(*ω*) are derived by evaluating **M**_{±} and **T**_{±} up to the *n*th order in perturbation theory using the intermediate state representation approach^{32–34} or via the formalism of the effective Liouvillean theory.^{46,53,62}

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,

where |Φ⟩ 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*, … and *a*, *b*, *c*, *d*, … indices, respectively,

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 us to determine contributions to **M**_{±} and **T**_{±} at each order in perturbation theory.

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

where operators $H\u0303(k)$ and $\xe3p(k)$ are the *k*th-order contributions to the effective Hamiltonian $H\u0303=e\u2212(T\u2212T\u2020)He(T\u2212T\u2020)$ and observable $\xe3p=e\u2212(T\u2212T\u2020)ape(T\u2212T\u2020)$ operators, while [$\cdots \u2009$] and [$\cdots \u2009$]_{+} denote the commutator and anticommutator, respectively. The operators $h+\mu (k)\u2020$ in Eqs. (7) and (8) are used to construct a set of basis states $|\Psi +\mu (k)\u2009=h+\mu (k)\u2020|\Phi \u2009$ 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

where the ionization operators $h\u2212\mu (k)\u2020$ define the (*N* − 1)-electron basis states $|\Psi \u2212\mu (k)\u2009=h\u2212\mu (k)\u2020|\Phi \u2009$. For the low-order EA/IP-ADC(*n*) approximations (*n* ≤ 3), only the zeroth- and first-order components of $h\xb1\mu (k)\u2020$ are required and are given as follows: $h+\mu (0)\u2020=aa\u2020$, $h+\mu (1)\u2020=ab\u2020aa\u2020ai$, $h\u2212\mu (0)\u2020=ai$, and $h\u2212\mu (1)\u2020=aa\u2020ajai$. The $H\u0303(k)$ and $\xe3p(k)$ operators in Eqs. (7)–(10) are obtained by expanding $H\u0303$ and $\xe3p$ using the Baker–Campbell–Hausdorff (BCH) formula,

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* ≤ 3) approximations require calculating up to the *n*th order single-excitation amplitudes ($tia(k)$) and up to the (*n* − 1)th order double-excitation amplitudes ($tijab(k\u22121)$). At each order *k*, these amplitudes are computed by solving a system of projected amplitude equations,

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 **M**_{±} and **T**_{±} are obtained, the EA/IP-ADC(*n*) transition energies are computed by solving the eigenvalue problem,

where **Ω**_{±} is a diagonal matrix of eigenvalues and **Y**_{±} is a matrix of eigenvectors that are used to compute the spectroscopic amplitudes,

which provide information about transition intensities. The elements of **X**_{±} can be used to calculate the spectroscopic factors,

and the ADC propagator and density of states,

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

### A. 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.^{55} 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 given as follows:

Transform two-electron repulsion integrals (

*pq*|*rs*) from the atomic to the RHF/UHF molecular orbital basis.Compute the amplitudes of the effective Hamiltonian ($tia(k)$ and $tijab(k\u22121)$,

*k*≤*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.Compute small blocks of the effective Hamiltonian matrix

**M**_{±}with respect to the zeroth-order operators ($h\xb1\mu (0)\u2020$, Sec. II A). These are the 1p–1p [*M*_{+ab}, Eq. (7)] and 1h–1h [*M*_{−ij}, Eq. (9)] blocks of**M**_{+}and**M**_{−}for EA- and IP-ADC(*n*), respectively.Define a function for calculating a matrix-vector product

*σ*_{±}=**M**_{±}**Y**_{±}using the approximate eigenvectors**Y**_{±}and the small blocks of the effective Hamiltonian matrix precomputed in step 3. Solve the eigenvalue problem (15) by optimizing the eigenvectors**Y**_{±}until convergence using an iterative diagonalization algorithm to obtain low-energy EA/IPs.Using the converged eigenvectors

**Y**_{±}, compute spectroscopic amplitudes**X**_{±}[Eq. (16)] and spectroscopic factors*P*_{±μ}[Eq. (17)]. If desired, compute spectral density*A*(*ω*) [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 (BLASs) 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.^{63} 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 a disk, and (iii) *density-fitted* (DF) algorithm that significantly reduces disk usage and enables memory savings by approximating two-electron integrals, as discussed in Sec. III C. We use the multiroot Davidson algorithm^{64,65} available in PySCF for iterative solution of the ADC eigenvalue problem.

Table I 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 $O(O2V3)$ and $O(O2V4)$, respectively. Importantly, precomputing the 1p–1p and 1h–1h blocks of the effective Hamiltonian matrix **M**_{±} at step 3 and reusing them in the Davidson diagonalization at step 4 lowers the cost of computing the *σ*_{±} = **M**_{±}**Y**_{±} matrix-vector products from $O(O2V4)$ to $O(O2V3)$ and $O(O3V2)$ 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 $O(O2V4)$ contractions, respectively, independent of the number of iterations in the Davidson algorithm. A similar approach that allows us 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).^{27}

. | EA . | IP . | Relative prefactor . | |||
---|---|---|---|---|---|---|

Component . | ADC(2) . | ADC(3) . | ADC(2) . | ADC(3) . | RADC(3) . | UADC(3) . |

Amplitudes and MPn energy | O^{2}V^{3} | O^{2}V^{4} | O^{2}V^{3} | O^{2}V^{4} | 1 | 3 |

M_{ab}(1p–1p) or M_{ij}(1h–1h) | O^{2}V^{3} | O^{2}V^{4} | O^{3}V^{2} | O^{2}V^{4} | 1 | 3 |

$\sigma \xb1(i+1)=M\xb1Y\xb1(i)$ | OV^{3} | O^{2}V^{3} | O^{3}V | O^{3}V^{2} | 1 | 4 |

Spectroscopic factors (P_{±μ}) | 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 |

. | EA . | IP . | Relative prefactor . | |||
---|---|---|---|---|---|---|

Component . | ADC(2) . | ADC(3) . | ADC(2) . | ADC(3) . | RADC(3) . | UADC(3) . |

Amplitudes and MPn energy | O^{2}V^{3} | O^{2}V^{4} | O^{2}V^{3} | O^{2}V^{4} | 1 | 3 |

M_{ab}(1p–1p) or M_{ij}(1h–1h) | O^{2}V^{3} | O^{2}V^{4} | O^{3}V^{2} | O^{2}V^{4} | 1 | 3 |

$\sigma \xb1(i+1)=M\xb1Y\xb1(i)$ | OV^{3} | O^{2}V^{3} | O^{3}V | O^{3}V^{2} | 1 | 4 |

Spectroscopic factors (P_{±μ}) | 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 |

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 **M**_{±} at step 3, the contributions of these blocks to the matrix-vector products (i.e., *∑*_{b}*M*_{+ab}*Y*_{b} and *∑*_{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 $O(O2V2)$ for EA/IP-ADC(2) and to $O(O2V3)$/$O(OV4)$ 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 (Sec. III B) and give an overview of the density fitting approximation used for the two-electron integrals (Sec. III C). The efficiency of our implementation and the accuracy of the density fitting approximation are analyzed in Sec. V A.

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

In our earlier work,^{53} we presented a general spin-orbital derivation and implementation of EA/IP-ADC(*n*) (*n* ≤ 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 ionization^{44} and neutral^{38,66} 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 (**M**_{±}) and effective transition moments (**T**_{±}) of matrix elements in Eqs. (7)–(10) by integrating over the spin variables. We denote each spin-orbital with an index *p*_{σ}, where *σ* = {*α*, *β*} 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*, … and *a*, *b*, *c*, *d*, … indices to label the occupied and virtual MOs, respectively, and reserve *p*, *q*, *r*, *s*, … to denote a general (occupied or virtual) MO.

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

where our notation for the **M**_{−} and **T**_{−} matrix elements matches that of Ref. 53.

To obtain working equations for **M**_{−} and **T**_{−} in IP-RADC, we write their spin-adapted matrix elements in terms of the spin-orbital Fock matrix eigenvalues ($\epsilon p\sigma $), antisymmetrized two-electron integrals ($\u2009p\sigma 1q\sigma 2\Vert r\sigma 3s\sigma 4\u2009$), and *k*th-order amplitudes of the effective Hamiltonian [$ti\sigma 1a\sigma 2(k)$ and $ti\sigma 1j\sigma 2a\sigma 3b\sigma 4(k)$, Eq. (6)] and perform spin-integration as follows:

where (*pq*|*rs*) denotes a chemist’s 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,

and iteratively optimizing the spin-adapted eigenvectors with elements **Y**_{−} = {*Y*_{−i}, *Y*_{−ija}}. The converged eigenvectors are then used to evaluate the spectroscopic factors, as described in the Appendix. Expressions for **M**_{+}, **T**_{+}, and *σ*_{+} 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 Sec. III A, with a threefold reduction in the formal computational prefactor, as shown in Table I. Most significant savings are achieved for the solution of the EA/IP-ADC eigenvalue problem, where the matrix-vector products *σ*_{±} are expressed in terms of a much smaller number of the spin-adapted **M**_{±} matrix elements [e.g., *M*_{−i,j}, *M*_{−i,jka}, and *M*_{−ilb,jka} in Eqs. (20)–(24)] compared to that of the spin-orbital implementation. We will analyze the computational efficiency of our EA/IP-RADC implementation in Sec. V A.

### C. Density fitting in EA/IP-ADC

In addition to using the conventional two-electron integrals [(*pq*|*rs*) or $\u2009p\sigma 1q\sigma 2|r\sigma 3s\sigma 4\u2009$], our EA/IP-RADC and UADC implementations feature the support of the density fitting (DF) approximation^{56–58,68} 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

where the $bpqQ$ coefficients are defined in terms of the MOs (*ϕ*_{p}) and the auxiliary basis functions (*χ*_{P}) as

The use of the DF approximation has virtually no effect on the results of the EA/IP-ADC calculations (see Sec. V A 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 the disk, and (ii) it dramatically lowers the computational scaling of the two-electron integral transformation from $O(N5)$ to $O(NauxN3)$ (where *N* and *N*_{aux} are the numbers of *ϕ*_{p} and *χ*_{P}, respectively), providing significant computational savings for large molecules and basis sets. In our implementation, we store the $bpqQ$ tensors either in memory or on the disk and recompute sectors of (*pq*|*rs*) on demand. Similar to other methods that employ the DF approximation,^{69–72} 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.^{54,73–75}

## IV. COMPUTATIONAL DETAILS

All EA/IP-ADC(*n*) calculations (*n* ≤ 3) were performed using the ADC module in PySCF^{53,55} 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.^{76–78} Specifically, in Sec. V A, 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 (Sec. V B) and DNA base pairs (Sec. V C) 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 Sec. V A were optimized using coupled cluster theory with single, double, and perturbative triple excitations [CCSD(T)].^{79,80} The EA/IP-ADC(*n*) calculations (*n* ≤ 3) of the TEMPO radical and anion used structures optimized using the *n*th-order Møller–Plesset perturbation theory^{81} (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 program,^{82} and the optimized structures are reported in the supplementary material. 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,^{60,83,84} where *X* corresponded to that of the main one-electron basis set.

Throughout this article, 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 (VEAs) were obtained directly from our EA-ADC(*n*) calculations, the corresponding adiabatic electron electron affinities (AEA) were determined as follows:

where $\Delta EMPn=EMPnN+1$(neutral) $\u2212\u2009EMPnN+1$(anion) is the difference of the MP*n* energies of the anionic system computed at the neutral and anion optimized geometries. The TEMPO radical AEAs have been corrected to account for the zero-point vibrational energy (ZPVE),

where *ZPVE*(*N*) and *ZPVE*(*N* + 1) are the ZPVEs 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,^{85} which are computed from the overlap of the ground state wavefunction and the electron-attached/ionized states,

where |*ϕ*_{p}⟩ are the reference Hartree–Fock orbitals and *X*_{±pμ} are the EA/IP-ADC spectroscopic amplitudes defined in Eq. (16). All Dyson orbitals were computed using the aug-cc-pVDZ basis set.

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

for a range of frequencies *ω*, where *ω*_{μ} denotes the ionization energies obtained by solving Eq. (15), *η* is a small broadening, and *P*_{−μ} 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.^{86}

## V. RESULTS

### A. 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 potential (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, the results are presented for sets of five closed-shell (H_{2}O, HF, N_{2}, SiH_{2}, H_{2}CO) and five open-shell (OH, NH, O_{2}, NH_{2}, CH_{3}) molecules computed using the RADC and UADC implementations, respectively (with data for individual molecules available in Tables S1 and S2 of the supplementary material). 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 (MAEs) of ∼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 MAEs less than 0.0005 eV. Incorporating the DF approximation into the SCF reference calculations increases the computed MAEs 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.^{71,72,87–89}

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 Sec. III A and is presented in Fig. 2.

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 3. Significant savings are observed for all steps of the ADC computation, with ∼2.5-fold speedup for calculating the amplitudes of the effective Hamiltonian and the 1p–1p/1h–1h blocks of the **M**_{±} matrix and more than a fourfold reduction in time for calculating the *σ*_{±} vector and spectroscopic factors, in good agreement with theoretical prefactors reported in Table I. 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 us to achieve nearly a fourfold speedup for both EA- and IP-ADC(3), pushing the limits of their applicability to larger chemical systems and basis sets. In Secs. V B and V C, we demonstrate this by applying EA/IP-RADC(*n*) (*n* = 2, 3) to the TEMPO radical and the DNA base pairs.

### B. TEMPO radical

We first study the (2,2,6,6-tetramethylpiperidin-1-yl)oxyl radical, commonly known as TEMPO (Fig. 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.^{90–92} TEMPO and its derivatives are also used in solar cells, lithium batteries, and radical sensors.^{93–98}

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.^{99–102} These studies demonstrated 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.^{103} 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 the 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.^{104} 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 an AEA of 0.44 eV.^{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.*^{103} Figure 4 shows the density of occupied states (DOSs) of TEMPO computed using IP-UADC(*n*) (*n* = 2, 3) with a broadening of 0.15 eV. For each method, the 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 IPs due to its higher-level description of electron correlation effects. The shifted DOS computed using both IP-UADC(2) and IP-UADC(3) agrees well with the experimental photoelectron spectrum, accurately reproducing peak spacings and relative peak intensities, with IP-UADC(3) showing better agreement than IP-ADC(2) for ionization energies greater than 10 eV. Upon shifting, both methods predict a signal at ∼7.2 eV, corresponding to the *X*^{2}*A*′ → 1^{1}*A*′ photoelectron transition, and two peaks at ∼9.0 eV and 9.3 eV that can be assigned as the *X*^{2}*A*′ → 1^{3}*A*″ and *X*^{2}*A*′ → 1^{1}*A*″ 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 (Fig. 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*′ → 1^{3}*A*″ and *X*^{2}*A*′ → 1^{1}*A*″ transitions to have similar spectral density, IP-UADC(3) yields a smaller spectral density for the higher-energy *X*^{2}*A*′ → 1^{1}*A*″ transition, in good agreement with the experiment. The IP-UADC(3) DOS also shows a weak satellite transition at ∼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 it in Table II. In agreement with a previous DFT study,^{103} all MP*n* and EA-UADC(*n*) methods predict a negative VEA ranging between −0.72 eV [EA-UADC(3)] and −0.21 eV (MP2), indicating an endothermic vertical electron attachment. The VEA computed using MP*n* and EA-UADC(*n*) is close to each other at each order in perturbation theory (within 0.15 eV) and becomes more negative as the order increases from *n* = 2–3.

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.3^{a} | 0.64^{b}, 0.44^{c} | 0.14^{a} |

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.3^{a} | 0.64^{b}, 0.44^{c} | 0.14^{a} |

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),^{105} 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 good agreement with experimental and theoretical findings of Kubala *et al.*^{103} Table II also reports AEA corrected by the zero-point vibrational energies (ZPVE) of the neutral and anionic states (see Sec. 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.

Figures 6(a) and 6(b) 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 [Fig. 6(a)]. Once the TEMPO geometry is relaxed to that of its anion, the lowest-energy electron-attached state acquires the valence-bound character [Fig. 6(b)] where most of the orbital is localized on the NO group of the radical.

### C. DNA 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 Fig. 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.^{106–108} In our previous work,^{53} we reported VEAs 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 VEAs of the nucleobases decreased in the order G > U ≳ C ≳ T > A, in good agreement with the G > U ≳ T ≳ C > A trend from the experiment.^{109} 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.^{110,111}

Although no experimental electron affinities have been reported for the isolated base pairs, several theoretical studies have been performed.^{110–117} Early computational investigations using Hartree–Fock and density functional theories (DFTs)^{110–112} 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 EAs. Another DFT study reported A to be the primary site for the electron localization in the AT anion.^{113} Electron attachment to GC and AT was studied by the Adamowicz^{114,115} and Schaefer^{116} groups using the MP2 method with the double-zeta basis sets incorporating additional diffuse functions. They reported AEAs 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.^{117} Using the aug-cc-pVTZ basis set with additional 5s5p4d diffuse functions, they reported the computed VEAs for the GC/AT pair to be 0.099/0.003 eV and the AEAs 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.

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.*,^{117} 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 [Fig. 8(b)] and T of AT [Fig. 9(b)], which agrees well with the results of other post-Hartree–Fock calculations reported in the literature.^{114–117}

Table III 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 VEAs 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 VEAs and AEAs calculated using the EOM-CCSD method with the aug-cc-pV*X*Z+5s5p4d (*X* = D, T) basis sets from Ref. 117.

. | Guanine–cytosine . | Adenine–thymine . | ||
---|---|---|---|---|

Method . | 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.14^{a} |

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+5s5p4d^{b} | 0.10 | 0.18 | 0.00 | −0.10 |

EOM-CCSD/aug-cc-pVTZ+5s5p4d^{b} | 0.10 | 0.24 | 0.00 | −0.01 |

. | Guanine–cytosine . | Adenine–thymine . | ||
---|---|---|---|---|

Method . | 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.14^{a} |

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+5s5p4d^{b} | 0.10 | 0.18 | 0.00 | −0.10 |

EOM-CCSD/aug-cc-pVTZ+5s5p4d^{b} | 0.10 | 0.24 | 0.00 | −0.01 |

As expected, the computed VEAs show a strong dependence on the basis set. For both GC and AT base pairs, the MP*n* and ADC(*n*) VEAs computed using the d-aug-cc-pVDZ basis set are in 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 VEAs of these molecules. The best results are demonstrated by the EA-ADC(3)/d-aug-cc-pVDZ method, which predicts VEAs of 0.09/−0.03 eV for the GC/AT pair, in close agreement with the reference VEAs 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 VEAs of 0.14/0.02 eV. The MP*n*/d-aug-cc-pVDZ (*n* = 2, 3) methods exhibit larger errors with VEAs of ∼ 0.06/−0.07 eV.

The computed AEAs 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 best agreement with the reference AEA (0.18 eV–0.24 eV), predicting AEAs of 0.15 eV 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 AEAs 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 AEAs that are close to those of EA-ADC(*n*) but further away from the reference. The fact that the VEAs and AEAs 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 potential (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 fourfold 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 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 eV–0.49 eV) of the adiabatic electron affinities (AEAs) 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.^{117} The VEAs 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) VEAs are in good agreement with the available reference data. The computed AEAs show a weaker basis set dependence. For GC, both EA-ADC(2) and EA-ADC(3) predict accurate AEAs (∼ 0.14 eV), which agree well with AEA from EOM-CCSD (0.18 eV–0.24 eV). The computed AEA of AT shows larger errors (∼ 0.15 eV–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.

## SUPPLEMENTARY MATERIAL

See the 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).

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material. Additional data that support the findings of this study are available from the corresponding author upon reasonable request.

### 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 Sec. III B, the solution of the RADC eigenvalue problem (15) delivers the spin-adapted eigenvectors **Y**_{±}. 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 $Y\u0303\u2212$ can be obtained from **Y**_{−} as follows:

where $Y\u2212\mu norm$ is a norm of root *μ*, defined as

Once computed, the renormalized eigenvectors $Y\u0303\u2212$ are contracted with the elements of the spin-adapted transition moment matrix **T**_{−} to obtain the spin-adapted spectroscopic amplitudes (**X**_{−}),

For EA-RADC, expressions for the elements of the renormalized eigenvectors $Y\u0303+$ and spectroscopic amplitudes **X**_{+} can be obtained from Eqs. (A1)–(A4) 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 **X**_{±}, the spectroscopic factors can be obtained as follows: