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.

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) theory30–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/Vis30,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-workers40,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 approximation56–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.

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

Gpq(ω)=Gpq+(ω)+Gpq(ω)=Ψ0N|ap(ωH+E0N)1aq|Ψ0N+Ψ0N|aq(ω+HE0N)1ap|Ψ0N.
(1)

Here, the forward (Gpq+(ω)) and backward (Gpq(ω)) components of the propagator describe electron attachment (EA) and ionization potential (IP) processes, respectively, |Ψ0N and E0N are the eigenfunction and eigenvalue of the electronic Hamiltonian H of the N-electron system, and a frequency ωω′ + is expressed in terms of its real (ω′) and infinitesimal imaginary () components. The operators ap and ap are the usual fermionic creation and annihilation operators, respectively.

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

G±(ω)=X̃±(ω1Ω̃±)1X̃±,
(2)

where Ω̃± are the diagonal matrices of the exact vertical attachment (ωn=EnN+1E0N) and ionization (ωn=E0NEnN1) energies and X̃± are the matrices of the spectroscopic amplitudes with elements X̃+pn=Ψ0N|ap|ΨnN+1 and X̃qn=Ψ0N|aq|ΨnN1, where |ΨnN+1 and |ΨnN1 are the exact eigenstates and EnN+1 and EnN1 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 |ΨnN+1 and |ΨnN1 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,

G±(ω)=T±(ωS±M±)1T±,
(3)

where the M± matrices are no longer diagonal, in contrast to Ω̃± 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,

G±(ω)=G±(0)(ω)+G±(1)(ω)++G±(n)(ω)+.
(4)

Truncating this expansion at the nth order defines the propagator of the nth 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 nth order in perturbation theory using the intermediate state representation approach32–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,

|Ψ0N=eTT|Φ,
(5)

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,

T=mNTm,Tm=1(m!)2ijabtijabaaabajai.
(6)

To define perturbative series (4), the electronic Hamiltonian H is partitioned into the zeroth-order part H(0) and a perturbation V = HH(0). This allows us to determine contributions to M± and T± at each order in perturbation theory.

The nth-order contributions to the EA-ADC matrices have the following form:

M+μν(n)=klmk+l+m=nΦ|[h+μ(k),[H̃(l),h+ν(m)]]+|Φ,
(7)
   T+pν(n)=klk+l=nΦ|[ãp(k),h+ν(l)]+|Φ,
(8)

where operators H̃(k) and ãp(k) are the kth-order contributions to the effective Hamiltonian H̃=e(TT)He(TT) and observable ãp=e(TT)ape(TT) operators, while [] and []+ denote the commutator and anticommutator, respectively. The operators h+μ(k) in Eqs. (7) and (8) are used to construct a set of basis states |Ψ+μ(k)=h+μ(k)|Φ necessary for expanding the eigenstates of the (N + 1)-electron system in the EA-ADC equations. Similarly, the nth-order contributions to the IP-ADC matrices are given as

Mμν(n)=klmk+l+m=nΦ|[hμ(k),[H̃(l),hν(m)]]+|Φ,
(9)
   Tpν(n)=klk+l=nΦ|[ãp(k),hν(l)]+|Φ,
(10)

where the ionization operators hμ(k) define the (N − 1)-electron basis states |Ψμ(k)=hμ(k)|Φ. For the low-order EA/IP-ADC(n) approximations (n ≤ 3), only the zeroth- and first-order components of h±μ(k) are required and are given as follows: h+μ(0)=aa, h+μ(1)=abaaai, hμ(0)=ai, and hμ(1)=aaajai. The H̃(k) and ãp(k) operators in Eqs. (7)(10) are obtained by expanding H̃ and ãp using the Baker–Campbell–Hausdorff (BCH) formula,

H̃=H(0)+V+[H(0),T(1)T(1)]+[H(0),T(2)T(2)]+12![V+(V+[H(0),T(1)T(1)]),T(1)T(1)]+,
(11)
ãp=ap+[ap,T(1)T(1)]+[ap,T(2)T(2)]+12![[ap,T(1)T(1)],T(1)T(1)]+,
(12)

