The ionization potential (IP) of a molecule quantifies the energy required to remove an electron from the system. As such, it is a fundamental quantity in the context of redox chemistry, charge transfer, and molecular electronics. The accurate theoretical prediction of this property is therefore highly desirable for virtual materials design. Furthermore, vertical IPs are of interest in the development of many-body Green’s function methods like the GW formalism, as well as density functionals and semiempirical methods. In this contribution, we report over 1468 vertical valence IPs calculated with the IP variant of equation-of-motion coupled cluster theory with singles and doubles (IP-EOM-CCSD) covering 155 molecules. The purpose of this is two-fold: First, the quality of the predicted IPs is compared with respect to experiments and higher-order coupled cluster theory. This confirms the overall high accuracy and robustness of this method, with some outliers which are discussed in detail. Second, a large set of consistent theoretical reference values for vertical valence IPs are generated. This addresses a lack of reliable reference data for lower-lying valence IPs, where experimental data are often unavailable or of dubious quality. The benchmark set is then used to assess the quality of the eigenvalues predicted by different density functional approximations (via Bartlett’s IP-eigenvalue theorem) and the extended Koopmans’ theorem approach. The QTP family of functionals are found to be remarkably accurate, low-cost alternatives to IP-EOM-CCSD.

Ionization potentials (IPs)96 and electron affinities (EAs) have long played a leading role in chemical thinking, e.g., in such central concepts as frontier molecular orbitals, hardness, and electronegativity.1–5 In early semiempirical quantum chemistry approaches like the (extended) Hückel6,7 and Pariser–Parr–Pople methods,8,9 the use of IPs to determine molecular integrals10 was critical to their success. IPs are also an important molecular signature for the analysis of chemical compounds through photoelectron spectroscopy.

With the rise of molecular electronics, the prediction of these properties has furthermore garnered great interest in applied research.11,12 IPs and EAs are the fundamental observables that govern the driving forces in organic (opto-)electronic devices like transistors, solar cells, and LEDs.13–17 Because the chemical space of possible molecules for these applications is far too great to be fully explored experimentally, virtual screening of potential candidates is a highly attractive prospect.18,19 Hence, the trustworthy computation of these properties is of great interest.

Arguably the most useful, accurate technique to this end is equation-of-motion coupled cluster theory (EOM-CC),20–27 specifically its IP and EA-EOM variants.28–30 However, the high computational demand of this method (in particular owing to the required ground-state CC calculation) so far precludes its use in truly large-scale applications, such as screening thousands of molecules. In principle, a much more efficient approximation for vertical IPs (VIPs) can be obtained from the orbital eigenvalues of mean-field methods like Hartree-Fock (HF) or Kohn-Sham density functional theory (KS-DFT) (through Bartlett’s IP-eigenvalue theorem).31 Alternatively, electron-propagator and Green’s function based methods including the popular GW method offer alternatives of intermediate computational cost (and accuracy).32–34 Interestingly, since both CC and GW are rooted in many-body theory, they are in fact intimately connected as recently discussed by Lange and Berkelbach.35 

DFT methods can be quite accurate for principal IPs (i.e., those corresponding to the highest occupied molecular orbital, HOMO), when the quantity is calculated as a difference between the total energies of the neutral and cationic systems (the ∆SCF method). This is, however, not generally applicable to lower-lying states as converging the (excited) cation KS calculation is not trivial in many cases. Furthermore, the ∆SCF method requires a separate KS calculation for each state of interest, making the method impractical for the computation of spectra (or the large number of states considered in this paper).

Meanwhile, if KS orbital eigenvalues are used, accurate predictions for specific systems are usually the result of a fair amount of calibration with respect to how the functional is defined (i.e., the underlying density functional approximation, the amount of HF exchange, and other effects like the choice of range-separation parameters).36–38 This also reflects on GW calculations, where the DFT starting point strongly influences the results in the most commonly used non-self-consistent implementation (G0W0).32,38,39

Beyond predicting VIPs, there are also methodological concerns regarding the accuracy of KS orbital eigenvalues.31,37,40,41 It was shown by Bartlett et al. that to be able to consistently treat ground and excited states (e.g., via adiabatic time-dependent DFT), the KS orbital energies must necessarily correspond to vertical ionization energies.41,42 In this paper, we further insist on the direct determination of IPs as eigenvalues in KS-DFT instead of functional differences. We consider this a fundamental property in its own right.43 Interestingly, virtually no popular functional complies with this IP theorem.40 Meanwhile, functionals that were purposefully designed to comply with it (like those from the QTP family44–46) show promising performance in cases where DFT routinely fails, e.g., for Rydberg excitations and self-interaction error-prone cases.44,45,47 We have also recently worked on translating these concepts into semiempirical methods beyond DFT.48 In a similar vein, the localized orbital scaling correction of Yang and co-workers appears to improve both delocalization errors and VIPs.49 

For both application and method development, the validation of computed VIPs is therefore crucial. As always, the experiment is in principle the ideal benchmark. However (particularly for lower-lying valence states), experimental references are often not available. Even when experimental values are known, the assignment of a spectroscopic peak to the correct vertical IP can be ambiguous as vibronic coupling and zero-point contributions cannot be disentangled from the purely “electronic” VIP that is obtained from theory (as the difference between cation and neutral energies at the minimum of the ground state potential energy surface).50–53 The comparison between experimental VIPs and electronic structure calculations should thus be taken with a grain of salt.

Recently, HOMO IPs for 100 molecules predicted with the GW method were reported.54 However, this was explicitly to compare the numerical precision of different implementations, not to provide benchmark values with respect to experiments. There is also a benchmark set of the HOMO VIPs and VEAs of 24 molecules [using ∆CCSD(T) reference energies] which has been used to test a large range of electronic structure methods.55,56 However, to the best of our knowledge, no comprehensive collection of highly accurate theoretical VIPs (including lower lying valence states) has so far been reported.

In this contribution, this deficiency is addressed by presenting a dataset of 1468 valence VIPs of 155 organic molecules at the IP variant of equation-of-motion coupled cluster theory with singles and doubles (IP-EOM-CCSD) level. The accuracy of this “workhorse” method is established by comparison with higher level methods including connected triple and quadruple excitations (IP-EOM-CCSDT and IP-EOM-CCSDTQ). Where available, experimental reference values are also provided. Finally, the accuracy of computationally efficient approximations based on (extended) Koopmans’ theorem is also tested.

The EOM-CC of Bartlett26 method is a time-independent generalization of the CC linear response theory of Monkhorst, and others.20,23,27,57 As it is commonly used today, EOM-CC was first formulated and implemented by Stanton and for the description of excited states.20 The IP and EA variants allow describing systems with more or less electrons than the reference state and are closely related to the Fock-space CC methods.29,58–61 For completeness, we will provide a short overview on (IP)EOM-CC; further details can be found in the original references.

EOM-CC calculations start with a ground-state CC calculation, where the wavefunction is defined as

Ψ0=eT^ϕ0,
(1)

with a single-determinant reference ϕ0 and the cluster operator T^. The latter is defined via a set of t-amplitudes (th1hmp1pm) and normal-ordered second-quantized creation (pm) and annihilation (hm) operators as

T^=T^1++T^n,
(2)
T^m=1m!2p1pmvirt.h1hmocc.th1hmp1pmp1pmhmh1.
(3)

By truncating T^ at double (m = 2), triple (m = 3), or quadruple (m = 4) excitations, one obtains specific CC methods, abbreviated as CCSD, CCSDT, and CCSDTQ, respectively.

In EOM-CC, excited or ionized states are obtained using a linear CI-like operator ,

Ψx=^eT^ϕ0,
(4)
^=^1++^m,
(5)
^m=1m!n!p1pnvirt.h1hmocc.rh1hmp1pnp1pnhmh1.
(6)

The form of the and T operators is thus very similar, with the difference that the number of creation and annihilation operators is not necessarily balanced in (hence the separate indices n and m). If n = m, charge neutral excitations are obtained (i.e., the original EOM-CC, or EE-EOM-CC framework). Alternatively, for n = m − 1, the IP-EOM-CC equations are obtained and ^ is a net annihilation operator. As before, the level of truncation defines the specific methods. In this paper, we apply truncations at the two-holes-one-particle (2h1p), three-holes-two-particles (3h2p), or four-holes-three-particles (4h3p) levels (abbreviated as IP-EOM-CCSD, IP-EOM-CCSDT, and IP-EOM-CCSDTQ respectively).

By inserting eq. (4) into the Schrödinger equation and rearranging, one obtains

eT^HeT^^ϕ0=ωϕ0.
(7)

This can be rewritten as an eigenvalue equation with the similarity transformed Hamiltonian H¯,

H¯=eT^HeT^.
(8)

The final steps of an EOM-CC calculation thus superficially resemble the configuration interaction (CI) approach. Specifically, the H¯ matrix is constructed in a basis of excited determinants relative to ϕ0 and solved for the desired number of eigenvalues. Two important differences to CI are, however, that (1) higher excitations are included in H¯ via the nonlinear operator eT^ and (2) H¯ is not hermitian since eT^ is not unitary. As a consequence of the latter, the final states have a biorthogonal form with different left- and right-hand eigenvectors (𝔏 and ) but a single set of eigenvalues. It is therefore necessary to obtain both 𝔏 and when calculating properties and energy derivatives.

Canonically, the excitation level of the and T operators is chosen to be the same. However, for practical purposes there is no violation of fundamental concepts if one chose to define the ground state effective Hamiltonian (H¯) at a particular level (say CCSD) and consider the EOM problem at any higher level (say EOM-CCSDT). In fact, there have been several studies of such nature.62 The opposite where the ground state is at a higher level of theory is also possible but may not have any practical value.

Beyond the approaches outlined above, higher ionization levels (DIP/TIP-EOM-CC) and electron-attached states (DEA/TEA-EOM-CC) can be treated analogously, and further modifications (e.g., the spin-flip or similarity transformed variants) of EOM-CC exist.58,63–67 As an aside, it should, furthermore, be noted that an excitation energy implementation of EOM-CC can be used to calculate ionized states by introducing an ultra-diffuse, continuum-like basis function.29,30 Importantly, this also means that there is a seamless transition from EE-EOM-CC to IP/EA-EOM-CC. This feature does not apply to most other quantum chemical methods.

All CC calculations were performed with the ACESII,68 ACESIII,69–71 and CFOUR72 program packages. Two benchmark datasets were compiled. First, experimental reference values for the valence IPs of 63 small organic molecules were taken from the work of Chong, Gritsenko, and Baerends, and references therein (referred to as the CGB-set in the following).73 Second, the HOMO IPs for 92 further organic molecules were taken from the NIST webbook (referred to as the NIST-set in the following). There is no overlap between the two sets.

For the CGB-set, geometries were optimized at the CCSD(T)/aug-cc-pVTZ level. IPs for the full set were calculated at the IP-EOM-CCSD level with cc-pVDZ, cc-pVTZ, and cc-pVQZ basis sets.74 Where computationally feasible, IP-EOM-CCSDT and IP-EOM-CCSDTQ calculations were also performed (using the cc-PVTZ basis set). This applies to 201 VIPs (from 41 molecules) for IP-EOM-CCSDT and 42 VIPs (14 molecules) for IP-EOM-CCSDTQ. For the NIST-set, geometry optimizations were carried out at the second-order many-body perturbation theory level [MBPT(2), MP2] with the def2-TZVP75 basis set. Here, we report 1067 VIPs at the IP-EOM-CCSD/cc-pVTZ level, bringing the total number of VIPs to 1468 covering 155 species.