and collecting terms at the kth 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 nth order single-excitation amplitudes (tia(k)) and up to the (n − 1)th order double-excitation amplitudes (tijab(k1)). At each order k, these amplitudes are computed by solving a system of projected amplitude equations,

Φ|aiaaH̃(k)|Φ=0,
(13)
Φ|aiajabaaH̃(k)|Φ=0.
(14)

The resulting amplitudes are equivalent to those that define the kth-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,

M±Y±=Y±Ω±,
(15)

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

X±=T±Y±,
(16)

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

P±μ=p|X±pμ|2,
(17)

and the ADC propagator and density of states,

G±(ω)=X±ω1Ω±1X±,
(18)
A(ω)=1πImTr G±(ω).
(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.

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:

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

  2. Compute the amplitudes of the effective Hamiltonian (tia(k) and tijab(k1), kn) by solving Eqs. (13) and (14), and evaluate the nth-order Møller–Plesset (MPn) correlation energy, where n is the order of the ADC(n) approximation.

  3. Compute small blocks of the effective Hamiltonian matrix M± with respect to the zeroth-order operators (h±μ(0), Sec. II A). These are the 1p–1p [M+ab, Eq. (7)] and 1h–1h [Mij, Eq. (9)] blocks of M+ and M for EA- and IP-ADC(n), respectively.

  4. 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.

  5. 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 algorithm64,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 

TABLE I.

Comparison of formal computational scaling of the EA/IP-ADC(n) methods (n ≤ 3) for each step of the algorithm outlined in Sec. III A. 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).

EAIPRelative prefactor
ComponentADC(2)ADC(3)ADC(2)ADC(3)RADC(3)UADC(3)
Amplitudes and MPn energy O2V3 O2V4 O2V3 O2V4 
Mab(1p–1p) or Mij(1h–1h) O2V3 O2V4 O3V2 O2V4 
σ±(i+1)=M±Y±(i) OV3 O2V3 O3V O3V2 
Spectroscopic factors (P±μO2V3 O2V3 O3V2 O3V2 
Overall O2V3 O2V4 O2V3 O2V4 
EAIPRelative prefactor
ComponentADC(2)ADC(3)ADC(2)ADC(3)RADC(3)UADC(3)
Amplitudes and MPn energy O2V3 O2V4 O2V3 O2V4 
Mab(1p–1p) or Mij(1h–1h) O2V3 O2V4 O3V2 O2V4 
σ±(i+1)=M±Y±(i) OV3 O2V3 O3V O3V2 
Spectroscopic factors (P±μO2V3 O2V3 O3V2 O3V2 
Overall O2V3 O2V4 O2V3 O2V4 

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., bM+abYb and jMijYj) 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.

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 ionization44 and neutral38,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

Miα,jα=Miβ,jβ=Mi,j,
(20)
Miα,jβkαaβ=Miβ,jαkβaα=Mi,jka,
(21)
Miα,jαkαaα=Miβ,jβkβaβ=Mi,jkaMi,kja,
(22)
Miβlαbβ,jβkαaβ=Miαlβbα,jαkβaα=Milb,jka,
(23)
Miαlαbα,jαkαaα=Miβlβbβ,jβkβaβ=Milb,jkaMlib,jkaMilb,kja+Mlib,kja,
(24)
Tpα,iα=Tpβ,iβ=Tp,i,
(25)
Tpα,jβkαaβ=Tpβ,jαkβaα=Tp,jka,
(26)
Tpα,jαkαaα=Tpβ,jβkβaβ=Tp,jkaTp,kja,
(27)

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 (εpσ), antisymmetrized two-electron integrals (pσ1qσ2rσ3sσ4), and kth-order amplitudes of the effective Hamiltonian [tiσ1aσ2(k) and tiσ1jσ2aσ3bσ4(k), Eq. (6)] and perform spin-integration as follows:

εpα=εpβ=εp,
(28)
pαqβrαsβ=pβqαrβsα=(pr|qs),
(29)
pαqβrβsα=pβqαrαsβ=(ps|qr),
(30)
pαqαrαsα=pβqβrβsβ=(pr|qs)(ps|qr),
(31)
tiαaα(k)=tiβaβ(k)=tia(k),
(32)
tiαjβaαbβ(k)=tiβjαaβbα(k)=tijab(k),
(33)
tiαjβaβbα(k)=tiβjαaαbβ(k)=tijba(k),
(34)
tiαjαaαbα(k)=tiβjβaβbβ(k)=tijab(k)tijba(k),
(35)

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,

σi=jMi,jYj+2jkaMi,jkaYjkajkaMi,jkaYkja,
(36)
σjka=iMi,jkaYi+ilbMjka,ilbYilb,
(37)

and iteratively optimizing the spin-adapted eigenvectors with elements Y = {Yi, Yija}. 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., Mi,j, Mi,jka, and Milb,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.

In addition to using the conventional two-electron integrals [(pq|rs) or pσ1qσ2|rσ3sσ4], our EA/IP-RADC and UADC implementations feature the support of the density fitting (DF) approximation56–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

(pq|rs)QNauxbpqQbrsQ,
(38)

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

bpqQ=PNaux(pq|P)J12PQ,
(39)
(pq|P)=ϕp(1)ϕq(1)1r12χP(2)dr1dr2,
(40)
JPQ=χP(1)1r12χQ(2)dr1dr2.
(41)

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 Naux 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

All EA/IP-ADC(n) calculations (n ≤ 3) were performed using the ADC module in PySCF53,55 and used the singly or doubly augmented Dunning’s correlation consistent basis sets denoted as aug-cc-pVXZ or d-aug-cc-pVXZ, 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 nth-order Møller–Plesset perturbation theory81 (MPn) 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 MPn calculations employed the density fitting approximation. The correlation MPn energies and ADC(n) excitation energies were computed using the resolution-of-identity (aug-cc-pVXZ-RI) auxiliary basis sets, while the SCF energies and molecular orbitals were obtained using the aug-cc-pVXZ-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 = ENEN+1), whereas a positive ionization energy corresponds to an endothermic process (IP = EN−1EN). 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:

AEA=VEA+ΔEMPn,
(42)

where ΔEMPn=EMPnN+1(neutral) EMPnN+1(anion) is the difference of the MPn 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),

AEA(ZPVE)=AEA+ZPVE(N)ZPVE(N+1),
(43)

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,

ϕ±μDyson=pX±pμ|ϕp,
(44)

where |ϕp⟩ are the reference Hartree–Fock orbitals and X± 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

A(ω)=1πImμPμωωμ+iη
(45)

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 

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 (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 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.

FIG. 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 of the supplementary material for data on individual molecules.

FIG. 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 of the supplementary material for data on individual molecules.

Close modal

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.

FIG. 2.

Comparison of the computational wall time (in s) 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 are 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 Sec. III A is shown with a different color. All calculations were performed on a single computer with a 12-core Intel Xeon Gold 6136 3.00 GHz CPU.

FIG. 2.

Comparison of the computational wall time (in s) 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 are 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 Sec. III A is shown with a different color. All calculations were performed on a single computer with a 12-core Intel Xeon Gold 6136 3.00 GHz CPU.

Close modal

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.

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 

FIG. 3.

Molecular structure of the TEMPO radical.

FIG. 3.

Molecular structure of the TEMPO radical.

Close modal

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.103Figure 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 X2A′ → 11A′ photoelectron transition, and two peaks at ∼9.0 eV and 9.3 eV that can be assigned as the X2A′ → 13A″ and X2A′ → 11A″ 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 X2A′ → 13A″ and X2A′ → 11A″ transitions to have similar spectral density, IP-UADC(3) yields a smaller spectral density for the higher-energy X2A′ → 11A″ 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.

FIG. 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.

FIG. 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.

Close modal
FIG. 5.

Dyson orbitals for the three lowest-energy ionization transitions of the neutral TEMPO radical corresponding to the 11A′ (a), 13A″ (b), and 11A″ (c) states of the TEMPO cation computed using the IP-UADC(2) method with the aug-cc-pVTZ basis set.

FIG. 5.

Dyson orbitals for the three lowest-energy ionization transitions of the neutral TEMPO radical corresponding to the 11A′ (a), 13A″ (b), and 11A″ (c) states of the TEMPO cation computed using the IP-UADC(2) method with the aug-cc-pVTZ basis set.

Close modal

We now turn our attention to VEA and AEA of the TEMPO radical computed using the MPn and EA-UADC(n) methods (n = 2, 3) and presented it in Table II. In agreement with a previous DFT study,103 all MPn 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 MPn 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.

TABLE II.

Vertical (VEA) and adiabatic (AEA) electron attachment energies (in eV) of the TEMPO radical computed using MPn 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.

MethodVEAAEAAEA (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 
MethodVEAAEAAEA (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.

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 MPn 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.103Table 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.

FIG. 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.

FIG. 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.

Close modal

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

FIG. 7.

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

FIG. 7.

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

Close modal

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 Adamowicz114,115 and Schaefer116 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 

FIG. 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.

FIG. 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.

Close modal
FIG. 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.

FIG. 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.

Close modal

Table III presents the vertical (VEA) and adiabatic (AEA) electron affinities of GC and AT computed using the second- and third-order MPn 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-pVXZ+5s5p4d (X = D, T) basis sets from Ref. 117.

TABLE III.

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

Guanine–cytosineAdenine–thymine
MethodVEAAEAVEAAEA
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 
Guanine–cytosineAdenine–thymine
MethodVEAAEAVEAAEA
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 ΔEMPn was computed using MP2/aug-cc-pVTZ.

b

Reference 117.

As expected, the computed VEAs show a strong dependence on the basis set. For both GC and AT base pairs, the MPn 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 MPn/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 MPn 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 MPn 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.

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.

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

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 

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.

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 Ỹ can be obtained from Y as follows:

Ỹi,μ=Yi,μYμnorm,
(A1)
Ỹija,μ=Yija,μYμnorm,
(A2)

where Yμnorm is a norm of root μ, defined as

Yμnorm=iYi,μYi,μ+ija(2Yija,μYija,μYija,μYjia,μ).
(A3)

Once computed, the renormalized eigenvectors Ỹ are contracted with the elements of the spin-adapted transition moment matrix T to obtain the spin-adapted spectroscopic amplitudes (X),

Xpμ=iTp,iỸi,μ+ija(2Tp,ijaỸija,μTp,jiaỸija,μ).
(A4)

For EA-RADC, expressions for the elements of the renormalized eigenvectors Ỹ+ 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:

P±μ=2p|X±pμ|2.
(A5)
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
,
Wiley Interdiscip. Rev.: 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
).
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. D.
Eland
,
Chem. Phys.
343
,
270
(
2008
).
14.
L.
Streit
,
F. B. C.
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.
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
,
Wiley Interdiscip. Rev.: 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
).
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
,
Wiley Interdiscip. Rev.: 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.
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
(
Moscow State University
,
2011
), p.
77
.
52.
M.
Schneider
, “
Weiterentwicklung und Implementierung quantenchemischer Methoden zur direkten Berechnung von Ionisationspotentialen und Elektroaffinitäten
,” Ph.D. thesis, Heidelberg University,
2015
.
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.
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
).
62.
D.
Mukherjee
and
W.
Kutzelnigg
,
Many-Body Methods in Quantum Chemistry
(
Springer
,
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érard-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 No. LBL-8158,
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
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
10
,
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. D.
Watts
,
S. A.
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.
DePrince
 III
,
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
.
91.
C.
Galli
,
Z.
Rappoport
, and
J. F.
Liebman
,
Nitroxyl Radicals in the Chemistry of Hydroxylamines, Oximes and Hydroxamic Acids, Volume 1
(
John Wiley & Sons
,
2008
), Vol. 175.
92.
G. I.
Likhtenshtein
,
J.
Yamauchi
,
S.
Nakatsuji
,
A.
Smirnov
, and
R.
Tamura
,
Nitroxides in Physicochemistry
(
Wiley Online Library
,
2008
), Chap. 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. M.
Moshurchak
,
R. L.
Wang
, and
J. R.
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. L.
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. A.
Gallup
, and
P. D.
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.
See http://osc.edu/ark:/19495/f5s1ph73 for Ohio Supercomputer Center.

Supplementary Material