We also include benchmarks for a small number of more approximate methods, tested on the CGB-set. Specifically, we consider the KS-eigenvalues of several DFT functionals, namely, LDA, B3LYP, CAM-B3LYP,76 CAM-QTP00,44 CAM-QTP01,45 CAM-QTP02,46 and QTP17,77 where the QTP functionals are specially designed to provide VIPs (based on Bartlett’s IP eigenvalue-theorem). Furthermore, we consider HF orbital eigenvalues (Koopmans’ theorem)78 and the extended Koopmans’ theorem (EKT)79–81 based on CCSD wavefunctions. Note that the CAM-QTP functionals are reparameterized versions of CAM-B3LYP76 and QTP17 is a re-parameterization of B3LYP, geared toward core ionizations and excitations, not valence IPs. For convenience, Table I summarizes the parameters for each functional. EKT and DFT calculation were conducted using ACESIII and NWChem,82 respectively.

TABLE I.

Definitions of the range-separated and global hybrid functionals used.

Range-separated:aαβμγ
CAM-B3LYP 0.19 0.46 0.33 0.81 
CAM-QTP00 0.54 0.37 0.29 0.80 
CAM-QTP01 0.23 0.77 0.31 0.80 
CAM-QTP02 0.28 0.72 0.34 1.00 
Global: α β  γ 
B3LYPb 0.20 0.72  0.81 
QTP17c 0.62 0.38  0.80 
Range-separated:aαβμγ
CAM-B3LYP 0.19 0.46 0.33 0.81 
CAM-QTP00 0.54 0.37 0.29 0.80 
CAM-QTP01 0.23 0.77 0.31 0.80 
CAM-QTP02 0.28 0.72 0.34 1.00 
Global: α β  γ 
B3LYPb 0.20 0.72  0.81 
QTP17c 0.62 0.38  0.80 
a

Exc=αExHF+1αβExB88,SR+βExHF,LR+γEcLYP+1γEcVWN5

b

Exc=αExHF+βExB88+(1β)ExLDA+γEcLYP+1γEcVWN5

c

Exc=αExHF+βExslater+γEcLYP+1γEcVWN5

In Sec. IV, we discuss the performance of each method in terms of mean signed and mean absolute deviations (MSDs/MADs). The full set of results is provided in the supplementary material, including the percentage of singles configurations that contribute to a given state in IP-EOM-CCSD. This is an important metric as IP-EOM-CCSD better describes states that are not primarily made up of singles configurations (see below).

The goal of this work is to provide a large set of IP-EOM-CC benchmark IPs. To this end, it is important to first establish how sensitive IP-EOM-CC results are relative to the neglect of higher excitation levels and basis-set incompleteness. We found that IP-EOM-CCSDT and IP-EOM-CCSDTQ results are in excellent agreement with each other (see Fig. 1, left). Specifically, the MAD between the two is 0.03 eV for a subset of 42 IPs from the CGB-set, indicating that effects beyond connected triples are essentially negligible in these calculations (see Tables S1 and S2 for statistics on all tested methods and basis sets).

FIG. 1.

Correlation plots between IPs computed with IP-EOM-CCSDT vs. IP-EOM-CCSDTQ (left) and IP-EOM-CCSD vs. IP-EOM-CCSDT (right).

FIG. 1.

Correlation plots between IPs computed with IP-EOM-CCSDT vs. IP-EOM-CCSDTQ (left) and IP-EOM-CCSD vs. IP-EOM-CCSDT (right).

Close modal

IP-EOM-CCSDT/cc-pVTZ calculations were performed for a total of 201 VIPs (41 molecules). For this subset, the much more efficient IP-EOM-CCSD method is also in good agreement with a MAD of 0.24 eV. If only the principal IPs (corresponding to the highest occupied molecular orbital, HOMO) are considered, the agreement is even better with a MAD of only 0.1 eV. Comparatively large deviations are only observed for the 4s state in HCl, the 1b1 state in ozone, and the 4b2 state in cis-C2H2F2 (i.e., the outliers in Fig. 1, right). This can be attributed to the low overall contribution of single excitations to these states (see below).

All tested IP-EOM-CC methods also show good agreement with the experimental reference values provided in the CGB paper. Specifically, IP-EOM-CCSDT and IP-EOM-CCSDTQ show very low MSDs (0.05 and 0.00 eV) and MADs of ∼0.2 eV. As discussed above, these results are converged with respect to the excitation level. This means that the remaining discrepancy between the theory and experiment is likely related to basis-set incompleteness, zero-point and non-adiabatic effects, as well as the uncertainty associated with the determination of vertical IPs from experiments. Some specific cases will be discussed further below.

With the same basis set (cc-pVTZ), the more approximate IP-EOM-CCSD method has a somewhat larger MSD of 0.19 eV and a MAD of 0.26 eV. As a matter of fact, the larger MSD indicates that IPs are systematically overestimated at this level of theory, a trend which is further aggravated when using the larger cc-pVQZ basis (MSD = 0.31 eV) (Fig. 2). While the higher-level methods display approximately the same performance for both lower-lying valence IPs and HOMO IPs, this is not the case at the CCSD level. In fact, IP-EOM-CCSD/cc-pVTZ is actually in better agreement with experiments for HOMO IPs than IP-EOM-CCSDTQ, with a MSD and MAD of 0.01 and 0.12 eV, respectively. This is due to fortuitous cancellation of errors with that particular basis set, as reflected in a somewhat worse performance of IP-EOM-CCSD at the cc-pVQZ level (MAD = 0.16 eV).

FIG. 2.

Correlation plots between experimental IPs and those computed with IP-EOM-CCSD using cc-pVTZ (black) and cc-pVQZ (red) basis sets. The outliers at 25 and 16.5 eV correspond to the 5a1 state in pyridine and the 1b1 state of ozone (see main text).

FIG. 2.

Correlation plots between experimental IPs and those computed with IP-EOM-CCSD using cc-pVTZ (black) and cc-pVQZ (red) basis sets. The outliers at 25 and 16.5 eV correspond to the 5a1 state in pyridine and the 1b1 state of ozone (see main text).

Close modal

In general, IP-EOM-CCSD performs better for HOMO IPs compared to lower-lying states, for all basis sets. This is because the HOMO ionizations display strong single excitation character. In Fig. 3, we plot the relative contribution of single excitations to all VIPs in the CGB-set (at the IP-EOM-CCSD/cc-pVQZ level). Out of 401 IPs, 31 consist of less than 85% single excitations and 9 of less than 60%. As shown in the figure, lower-lying states (higher VIPs) tend to display the lower single excitation character. Such ionizations are associated with a larger amount of orbital relaxation, which is not always accurately captured at the IP-EOM-CCSD level. In this respect, low-lying valence states are similar to core-ionized states, where the description of orbital relaxation is indeed the primary challenge for IP-EOM.83 In our experience, the threshold of 85% single excitation character applied above is a useful (conservative) rule of thumb for estimating the reliability of IP-EOM-CCSD in such cases.

FIG. 3.

Singles contribution for IPs at the IP-EOM-CCSD/cc-pVQZ level.

FIG. 3.

Singles contribution for IPs at the IP-EOM-CCSD/cc-pVQZ level.

Close modal

The correlation plot between IP-EOM-CCSD and experimental values for the NIST set is shown in Fig. 4. As these are exclusively HOMO VIPs, the predictions are generally of high quality, in line with the above discussion. Specifically, the MAD for all HOMO IPs is 0.15 eV and the MSD is 0.06 eV.

FIG. 4.

Correlation plots between experimental IPs and those computed at the IP-EOM-CCSD/cc-pVTZ level. Outliers at 9 eV correspond to the dimethylhydrazine isomers (see main text).

FIG. 4.

Correlation plots between experimental IPs and those computed at the IP-EOM-CCSD/cc-pVTZ level. Outliers at 9 eV correspond to the dimethylhydrazine isomers (see main text).

Close modal

Beyond this statistical analysis, it is worth discussing some cases in more detail. For our most accurate calculations (IP-EOM-CCSDTQ/cc-pVTZ), only four VIPs show deviations larger than 0.5 eV with respect to experimental data. These are (in descending order) the 1t2 state of CH4 (∆E = 0.76 eV), the 4σ state of HCl (∆E = 0.69 eV), the 4σ state of HCN (∆E = 0.50 eV), and the 1ε state of NH3 (∆E = 0.50 eV). In the case of CH4, the deviation is likely caused by the Jahn-Teller distortion of the ionized CH4 molecule, as reported by Porter et al.84 The Born-Oppenheimer approximation breaks down close to the triply degenerate tetrahedral ground state of CH4+, and the vertical approximation is thus not strictly valid. For the other cases, the cause of the discrepancies is not as clear, potentially pointing to inaccuracies in the experimental values. Indeed, highly accurate CC calculations have previously been used to clarify the assignment of experimental VIPs, e.g., for the proto-typical molecule, ethylene.53,85

There are also some interesting discrepancies for IP-EOM-CCSD. In particular, there is a large deviation for the 1b1 state of O3, with a predicted IP of 18.8 eV and a value of 19.99 eV reported in CGB. Counterintuitively, this deviation increases for IP-EOM-CCSDT, where the IP lies at 16.99 eV. Upon inspection of the original source cited in CGB,86 we conclude that the correct value should in fact be 16.50 eV. Another interesting case in the CGB-set is the low lying 5a1 state of pyridine. Here the predicted ionization energy is off by 6.8 eV at the IP-EOM-CCSD/cc-pVTZ level of theory. Surprisingly, this error is reduced to 1.9 eV when the larger cc-pVQZ basis set is used. Such a drastic basis-set effect is uncharacteristic for IPs, but it serves as a warning that the fortuitous cancellation of errors observed for the HOMOs with the TZ basis set does not necessarily translate to lower-lying states. In this case, the remaining error with the cc-pVQZ basis is still fairly large and can be attributed to the very low single excitation character of this state (47.6%). For ground-state energies, basis-set effects are routinely mitigated by extrapolation to the complete-basis-set (CBS) limit.87,88 This is less common for ionizations and excitation energies. In particular, we are not aware of a systematic study on CBS extrapolation for IP-EOM-CC although this would certainly be possible with the help of explicitly correlated reference calculations.89,90

In the NIST test set, the largest errors for IP-EOM-CCSD/cc-pVTZ are found for the HOMOs of dimethylhydrazine isomers (0.77 and 0.97 eV). Here, it is revealing to consider the experimental values for hydrazine and methylhydrazine. In both cases, different sources collected in the NIST webbook report vertical IPs that differ by 1-2 eV. Specifically, the hydrazine IP is reported to be 8.98, 10.68, and 9.90 eV, while for methylhydrazine 8.4 and 9.34 eV are given. This indicates that the experimental determination of vertical IPs for this class of compound is particularly difficult. Accordingly, the values for the dimethylhydrazines may also not be accurate to within more than the observed deviation with IP-EOM-CCSD.

A further interesting case is the splitting between the e2g(σ) and a1u(π) states (HOMO-1 and HOMO-2) in benzene. Here, GW methods predict a vanishingly small separation irrespective of the reference used. Meanwhile, the inclusion of (screened or bare) second-order exchange in GW leads to a very strong splitting of 1-2 eV.91 By contrast, the splitting is predicted to be ∼0.5 eV at the IP-EOM-CCSD/cc-pVTZ, in good agreement with experiments.

Having established the high accuracy and robustness of the IP-EOM-CC methods, we briefly consider the performance of alternative (more efficient) methods for the CGB-set. Specifically, we used the orbital eigenvalues of different DFT functionals and HF (Koopmans’ theorem), as well as the Extended Koopmans’ Theorem (EKT)92–94 method with CCSD densities (Table II).

TABLE II.

Mean signed deviation (MSD) and mean absolute deviation (MAD) of eigenvalues of LDA, B3LYP, CAM-B3LYP, CAM-QTP00, CAM-QTP01, CAM-QTP02, QTP17, HF, and EKT (CCSD) compared to experiments and IP-EOM-CCSD/cc-pVTZ basis sets. Results are reported in eV.

All valence satesHOMO
Vs. experimentNMSDMADNMSDMAD
LDA 401 −4.66 4.66 63 −4.33 4.33 
B3LYP 401 −3.26 3.26 63 −3.23 3.23 
CAM-B3LYP 401 −1.37 1.37 63 −1.50 1.50 
CAM-QTP00 401 0.76 0.81 63 0.16 0.32 
CAM-QTP01 401 0.13 0.40 63 −0.15 0.30 
CAM-QTP02 401 0.38 0.49 63 0.00 0.25 
QTP17 401 −0.18 0.49 63 −0.78 0.78 
HF 401 1.84 1.86 63 0.63 0.74 
EKT (CCSD) 401 0.97 0.97 63 0.34 0.35 
Vs. IPEOM-CCSD 
LDA 401 −4.84 4.84 63 −4.34 4.34 
B3LYP 401 −3.44 3.44 63 −3.23 3.23 
CAM-B3LYP 401 −1.55 1.55 63 −1.51 1.51 
CAM-QTP00 401 0.58 0.65 63 0.16 0.32 
CAM-QTP01 401 −0.05 0.31 63 −0.16 0.29 
CAM-QTP02 401 0.20 0.34 63 0.00 0.24 
QTP17 401 −0.36 0.47 63 −0.79 0.79 
HF 401 1.66 1.69 63 0.63 0.76 
EKT (CCSD) 401 0.79 0.80 63 0.34 0.34 
All valence satesHOMO
Vs. experimentNMSDMADNMSDMAD
LDA 401 −4.66 4.66 63 −4.33 4.33 
B3LYP 401 −3.26 3.26 63 −3.23 3.23 
CAM-B3LYP 401 −1.37 1.37 63 −1.50 1.50 
CAM-QTP00 401 0.76 0.81 63 0.16 0.32 
CAM-QTP01 401 0.13 0.40 63 −0.15 0.30 
CAM-QTP02 401 0.38 0.49 63 0.00 0.25 
QTP17 401 −0.18 0.49 63 −0.78 0.78 
HF 401 1.84 1.86 63 0.63 0.74 
EKT (CCSD) 401 0.97 0.97 63 0.34 0.35 
Vs. IPEOM-CCSD 
LDA 401 −4.84 4.84 63 −4.34 4.34 
B3LYP 401 −3.44 3.44 63 −3.23 3.23 
CAM-B3LYP 401 −1.55 1.55 63 −1.51 1.51 
CAM-QTP00 401 0.58 0.65 63 0.16 0.32 
CAM-QTP01 401 −0.05 0.31 63 −0.16 0.29 
CAM-QTP02 401 0.20 0.34 63 0.00 0.24 
QTP17 401 −0.36 0.47 63 −0.79 0.79 
HF 401 1.66 1.69 63 0.63 0.76 
EKT (CCSD) 401 0.79 0.80 63 0.34 0.34 

HF eigenvalues tend to systematically overestimate VIPs as the orbitals are optimized for the neutral state (i.e., orbital relaxation in the ionized state is neglected). While there is some cancellation of errors due to the neglect of electron correlation, this only works well for HOMO IPs. Meanwhile, the eigenvalues of conventional KS-DFT methods systematically underestimate VIPs. By contrast, the range-separated and global hybrid functionals of the QTP family actually outperform the more elaborate EKT method. Out of the tested density functionals, CAM-QTP01 performed best, followed by CAM-QTP00. In contrast to most other methods, the QTP17 functional shows larger deviations for HOMOs than for lower-lying states. This is a consequence of the QTP17 parameterization procedure, where core states were more strongly weighted. Nevertheless, it performs very well for a global hybrid functional (compared with B3LYP). Considering the computational cost and accuracy, the QTP functionals are thus an excellent choice for screening in cases where IP-EOM-CC is prohibitively expensive to apply or might be overkill.

We have presented a benchmark database of 1468 vertical IPs covering 155 species at the IP-EOM-CCSD level of theory. We expect that this database will be a valuable resource for developers of electronic structure methods in Green’s function, GW, density functional, and semiempirical communities. The accuracy of IP-EOM-CCSD was established by comparing with experiment and higher-level CC (IP-EOM-CCSDT and IP-EOM-CCSDTQ) reference values.

Our workhorse method is found to have a MAD of ∼0.1 eV for HOMOs and ∼0.25 eV for all valence IPs. Particularly, the combination IP-EOM-CCSD/cc-pVTZ displays favorable cancellation of errors for HOMO states, where it actually outperforms IPEOM-CCSDTQ. In a few cases, uncharacteristically large deviations (0.5-1 eV) from experiments were found. This is the case for certain lower-lying states, where orbital relaxation is large. Additionally, the vertical approximation can break down in certain cases, such as the CH4 HOMO. For two states (HCl 4σ and O3 1b1), even IP-EOM-CCSDTQ shows large errors, potentially indicating that the experimental references are unreliable. For high-throughput or large-scale applications, the density functionals of the QTP family are found to be an efficient and reasonably accurate alternative to IPEOM-CCSD.

This paper does not directly assess core ionization, but prior work has shown that IP-EOM-CCSD in a cc-pCVTZ basis provides an accuracy of ∼2 eV (out of hundreds of eV) for core ionizations of light atoms in molecules95 (see also Table S6). Adding further orbital relaxation to the calculations reduces the error to ∼1 eV. QTP17 and CAM-QTP00 have been shown to be exceptionally accurate compared to other KS-DFT approximations for core ionizations with results ∼1 eV.46 These results account for the massive contribution of orbital relaxation to core-hole ionization. It is certainly the simplest way to obtain accurate core ionization for ESCA or XPS from a single calculation on the neutral. Among other applications, this should mean that relaxation errors that occur in some of the lower-lying VIPs might actually be improved compared to the CC results for outliers. For heavier, second-row atoms, core excitation values are discussed in Ref. 42.

See supplementary material for geometries, vertical ionization computed using IP-EOM-CCSD, IP-EOM-CCSDT, and IP-EOM-CCSDTQ for species in the CGB test set, vertical ionization computed using IPEOM-CCSD for species in the NIST test set, vertical ionization computed using LDA, B3LYP, CAM-B3LYP, CAM-QTP00, CAM-QTP01, CAM-QTP02, QTP17, HF, and EKT (CCSD) for species in the CGB test set, core ionization calculation computed using QTP17, CAM-QTP00, and IP-EOM CC methods, and mean signed deviation (MSD) and mean absolute deviation (MAD) of eigenvalues of density functionals, HF, and EKT (CCSD) compared to IP-EOM-CCSDT and IP-EOM-CCSDTQ methods.

This work was sponsored by the U.S. Air Force Office of Scientific Research, Under Grant No. FA9550-11-1-0065. The authors thank Professor John Stanton for providing access to CFOUR program package and the University of Florida Research Computing for providing computational resources and support that have contributed to the research results reported in this publication. URL: http://researchcomputing.ufl.edu.

1.
L.
Pauling
,
J. Am. Chem. Soc.
54
,
3570
(
1932
).
2.
R. S.
Mulliken
,
J. Chem. Phys.
3
,
573
(
1935
).
3.
R. G.
Parr
and
R. G.
Pearson
,
J. Am. Chem. Soc.
105
,
7512
(
1983
).
4.
R. G.
Parr
,
L. V.
Szentpály
, and
S.
Liu
,
J. Am. Chem. Soc.
121
,
1922
(
1999
).
5.
R. G.
Parr
and
W.
Yang
,
J. Am. Chem. Soc.
106
,
4049
(
1984
).
6.
E.
Hückel
,
Z. Phys.
70
,
204
(
1931
).
7.
R.
Hoffmann
,
J. Chem. Phys.
39
,
1397
(
1963
).
8.
R.
Pariser
and
R. G.
Parr
,
J. Chem. Phys.
21
,
466
(
1953
).
9.
J. A.
Pople
,
Trans. Faraday Soc.
49
,
1375
(
1953
).
10.
R.
Pariser
,
J. Chem. Phys.
21
,
568
(
1953
).
11.
J. L.
Bredas
,
Mater. Horiz.
1
,
17
(
2014
).
12.
J. E.
Anthony
,
Chem. Rev.
106
,
5028
(
2006
).
13.
A.
Nitzan
and
M. A.
Ratner
,
Science
300
,
1384
(
2003
).
14.
A.
Aviram
and
M. A.
Ratner
,
Chem. Phys. Lett.
29
,
277
(
1974
).
15.
C. J.
Brabec
,
N. S.
Sariciftci
, and
J. C.
Hummelen
,
Adv. Funct. Mater.
11
,
15
(
2001
).
16.
W.
Shockley
and
H. J.
Queisser
,
J. Appl. Phys.
32
,
510
(
1961
).
17.
M. C.
Scharber
,
D.
Mühlbacher
,
M.
Koppe
,
P.
Denk
,
C.
Waldauf
,
A. J.
Heeger
, and
C. J.
Brabec
,
Adv. Mater.
18
,
789
(
2006
).
18.
J.
Hachmann
,
R.
Olivares-Amaya
,
A.
Jinich
,
A. L.
Appleton
,
M. A.
Blood-Forsythe
,
L. R.
Seress
,
C.
Román-Salgado
,
K.
Trepte
,
S.
Atahan-Evrenk
,
S.
Er
,
S.
Shrestha
,
R.
Mondal
,
A.
Sokolov
,
Z.
Bao
, and
A.
Aspuru-Guzik
,
Energy Environ. Sci.
7
,
698
(
2014
).
19.
C.
Schober
,
K.
Reuter
, and
H.
Oberhofer
,
J. Phys. Chem. Lett.
7
,
3973
(
2016
).
20.
J. F.
Stanton
and
R. J.
Bartlett
,
J. Chem. Phys.
98
,
7029
(
1993
).
21.
D. C.
Comeau
and
R. J.
Bartlett
,
Chem. Phys. Lett.
207
,
414
(
1993
).
22.
J.
Geertsen
,
M.
Rittby
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
164
,
57
(
1989
).
23.
J.
Geertsen
and
J.
Oddershede
,
J. Chem. Phys.
85
,
2112
(
1986
).
24.
H.
Sekino
and
R. J.
Bartlett
,
Int. J. Quantum Chem.
26
,
255
(
1984
).
25.
K.
Emrich
,
Nucl. Instrum. Methods Phys. Res., Sect. A
351
,
379
(
1981
).
26.
R. J.
Bartlett
,
Wiley Interdiscip. Rev. Comput. Mol. Sci.
2
,
126
(
2012
).
27.
H. J.
Monkhorst
,
Int. J. Quantum Chem.
12
,
421
(
1977
).
28.
J. F.
Stanton
and
J.
Gauss
,
J. Chem. Phys.
101
,
8938
(
1994
).
29.
M.
Nooijen
and
R. J.
Bartlett
,
J. Chem. Phys.
102
,
3629
(
1995
).
30.
J. F.
Stanton
and
J.
Gauss
,
J. Chem. Phys.
111
,
8785
(
1999
).
31.
R. J.
Bartlett
,
V. F.
Lotrich
, and
I. V.
Schweigert
,
J. Chem. Phys.
123
,
062205
(
2005
).
32.
M. J.
Van Setten
,
F.
Caruso
,
S.
Sharifzadeh
,
X.
Ren
,
M.
Scheffler
,
F.
Liu
,
J.
Lischner
,
L.
Lin
,
J. R.
Deslippe
,
S. G.
Louie
,
C.
Yang
,
F.
Weigend
,
J. B.
Neaton
,
F.
Evers
, and
P.
Rinke
,
J. Chem. Theory Comput.
11
,
5665
(
2015
).
33.
X.
Ren
,
P.
Rinke
,
V.
Blum
,
J.
Wieferink
,
A.
Tkatchenko
,
A.
Sanfilippo
,
K.
Reuter
, and
M.
Scheffler
,
New J. Phys.
14
,
053020
(
2012
).
34.
T.
Rangel
,
S. M.
Hamed
,
F.
Bruneval
, and
J. B.
Neaton
,
J. Chem. Theory Comput.
12
,
2834
(
2016
).
35.
M. F.
Lange
and
T. C.
Berkelbach
,
J. Chem. Theory Comput.
14
,
4224
(
2018
).
36.
T.
Körzdörfer
,
R. M.
Parrish
,
N.
Marom
,
J. S.
Sears
,
C. D.
Sherrill
, and
J.-L.
Brédas
,
Phys. Rev. B
86
,
205110
(
2012
).
37.
R.
Baer
,
E.
Livshits
, and
U.
Salzner
,
Annu. Rev. Phys. Chem.
61
,
85
(
2010
).
38.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
N.
Govind
,
J.
Autschbach
,
J. B.
Neaton
,
R.
Baer
, and
L.
Kronik
,
Phys. Rev. Lett.
109
,
226405
(
2012
).
39.
M.
Dauth
,
F.
Caruso
,
S.
Kümmel
, and
P.
Rinke
,
Phys. Rev. B
93
,
121115
(
2016
).
40.
P.
Verma
and
R. J.
Bartlett
,
J. Chem. Phys.
137
,
134102
(
2012
).
41.
R. J.
Bartlett
,
Chem. Phys. Lett.
484
,
1
(
2009
).
42.
U.
Salzner
and
R.
Baer
,
J. Chem. Phys.
131
,
231101
(
2009
).
43.
R. J.
Bartlett
and
D. S.
Ranasinghe
,
Chem. Phys. Lett.
669
,
54
(
2017
).
44.
P.
Verma
and
R. J.
Bartlett
,
J. Chem. Phys.
140
,
18A534
(
2014
).
45.
Y.
Jin
and
R. J.
Bartlett
,
J. Chem. Phys.
145
,
034107
(
2016
).
46.
Y.
Jin
and
R. J.
Bartlett
,
J. Chem. Phys.
149
,
064111
(
2018
).
47.
J. T.
Margraf
,
P.
Verma
, and
R. J.
Bartlett
,
J. Chem. Phys.
145
,
104106
(
2016
).
48.
J. T.
Margraf
,
D.
Claudino
, and
R. J.
Bartlett
,
Mol. Phys.
115
,
538
(
2017
).
49.
C.
Li
,
X.
Zheng
,
N. Q.
Su
, and
W.
Yang
,
Natl. Sci. Rev.
5
,
203
(
2018
).
50.
E. R.
Davidson
and
A. A.
Jarzȩcki
,
Chem. Phys. Lett.
285
,
155
(
1998
).
51.
L. S.
Cederbaum
and
W.
Domcke
,
J. Chem. Phys.
64
,
603
(
1976
).
52.
L. S.
Cederbaum
and
W.
Domcke
,
Adv. Chem. Phys.
36
,
205
(
1977
).
53.
A. D.
Yau
,
S. A.
Perera
, and
R. J.
Bartlett
,
Mol. Phys.
100
,
835
(
2002
).
54.
K.
Krause
,
M. E.
Harding
, and
W.
Klopper
,
Mol. Phys.
113
,
1952
(
2015
).
55.
J. W.
Knight
,
X.
Wang
,
L.
Gallandi
,
O.
Dolgounitcheva
,
X.
Ren
,
J. V.
Ortiz
,
P.
Rinke
,
T.
Körzdörfer
, and
N.
Marom
,
J. Chem. Theory Comput.
12
,
615
(
2016
).
56.
R. M.
Richard
,
M. S.
Marshall
,
O.
Dolgounitcheva
,
J. V.
Ortiz
,
J. L.
Brédas
,
N.
Marom
, and
C. D.
Sherrill
,
J. Chem. Theory Comput.
12
,
595
(
2016
).
57.
M.
Nooijen
and
J. G.
Snijders
,
Int. J. Quantum Chem.
44
,
55
(
1992
).
58.
M.
Nooijen
and
V.
Lotrich
,
J. Chem. Phys.
113
,
494
(
2000
).
59.
S.
Hirata
,
M.
Nooijen
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
328
,
459
(
2000
).
60.
M.
Musiał
,
A.
Perera
, and
R. J.
Bartlett
,
J. Chem. Phys.
134
,
114108
(
2011
).
61.
M.
Musial
and
R. J.
Bartlett
,
J. Chem. Phys.
119
,
1901
(
2003
).
62.
J. N.
Byrd
,
V.
Rishi
,
A.
Perera
, and
R. J.
Bartlett
,
J. Chem. Phys.
143
,
164103
(
2015
).
63.
J. D.
Watts
and
R. J.
Bartlett
,
Chem. Phys. Lett.
233
,
81
(
1995
).
64.
M.
Musiał
,
M.
Olszówka
,
D. I.
Lyakh
, and
R. J.
Bartlett
,
J. Chem. Phys.
137
,
174102
(
2012
).
65.
D.
Casanova
,
L. V.
Slipchenko
,
A. I.
Krylov
, and
M.
Head-Gordon
,
J. Chem. Phys.
130
,
044103
(
2009
).
66.
A.
Perera
,
R. W.
Molt
,
V. F.
Lotrich
, and
R. J.
Bartlett
,
Theor. Chem. Acc.
133
,
1514
(
2014
).
67.
M.
Nooijen
and
R. J.
Bartlett
,
J. Chem. Phys.
106
,
6441
(
1997
).
68.
J. F.
Stanton
,
J.
Gauss
,
A.
Perera
,
A.
Yau
,
J. D.
Watts
,
M.
Nooijen
,
N.
Oliphant
,
P. G.
Szalay
,
W. J.
Lauderdale
,
S. R.
Gwaltney
,
S.
Beck
,
A.
Balkov
,
D. E.
Bernholdt
,
K.-K.
Baeck
,
P.
Rozyczko
,
C.
Sekino
,
H.
Huber
,
J.
Pittner
,
R. J.
Bartlett
, and
P. R.
Taylor
, ACES II is a product of the Quantum Theory Project, University of Florida. Integral packages included are VMOL (J. Almlöf and P. R. Taylor), VPROPS (P. R. Taylor), and ABACUS (T. Helgaker, H. J. Aa. Jensen, P. Jørgensen, J. Olsen, and P. R. Taylor),
2009
.
69.
V. F.
Lotrich
,
J. M.
Ponton
,
A. S.
Perera
,
E.
Deumens
,
R. J.
Bartlett
, and
B. A.
Sanders
,
Mol. Phys.
108
,
3323
(
2010
).
70.
V.
Lotrich
,
N.
Flocke
,
M.
Ponton
,
A. D.
Yau
,
A.
Perera
,
E.
Deumens
, and
R. J.
Bartlett
,
J. Chem. Phys.
128
,
194104
(
2008
).
71.
E.
Deumens
,
V. F.
Lotrich
,
A.
Perera
,
M. J.
Ponton
,
B. A.
Sanders
, and
R. J.
Bartlett
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
895
(
2011
).
72.
J. F.
Stanton
,
J.
Gauss
,
M. E.
Harding
,
P. G.
Szalay
,
A. A.
Auer
,
R. J.
Bartlett
,
U.
Benedikt
,
C.
Berger
,
D. E.
Bernholdt
,
Y. J.
Bomble
,
L.
Cheng
,
O.
Christiansen
,
M.
Heckert
,
O.
Heun
,
C.
Huber
,
T.-C.
Jagau
,
D.
Jonsson
,
J.
Juselius
,
K.
Klein
,
W. J.
Lauderdale
,
T. M. D. A.
Matthews
,
L. A.
Muck
,
D. P.
ONeill
,
D. R.
Price
,
E.
Prochnow
,
C.
Puzzarini
,
W. S. K.
Ruud
,
F.
Schiffmann
,
C.
Simmons
,
S.
Stopkowicz
,
A.
Tajti
,
J.
Vazquez
, and
F.
Wang
, CFOUR, Coupled-Cluster Techniques for Computational Chemistry, a Quantum-Chemical Program Package,
2010
, http://www.cfour.de.
73.
D. P.
Chong
,
O. V.
Gritsenko
, and
E. J.
Baerends
,
J. Chem. Phys.
116
,
1760
(
2002
).
74.
T. H.
Dunning
, Jr
,
J. Chem. Phys.
90
,
1007
(
1989
).
75.
F.
Weigend
and
R.
Ahlrichs
,
Phys. Chem. Chem. Phys.
7
,
3297
(
2005
).
76.
T.
Yanai
,
D. P.
Tew
, and
N. C.
Handy
,
Chem. Phys. Lett.
393
,
51
(
2004
).
77.
R. L. A.
Haiduke
and
R. J.
Bartlett
,
J. Chem. Phys.
149
,
131101
(
2018
).
78.
T.
Koopmans
,
Physica
1
,
104
(
1934
).
79.
M. D.
Newton
,
J. Chem. Phys.
48
,
2825
(
1968
).
80.
D. W.
Smith
,
J. Chem. Phys.
62
,
113
(
1975
).
81.
O. W.
Day
,
D. W.
Smith
, and
R. C.
Morrison
,
J. Chem. Phys.
62
,
115
(
1975
).
82.
M.
Valiev
,
E. J.
Bylaska
,
N.
Govind
,
K.
Kowalski
,
T. P.
Straatsma
,
H. J. J.
Van Dam
,
D.
Wang
,
J.
Nieplocha
,
E.
Apra
,
T. L.
Windus
, and
W. A.
De Jong
,
Comput. Phys. Commun.
181
,
1477
(
2010
).
83.
T. J.
Watson
and
R. J.
Bartlett
,
Chem. Phys. Lett.
555
,
235
(
2013
).
84.
A. R.
Porter
,
O. K.
Al-Mushadani
,
M. D.
Towler
, and
R. J.
Needs
,
J. Chem. Phys.
114
,
7795
(
2001
).
85.
M.
Musia
and
R. J.
Bartlett
,
Chem. Phys. Lett.
384
,
210
(
2004
).
86.
S.
Katsumata
,
H.
Shiromaru
, and
T.
Kimura
,
Bull. Chem. Soc. Jpn.
57
,
1784
(
1984
).
87.
A.
Halkier
,
T.
Helgaker
,
P.
Jørgensen
,
W.
Klopper
,
H.
Koch
,
J.
Olsen
, and
A. K.
Wilson
,
Chem. Phys. Lett.
286
,
243
(
1998
).
88.
D. W.
Schwenke
,
J. Chem. Phys.
122
,
014107
(
2005
).
89.
D.
Bokhan
,
D. N.
Trubnikov
,
M.
Musiał
, and
R. J.
Bartlett
,
Chem. Phys. Lett.
610-611
,
173
(
2014
).
90.
D.
Bokhan
and
S.
Ten-no
,
J. Chem. Phys.
132
,
021101
(
2010
).
91.
X.
Ren
,
N.
Marom
,
F.
Caruso
,
M.
Scheffler
, and
P.
Rinke
,
Phys. Rev. B
92
,
81104
(
2015
).
92.
J.
Cioslowski
,
P.
Piskorz
, and
G.
Liu
,
J. Chem. Phys.
107
,
6804
(
1997
).
93.
M.
Ernzerhof
,
J. Chem. Theory Comput.
5
,
793
(
2009
).
94.
D.
Vanfleteren
,
D.
Van Neck
,
P. W.
Ayers
,
R. C.
Morrison
, and
P.
Bultinck
,
J. Chem. Phys.
130
,
194104
(
2009
).
95.
T. J.
Watson
, Jr.
,
Extending the Capabilities of Accurate Ab Initio Methods: Novel Algorithms and Massively Parallelizable Implementations for Feasible Calculations
(
University of Florida
,
2012
).
96.

Please note that the preferred IUPAC nomenclature is ionization energy (IE) instead of ionization potential. As the name of the IP-EOM-CC method is firmly established, however, we will use IP throughout this manuscript, to avoid confusion.

Supplementary Material