Towards Chemical Accuracy in Atomic Ionization Energies: The Case of H, C, N, O, F, P, and S Atoms

1. Introduction

The calculation of ionization energies of atoms may be useful for understanding their oxidation mechanism, for benchmarking quantum-chemical methods, and for modeling ionization processes in chemical and biological systems.

Two types of ionization energy (i.e., vertical and adiabatic) can be defined for atomic or molecular systems. The vertical ionization energy assumes that, during the removal of an electron, the orbitals of the system remain unchanged. By contrast, in the case of the adiabatic ionization energy, the orbitals relax during the ionization process as a result of the redistribution of the electronic density. Experimentally, ionization energies can be determined, for example, by photoelectron spectroscopy [1], from the difference between the energy of the incident photons and the kinetic energy of the emitted photoelectrons.

We computed the ionization energies of several atoms commonly found in proteins, namely H, C, N, O, and S, as well as two additional atoms, F and P. The atoms H, C, N, O, and S are present in essential amino acids, and the calculation of their ionization potentials may contribute to a better understanding of ionization processes in proteins. The atoms F and P may occur in synthetic amino acids and exhibit properties that can be exploited in the pharmaceutical industry [2,3]. Because fluorine is the most electronegative element, it forms strong bonds with carbon and can therefore be used in protein engineering for the synthesis of highly stable proteins [4]. Phosphorus-containing amino acids also have a variety of applications. For example, they may act as protecting groups, serve as antiviral agents incorporated into prodrugs, or be used for ¹⁸F labeling [5,6].

The methods employed in this work to compute ionization energies include the energy-difference approach and the electron propagator theory (EPT) methods, namely the Outer Valence Green’s Function (OVGF) and the renormalized partial third-order quasiparticle (P3+) methods. The energy-difference approach was combined with the B2PLYP DFT method and various post-Hartree–Fock methods. In addition, the coupled-cluster methods CCSD(T), CCSDT, and IP-EOM-CCSD were employed, in conjunction with aug-cc-pVXZ-DK, aug-cc-pVXZ, and ANO-RCC basis sets families, with all-electron correlation treatment, the second-order Douglas-Kroll-Hess (DKH2) scalar-relativistic Hamiltonian, and extrapolation to the complete basis set (CBS) limit using the two-point Helgaker formula. The aim of these calculations was to achieve chemical accuracy, defined here as a maximum deviation of 1 kcal/mol between the experimental and calculated ionization energies (1 kcal/mol = 4.184 kJ/mol = 0.0434 eV) [7].

2. Computational Methods

In the energy-difference approach, the ionization potential is obtained as the energy difference between the cation, generated by removing one electron, and the corresponding neutral system [8]:

$$I=E_{cation}-E_{neutral}$$

(1)

To obtain the adiabatic ionization energy, both the neutral system and the cation are optimized, and the difference between their total energies is then evaluated. By contrast, the vertical ionization energy is obtained by calculating the energy of the cation at the geometry of the neutral system. The approach was previously used by Jursic [9] to compute the ionization potentials of several second-row elements (C, N, O, F) employing DFT and ab initio methods.

The Electron Propagator Theory (EPT) [10] is a useful framework for calculation of ionization energies, electron affinities, one-electron properties, or spectral intensities. EPT is based on the concept of propagator. For an electron in a quantum system, the propagator is given by the expectation value of the time-evolution operator between the initial and final states,

$$G\left({\overrightarrow{r}}_i,{\overrightarrow{r}}_f,t,t^{\text{'}}\right)=⟨{\overrightarrow{r}}_f\left|\left.\widehat{U}\left(t,t^{\text{'}}\right)\right|{\overrightarrow{r}}_i\right.⟩$$

(2)

where $\widehat{U}\left(t,t^{\text{'}}\right)$ is the time evolution operator. The one-particle Green’s function [11] can be used to calculate electron propagators. EPT is based on the solution of the Dyson equation [12,13,14],

$$G=G^0+G^0\mathsf{\Sigma }G$$

(3)

where $G$ is the one-particle Green’s function, $G^0$ is the free Green’s function, and $\mathsf{\Sigma }$ is the self-energy which contains the correlation and final-state relaxation effects [15,16]. The Dyson equation may be transformed in the following one-electron equation [16,17,18]:

$$\left[F+\mathsf{\Sigma }\left(E\right)\right]{\varphi }^{Dyson}=E{\varphi }^{Dyson}$$

(4)

where $F$ is the one-electron Fock operator, $E$ is the energy, and ${\varphi }^{Dyson}$ is the Dyson orbital [19]. The Dyson orbital corresponding to the ionization from the initial state ${\mathsf{\Psi }}_N(x_1,\dots ,x_N)$ to the final state ${\mathsf{\Psi }}_{s,N-1}(x_2,\dots ,x_N)$ is defined as [17]:

$${\varphi }_s^{Dyson}\left(x_1\right)=\sqrt{N}\int {\mathsf{\Psi }}_N\left(x_1,x_2,\dots ,x_N\right){\mathsf{\Psi }}_{s,N-1}^{*}\left(x_2,x_3,\dots ,x_N\right)d{x}_{2}d{x}_3\dots d{x}_N$$

(5)

where $x_i$ stands for the space-spin coordinate of electron i.

Based on the Dyson orbital, the pole strength, which takes values in the [0,1] interval, can be computed as follows [18,19]:

$$p_s=\int {\left|{\varphi }_s^{Dyson}\left(x\right)\right|}^{2}dx$$

(6)

For strong correlation the pole strength is close to zero, whereas for weak correlation it approaches unity [16]. When the Koopmans picture provides a qualitatively valid description of the ${\mathsf{\Psi }}_N$ and ${\mathsf{\Psi }}_{s,N-1}$ states, the pole strength is greater than 0.85 [16].

In this paper we employ two diagonal approximations of the EPT, namely OVGF and P3+, in which the off-diagonal elements of the self-energy matrix are neglected [20]. The OVGF method is designed for outer valence electrons and is based on a perturbative expansion of the self-energy that is exact through third order, while higher order contributions are treated approximately [14]. There are three variants of the OVGF method, denoted A, B, and C, which differ in the form of the self-energy expression. The corresponding expressions for the self-energy are given in ref. [10]. Method A is applied when the ionization energy is greater than 15 eV, whereas method B is used when the ionization energy is smaller than 15 eV [20]. Approximations A and B correspond to the cases in which the second-order terms in the self-energy exceed the third-order terms, whereas approximation C is used when the second order terms are small [20]. An algorithm for selecting among methods A, B, and C is listed in ref. [20]. The computational costs of the OVGF method scales as OV⁴, where O is the number of occupied orbitals and V is the number of virtual orbitals [15,16,21]. For closed-shell systems, the absolute error in computed vertical ionization potential below 20 eV is about 0.25 eV [20,22].

The P3 method [23] neglects some terms from the self-energy expression used in OVGF. A renormalized version of P3, denoted P3+, has proven to be very effective for computing the ionization energies of challenging anions [17]. The computational cost of both P3 and P3+ scales as O³V², which represents an improvement over OVGF method because the number of occupied orbitals is usually much smaller than the number of virtual orbitals [16]. For P3+, mean absolute deviations of less than 0.2 eV between calculated and experimental ionization energies have been reported for anions with ionization potentials below 20 eV [17].

For the seven atoms considered in this work, the energy-difference approach was used in conjunction with the DFT functional B2PLYP [24], as well as the post-Hartree–Fock methods MP2 [25], CCSD(T) [26,27], QCISD(T) [27], G3 [28], and CBS-QB3 [29]. These methods were combined with the aug-cc-pVQZ basis set. The OVGF and P3+ calculations were performed with the 6-311+G(2df,p), cc-pVQZ, and aug-cc-pVQZ basis sets. All these calculations were carried out with Gaussian 16 Revision C.01 [30] on the High-Performance Computing Center of Babeș-Bolyai University [31].

In addition to the calculations described above, the coupled-cluster methods CCSD(T) [26,27], CCSDT [32], and IP-EOM-CCSD [33,34] were used for the calculation of the ionization energies of the C, N, O, F, P, and S atoms, as implemented in ORCA 6.1.0 [35]. The CCSD(T) method combines the iterative coupled-cluster treatment of single and double excitations with a perturbative correction of the triple excitations. The CCSDT method provides an iterative treatment of the triple excitations and was employed in order to evaluate the quality of the perturbative correction in the CCSD(T). The CCSDT calculations were performed using the AUTOCI module implemented in ORCA. The IP-EOM-CCSD method is an equation-of-motion (EOM) coupled-cluster approach for ionization energies. It is based on a neutral CCSD reference state, while the ionized states are described through electron removal operators in the EOM space. Therefore, unlike the energy-difference approach, IP-EOM-CCSD gives the ionization energies directly. For the hydrogen atom, which is a one-electron system, only Hartree-Fock (HF) calculations were performed. Very tight criterion was used for the convergence of SCF calculations.

Three families of basis sets were used: the aug-cc-pVXZ basis sets (X = Q, 5) [36,37], the corresponding relativistically recontracted aug-cc-pVXZ-DK basis sets (X = Q, 5) [38], and the relativistically contracted atomic natural orbitals basis sets (ANO-RCC) at the TZP, QZP, and Full contraction levels [39]. The aug-cc-pVXZ-DK and ANO-RCC calculations were performed in conjunction with the second order Douglas-Kroll-Hess scalar relativistic Hamiltonian (DKH2) [40,41] and the finite nucleus model, while the aug-cc-pVXZ calculations were done without relativistic corrections. All electrons were included in the correlation treatment through the NoFrozenCore option. For the CCSD(T) and CCSDT calculations, the ionization energies were obtained as the energy differences between the neutral and cationic atoms. For the IP-EOM-CCSD calculations, the ionization energies were obtained directly from the ionized roots of the neutral reference state. For the six atoms with more than one electron, the α and β electron removal roots were treated separately, and the first ionization energy was considered to be the lower value.

In order to reduce the basis set error, the ionization energies obtained with CCSD(T), CCSDT, and IP-EOM-CCSD were extrapolated to the complete basis set (CBS) limit using the two-point formula of Helgaker et al. [42,43]:

$$E_X=E_{CBS}+{\displaystyle \frac{A}{X^3}}$$

(7)

which results in

$$E_{CBS}={\displaystyle \frac{X^3{E}_X-Y^3{E}_Y}{X^3-Y^3}}$$

(8)

where X and Y are the cardinal numbers of the basis sets used in the extrapolation. For the aug-cc-pVXZ-DK and aug-cc-pVXZ basis sets, the (Q,5) pair was used, while for the ANO-RCC basis sets the (T,Q) pair was employed based on the TZP and QZP contractions. For the CCSD(T) and CCSDT methods, a dual-level CBS extrapolation scheme was used, in which the total energy at each basis set level was separated into a self-consistent field (SCF) component and a correlation component:

$$E_{tot}\left(X\right)=E_{SCF}\left(X\right)+E_{corr}\left(X\right)$$

(9)

This separation was justified by the fact that the two contributions exhibit different convergence patterns with increasing the basis set size [43,44]. The SCF energy converges exponentially with the cardinal number and is already close to the CBS limit at the 5Z level, while the correlation energy converges asymptotically with X^-3. Therefore, the ionization energy at the CBS limit was obtained as:

$${IE}_{SCF}\left(Y\right)=E_{SCF}^{cat}\left(Y\right)-E_{SCF}^{neu}\left(Y\right)$$

(10)

$${IE}_{corr}\left(X\right)=E_{corr}^{cat}\left(X\right)-E_{corr}^{neu}\left(X\right)$$

(11)

$${IE}_{corr,CBS}={\displaystyle \frac{Y^3{IE}_{corr}\left(Y\right)-X^3{IE}_{corr}\left(X\right)}{Y^3-X^3}}$$

(12)

$${IE}_{CBS}={IE}_{SCF}\left(Y\right)+{IE}_{corr,CBS}$$

(13)

where Y = X+1 and the SCF contribution is taken from the larger basis set, namely aug-cc-pV5Z, aug-cc-pV5Z-DK, and ANO-RCC-QZP. This scheme was used by Richard et al. [44] for the calculation of ionization potentials and electron affinities at the CCSD(T)/CBS limit. For the IP-EOM-CCSD method, which computes the ionization energy directly, the extrapolation was applied directly to the α and β electron removal ionization energies [45]:

$${IP}_{s,CBS}={\displaystyle \frac{Y^3{IE}_s\left(Y\right)-X^3{IE}_s\left(X\right)}{Y^3-X^3}},s=\alpha ,\beta$$

(14)

The first ionization energy at the CBS limit was then obtained as the lower of the two extrapolated values.

3. Results and Discussion

3.1. Ionization Energies for the C, N, O, F, P, and S Atoms Obtained with the OVGF and P3+ Methods

Table 1 summarizes the ionization energies of the six atoms C, N, O, F, P, and S calculated using the OVGF and P3+ methods. The hydrogen atom was not considered within these approaches because, having only one electron, electron correlation is absent. Three basis sets were used: cc-pVQZ, aug-cc-pVQZ, and 6-311+G(2df,p). The table reports both the experimental ionization energies and the values computed with the OVGF and P3+ methods. For OVGF, the A, B, and C variants were evaluated according to the algorithm implemented in Gaussian package, and the recommended values were used in the discussion. Beside these quantities, the table includes the spin multiplicities of the neutral atoms, the indexes of the HOMO and LUMO orbitals, and the orbital window employed in the calculations.

Figure 1 sketch the atomic orbital levels of the C, N, O, F, P, and S atoms used in the OVGF calculations, while Figure 2 exemplifies the orbital indexing for the OVGF method in case of the sulfur atom.

Table 2 contains the MAE and MSE values averaged over the six atoms for the three basis sets with the OVGF and P3+ methods. Chemical accuracy is not achieved for any combination of basis set and method on average. For OVGF, the smallest MAE is obtained with the aug-cc-pVQZ basis set, followed by cc-pVQZ and then by 6-311+G(2df,p), indicating that the inclusion of diffuse functions improves the description of the ionized state. Higher errors are obtained for the P3+ method, with the same order of the basis sets as for OVGF.

Figure 3 shows the deviations of the calculated ionization energy with respect to the experimental values for the six atoms for the OVGF and P3+ methods. For carbon atom, the OVGF overestimates and P3+ underestimates the experimental value. For nitrogen, the OVGF deviations are positive, while the P3+ underestimates the IE. None of the calculated values achieves chemical accuracy for C and N atoms. For oxygen, OVGF/aug-cc-pVQZ yields a deviation of only 0.0083 eV, within chemical accuracy, while the other two basis sets give larger negative deviations. Moreover, the P3+ method strongly underestimates the IE for oxygen, with deviations from the experiment in the range of -0.5191 eV to -0.3431 eV. For fluorine, OVGF/aug-cc-pVQZ achieves chemical accuracy with a deviation of -0.0346, while the P3+ underestimates the IE. For phosphorus, the OVGF deviations are below 0.0434 eV for all three basis sets, reaching chemical accuracy, while P3+ is close to the limit only with the aug-cc-pVQZ basis set. For sulfur, the OVGF underestimates the IE, while the P3+ method yields large deviations. In all cases, the 6-311+G(2df,p) basis set yields the largest deviations, while aug-cc-pVQZ gives the smallest, indicating that the inclusion of diffuse functions improves the description of the ionized state. The P3+ method systematically underestimates the experimental ionization energies for all atoms and basis sets, as also reflected in the negative MSE values reported in Table 2.

3.2. Ionization Energies of the H, C, N, O, F, P, and S Atoms Obtained with the Energy-Difference Approach and IP-EOM-CCSD

Since the OVGF and P3+ methods do not achieve chemical accuracy for the considered atoms (Table 1), higher level coupled-cluster calculations were performed using CCSD(T), CCSDT, and IP-EOM-CCSD methods as implemented in ORCA 6.1.0, combined with the basis set extrapolation to the complete basis set (CBS) limit and scalar relativistic corrections through the DKH2 Hamiltonian. In addition, a series of calculations were carried out with Gaussian 16 using B2PLYP, MP2, CCSD(T), QCISD(T), G3, and CBS-QB3 methods with the default frozen-core approximation, in which only the valence electrons are included in the correlation treatment.

Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 summarize the ionization energies of the H, C, N, O, F, P, and S atoms. The upper part of each table reports the results obtained with ORCA 6.1.0 using the CCSD(T), CCSDT, and IP-EOM-CCSD methods with three families of basis sets (aug-cc-pVXZ-DK, aug-cc-pVXZ, and ANO-RCC) together with the corresponding CBS extrapolations. For CCSD(T) and CCSDT the ionization energies were obtained as total energy differences between the neutral and cationic species, while for the IP-EOM-CCSD method they were obtained directly from the ionized roots of the neutral reference state. The lower part of the tables shows the results obtained with Gaussian 16 using the frozen-core default approximation and the energy-difference approach for the B2PLYP, MP2, CCSD(T), and QCISD(T) methods combined with aug-cc-pVQZ basis set, as well as the composite methods G3 and CBS-QB3. All Gaussian calculations employed the default frozen-core approximation. All deviations (ΔIE) are reported relative to the experimental values from NIST Atomic Spectra Database (ver. 5.12) [46]. In addition, the tables provide the spin multiplicities of both the neutral and cationic species. For each atom, the corresponding deviations are also presented graphically in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.

For the hydrogen atom (Table 3), which is a single-electron system, there is no electron correlation, and the ionization energy corresponds to the Hartree-Fock (HF) value. HF calculations were performed with three families of basis sets: the aug-cc-pVXZ (X = Q, 5), the relativistic aug-cc-pVXZ-DK (X = Q, 5) including the DKH2 Hamiltonian, and the ANO-RCC basis sets at the TZP, QZP, and Full contraction levels. All computed IE values are within the range of 13.6043 – 13.6057 eV, with deviations ranging from 0.0058 to 0.0073 eV relative to the experimental value of 13.5984 eV, which are clearly within chemical accuracy. The effect of scalar relativistic corrections (DKH2) is minimal for H, corresponding to only 0.0002 eV, as seen from the comparison between the values obtained with aug-cc-pVQZ-DK and aug-cc-pVQZ. The ANO-RCC basis sets already converge at the TZP contraction level (13.6051 eV), QZP and Full giving identical values (13.6054 eV). Moreover, the G3 and CBS-QB3 yield 13.6330 eV (ΔIE = 0.0345 eV) and 13.6007 eV (ΔIE = 0.0023 eV), respectively, CBS-QB3 providing the closest agreement with the experiment among all methods considered for this atom.

These results are also presented in Figure 4, which shows that IE values for the hydrogen atom are well within the chemical accuracy interval, with only minor differences between the basis sets.

In the case of carbon atom (Table 4), the CCSD(T)/aug-cc-pVXZ-DK ionization energy values converge from 11.2312 eV (QZ, ΔIE = -0.0291 eV) to 11.2476 eV (5Z, ΔIE = -0.0127 eV), while the CBS extrapolation yields 11.2635 eV (ΔIE = 0.0032 eV), remaining within the limits of chemical accuracy. CCSDT produces a CBS value of 11.2671 eV (ΔIE = 0.0068 eV), with a difference of only 0.0036 eV relative to CCSD(T), showing that the perturbative triples correction in CCSD(T) is enough for an accurate description of the ionization energy. For IP-EOM-CCSD, the first ionization is attributed to the α electron removal. This method already overestimates the IE at the QZ level with 0.0325 eV, and the CBS extrapolation further increases this deviation to 0.0778 eV, outside chemical accuracy, due to the missing triple excitations.

The aug-cc-pVXZ results are slightly higher, with the CCSD(T) CBS value at 11.2672 eV (ΔIE = 0.0069 eV) and CCSDT at 11.2708 eV (ΔIE = 0.0105 eV), both within chemical accuracy. The non-relativistic IP-EOM-CCSD CBS value of 11.3419 eV (ΔIE = 0.0816 eV) follows the same overestimation pattern. The DKH2 contribution is approximately -0.0037 eV across all methods and basis sets.

With ANO-RCC basis sets, the CCSD(T) values increase from 11.1962 eV (TZP) to 11.2163 eV (QZP) and 11.2274 eV (Full), the last value being within chemical accuracy (Full, ΔIE = -0.0329 eV). The CBS(T,Q) extrapolation yields 11.2337 eV (ΔIE = -0.0266 eV), also within chemical accuracy but further from experiment than the aug-cc CBS result. CCSDT CBS(T,Q) extrapolation value is 11.2384 eV (ΔIE = -0.0219 eV), while IP-EOM-CCSD CBS(T,Q) reaches 11.3029 eV (ΔIE = 0.0426 eV), at the chemical accuracy limit.

Among the Gaussian 16 results with the aug-cc-pVQZ basis set, B2PLYP overestimates the IE with 0.0893 eV, while MP2 gives a deviation of 0.0198 eV. QCISD(T) yields a deviation of -0.0531 eV. The results of the composite methods G3 and CBS-QB3 are 11.2114 eV (ΔIE = -0.0489 eV) and 11.1925 eV (ΔIE = -0.0678 eV), respectively.

A graphical summary of the deviations from the experimental value corresponding to the C atom is depicted in Figure 5. The four panels clearly show the good agreement with the experiment obtained for CCSD(T) and CCSDT, particularly with the CBS extrapolation, the systematic overestimation by the IP-EOM-CCSD method with the aug-cc basis sets, and the large negative deviations observed with the frozen-core methods.

For the nitrogen atom (Table 5), the CCSD(T)/aug-cc-pVXZ-DK results show a particularly fast convergence. At the 5Z level, the ionization energy of 14.5346 eV (ΔIE = 0.0005 eV), is in excellent agreement with the experimental value of 14.5341 eV. The CBS value of 14.5506 eV (ΔIE = 0.0165 eV) slightly overestimates the IE_exp value, but it remains within chemical accuracy. CCSDT yields identical results, with a CBS value of 14.5502 eV (ΔIE = 0.0161 eV) and a difference of only 0.0004 eV relative to CCSD(T). For IP-EOM-CCSD, the first ionization corresponds to the lowest α electron removal root. As observed for carbon, the method already overestimates the IE at the QZ level (ΔIE = 0.0100 eV), and the CBS extrapolation increases this deviation to 0.0588 eV, outside chemical accuracy. The aug-cc-pVXZ calculations follow the same pattern, with CCSD(T) CBS at 14.5567 eV (ΔIE = 0.0226 eV) and CCSDT CBS at 14.5563 eV (ΔIE = 0.0222 eV), both within chemical accuracy. The IP-EOM-CCSD CBS value of 14.5990 eV (ΔIE = 0.0649 eV) remains outside the threshold. The DKH2 contribution is -0.0061 eV, larger than for carbon (-0.0037 eV), consistent with the higher nuclear charge.

With the ANO-RCC basis sets, the CCSD(T) values increase from 14.4838 eV (TZP, ΔIE = -0.0504 eV) to 14.4976 eV (QZP, ΔIE = -0.0365 eV) and 14.5173 eV (Full, ΔIE = -0.0168 eV), the last two being within chemical accuracy. The CBS(T,Q) extrapolation yields 14.5114 eV (ΔIE = -0.0227 eV), also within chemical accuracy. CCSDT and IP-EOM-CCSD follow the same pattern, with ANO-RCC-Full values of 14.5176 eV (ΔIE = -0.0166 eV) and 14.5487 eV (ΔIE = 0.0146 eV), respectively.

Furthermore, B2PLYP exhibits a deviation of -0.0171 eV, lower than 0.0434 eV, while MP2 overestimates the value of IE with 0.0684 eV. The QCISD(T) yields 14.5001 eV (ΔIE = -0.0340 eV), whereas the G3 and CBS-QB3 exhibit deviations of -0.0273 eV and -0.0415 eV, respectively, all of them being within chemical accuracy.

Figure 6 presents the comparison of the ΔIE values for the nitrogen atom. The compact range of the CCSD(T) and CCSDT results around zero confirms the fast basis set convergence observed for this atom, while the IP-EOM-CCSD panel shows the progressive overestimation with increasing the basis set size.

In the case of the oxygen atom (Table 6), the CCSD(T) method with the aug-cc-pVXZ-DK basis sets yields IE values that converge systematically from 13.5204 eV at the QZ level (ΔIE = -0.0977 eV) to 13.5620 eV at the 5Z level (ΔIE = -0.0560 eV). The CBS value obtained with the SCF(5Z) + CBS(Q,5)-extrapolated correlation energy scheme is 13.6066 eV (ΔIE = -0.0115 eV), achieving chemical accuracy. The CCSDT results follow the same convergence pattern, with the CBS value of 13.6094 eV (ΔIE = -0.0087 eV), confirming that the perturbative triples in CCSD(T) is enough for an accurate description of the ionization energy, since the difference between the two methods at the CBS limit is only 0.0028 eV. For the IP-EOM-CCSD, the first ionization corresponds to the lowest β electron removal. The CBS value of 13.6483 eV (ΔIE = 0.0302 eV) remains within chemical accuracy, although the positive deviation indicates a slight overestimation of the experimental value. The best agreement for this method is obtained at the 5Z level (ΔIE = -0.0226 eV), where the basis set incompleteness error partially compensates for the overestimation associated with the absence of triple excitations.

The non-relativistic calculations with the aug-cc-pVXZ basis sets show a similar convergence but with higher ionization energy values. For CCSD(T), the CBS value is 13.6155 eV (ΔIE = -0.0026 eV), closer to the experiment than the DKH2 result. For CCSDT, the CBS value of 13.6183 eV (ΔIE = 0.0002 eV) is in excellent agreement with the experimental value. The IP-EOM-CCSD CBS value of 13.6571 eV (ΔIE = 0.0390 eV) follows the same trend and remains under the 0.0434 eV value. The contribution of DKH2 scalar relativistic correction to the IE of oxygen is constant, with a value of -0.0090 eV across all methods and basis sets.

The ANO-RCC basis sets show a significantly slower convergence. For CCSD(T), the ionization energy increases from 13.3830 eV (TZP, ΔIE = -0.2350 eV) to 13.4440 eV (QZP, ΔIE = -0.1741 eV) and 13.5149 eV (Full, ΔIE = -0.1031 eV), although none of them achieve chemical accuracy. The CBS(T,Q) extrapolation yields 13.4996 eV (ΔIE = -0.1185 eV), which is worse than the uncontracted Full result, indicating the absence of diffuse functions in the ANO-RCC family. The same trend is observed for CCSDT and IP-EOM-CCSD, with ANO-RCC-Full values of 13.5183 eV (ΔIE = -0.0998 eV) and 13.5454 eV (ΔIE = -0.0726 eV), respectively.

For the results obtained with Gaussian 16 using the aug-cc-pVQZ basis set, B2PLYP overestimates the IE with 0.1965 eV, while MP2 underestimates it with 0.1570 eV. The QCISD(T) method yields a deviation from the experimental value of -0.1032 eV. None of these methods achieves chemical accuracy. Moreover, the G3 and CBS-QB3 methods yield values of 13.5478 eV (ΔIE = -0.0702 eV) and 13.5869 eV (ΔIE = -0.0311 eV), respectively. Among these composite methods, only CBS-QB3 achieves chemical accuracy.

The corresponding deviations are illustrated in Figure 7, where the convergence from negative ΔIE values at the QZ and 5Z levels toward the experimental value at the CBS limit is visible for both CCSD(T) and CCSDT. The large underestimation with the ANO-RCC basis sets is also visible.

For the fluorine atom (Table 7), the CCSD(T)/aug-cc-pVXZ-DK converges from 17.3488 eV (QZ, ΔIE = -0.0740 eV) to 17.3863 eV (5Z, ΔIE = -0.0365 eV), and the CBS value of 17.4283 eV (ΔIE = 0.0055 eV) achieves chemical accuracy. CCSDT produces a CBS value of 17.4275 eV (ΔIE = 0.0047 eV), following the same trend with CCSD(T). For IP-EOM-CCSD, the first ionization is attributed to the lowest β electron removal root. Unlike the other atoms where IP-EOM-CCSD overestimates IE, here both QZ (ΔIE = -0.1017 eV) and 5Z (ΔIE = -0.0498 eV) underestimate it, and the CBS extrapolation yields a result of 17.4276 eV (ΔIE = 0.0047 eV), in excellent agreement with the experiment and within chemical accuracy.

The aug-cc-pVXZ results have higher ionization energy values. The CCSD(T) CBS is 17.4410 eV (ΔIE = 0.0182 eV) and the CCSDT CBS is 17.4403 eV (ΔIE = 0.0175 eV), both within chemical accuracy. The IP-EOM-CCSD CBS value of 17.4402 eV (ΔIE = 0.0174 eV) also remains within the threshold. The DKH2 contribution is -0.0127 eV, the largest among the period 2 atoms, consistent with fluorine having the highest nuclear charge in this group.

The ANO-RCC basis sets show a slow convergence also for fluorine. CCSD(T) values range from 17.1979 eV (TZP, ΔIE = -0.2249 eV) to 17.2638 eV (QZP, ΔIE = -0.1591 eV) and 17.3443 eV (Full, ΔIE = -0.0785 eV), none achieving chemical accuracy. The CBS(T,Q) extrapolation yields 17.3251 eV (ΔIE = -0.0977 eV), again worse than the Full result. Moreover, the CCSDT values follow the same trend, underestimating the ionization energy with 0.2238 eV (TZP), 0.1586 (QZP), 0.0787 (Full), and 0.0977 (CBS(T,Q)). For IP-EOM-CCSD, the ANO-RCC-Full and CBS (T,Q) deviations of -0.0987 eV and -0.1309 eV are both far from chemical accuracy.

Furthermore, B2PLYP strongly overestimates IE with 1.4022 eV, the largest error observed for any method across all atoms studied. MP2 gives a deviation of -0.0270 eV, within chemical accuracy. QCISD(T) yields a value of 17.3486 eV, while G3 and CBS-QB3 methods exhibit deviations of -0.0344 eV and 0.0513 eV, respectively, only the G3 result being within chemical accuracy.

As shown in Figure 8, for the fluorine atom, the CCSD(T), CCSDT, and IP-EOM-CCSD values improve with the basis set, and the CBS extrapolations achieve chemical accuracy. The anomalous B2PLYP overestimation is clearly visible in the frozen-core panel. For clarity, the vertical axis was limited to 0.3 eV, even though the B2PLYP deviation is significantly higher, at 1.4022 eV.

In the case of the phosphorus atom (Table 8), the CCSD(T)/aug-cc-pVXZ-DK ionization energy of 10.4796 eV at the QZ level (ΔIE = -0.0071 eV) is already very close to the experimental value of 10.4867 eV. At the 5Z level, the value increases to 10.5003 eV (ΔIE = 0.0136 eV), and the CBS extrapolation further raises the overestimation to 10.5242 eV (ΔIE = 0.0375 eV), approaching the chemical accuracy limit. The CCSDT exhibits the same behavior, with a CBS value of 10.5245 eV (ΔIE = 0.0378 eV) and a difference of only 0.0003 eV relative to CCSD(T). For IP-EOM-CCSD, the first ionization corresponds to the lowest root associated with the α electron removal. The method overestimates the IE already at the QZ level (ΔIE = 0.0192 eV), and the CBS value of 10.5817 eV (ΔIE = 0.0951 eV) represents the largest IP-EOM-CCSD deviation among all atoms studied with the aug-cc-pVXZ-DK basis sets.

The aug-cc-pVXZ calculations yield systematically higher values. The CCSD(T) CBS of 10.5365 eV (ΔIE = 0.0498 eV) exceeds the chemical accuracy threshold, in contrast to the DKH2 result, which remains within the limit. This makes phosphorus the only atom for which the inclusion of DKH2 relativistic corrections improves the CCSD(T) CBS result from outside to within the chemical accuracy. The CCSDT CBS value of 10.5368 eV (ΔIE = 0.0501 eV) exhibits the same behavior. The IP-EOM-CCSD CBS value of 10.5939 eV (ΔIE = 0.1072 eV) confirms the strong overestimation pattern. The DKH2 contribution is -0.0121 eV, comparable to that of fluorine (-0.0127 eV), indicating that the scalar relativistic correction remains significant for the period 3 atom.

With the ANO-RCC basis sets, the CCSD(T) energy increase from 10.4497 eV (TZP, ΔIE = -0.0370 eV) to 10.4594 eV (QZP, ΔIE = -0.0273 eV) and 10.4890 eV (Full, ΔIE = 0.0023 eV), the latter being in excellent agreement with experiment. The CBS(T,Q) extrapolation yields 10.4672 eV (ΔIE = -0.0195 eV), within chemical accuracy but further from experiment than the Full value. CCSDT with ANO-RCC-Full basis set yields 10.4909 eV (ΔIE = 0.0043 eV), also within chemical accuracy. For IP-EOM-CCSD with the ANO-RCC basis sets, the ionization energy increases from 10.4708 eV (TZP, ΔIE = -0.0159 eV) to 10.4865 eV (QZP, ΔIE = -0.0002 eV), the latter being in excellent agreement with the experiment. However, the Full value of 10.5337 eV (ΔIE = 0.0471 eV) overestimates the experimental value and slightly exceeds the chemical accuracy. The CBS(T,Q) yields a value of 10.4980 eV (ΔIE = 0.0113 eV), within chemical accuracy.

For the Gaussian 16 results, B2PLYP underestimates the ionization energy with 0.1445 eV and MP2 yields a deviation of only 0.0029 eV, the closest to experiment among all Gaussian methods with the aug-cc-pVQZ basis set. QCISD(T) result overestimates the value with 0.0125 eV, whereas the G3 and CBS-QB3 yield deviations of -0.0230 eV and -0.0394 eV, respectively, both within chemical accuracy.

The performance of the tested methods for the phosphorus atom is illustrated in Figure 9. The figure highlights that the QZ values are already close to the experiment and the CBS extrapolation leads to a slight overestimation for CCSD(T) and CCSDT, while the IP-EOM-CCSD CBS values exceed the chemical accuracy limit.

For the sulfur atom (Table 9), the IE results obtained with the CCSD(T) and the aug-cc-pVXZ-DK basis sets show the largest basis set dependence among all atoms studied. The QZ value of 10.2414 eV (ΔIE = -0.1186 eV) and the 5Z value of 10.2842 (ΔIE = -0.0758 eV) significantly underestimate the experimental value of 10.3600 eV, while the CBS value of 10.3297 eV (ΔIE = -0.0304 eV) achieves chemical accuracy. CCSDT yields a CBS value of 10.3332 eV (ΔIE = -0.0268 eV), with a difference of 0.0035 eV relative to CCSD(T). For IP-EOM-CCSD, the first ionization corresponds to the lowest root associated with the β electron removal. The QZ value (ΔIE = -0.0415 eV) underestimates the IE, while the 5Z value (ΔIE = 0.0169 eV) shifts to the opposite side. The CBS extrapolation increases the overestimation to 0.0782 eV, outside chemical accuracy, a behavior similar to that observed for oxygen but with a larger final deviation.

The aug-cc-pVXZ results yield higher ionization energy values, with the CCSD(T) CBS result being 10.3441 eV (ΔIE = -0.0159 eV), closer to the experiment than the DKH2 result. The CCSDT CBS is 10.3476 eV (ΔIE = -0.0124 eV), also within chemical accuracy and providing the best overall agreement with the experiment among the ORCA methods. The IP-EOM-CCSD CBS of 10.4530 eV (ΔIE = 0.0930 eV) remains far from the threshold. The DKH2 contribution is -0.0142 eV, the largest among all atoms studied, consistent with sulfur having the highest nuclear charge.

The ANO-RCC basis sets show the weakest convergence for sulfur compared to all other atoms. The CCSD(T) values range from 10.1113 eV (TZP, ΔIE = -0.2487 eV) to 10.1488 eV (QZP, ΔIE = -0.2112 eV) and 10.2542 eV (Full, ΔIE = -0.1058 eV), all being far from chemical accuracy. The CBS(T,Q) extrapolation yields 10.1778 eV (ΔIE = -0.1823 eV), significantly worse than the Full result. CCSDT follows the same pattern, with ANO-RCC-Full at 10.2581 eV (ΔIE = -0.1019 eV) and CBS(T,Q) at 10.1778 eV (ΔIE = -0.1823 eV). For IP-EOM-CCSD, the ANO-RCC-Full value of 10.3493 eV (ΔIE = -0.0107 eV) is within chemical accuracy, although this result likely reflects the cancelation of errors between the missing triple excitations and basis set incompleteness.

Finally, the B2PLYP computation yields a deviation of only -0.0047 eV, the closest to the experimental value for this atom, while the MP2 underestimates the ionization energy with 0.2204 eV, the largest MP2 deviation among all atoms studied. QCISD(T) yields a deviation of -0.1051 eV, while G3 and CBS-QB3 underestimate the experimental value with 0.0909 eV and 0.1156 eV, respectively. Thus, sulfur and carbon are the only atoms for which none of the Gaussian composite methods achieve chemical accuracy.

Figure 10 shows the ΔIE values for the sulfur atom. The ANO-RCC basis sets show significant negative deviations and IP-EOM-CCSD goes from negative to positive ΔIE from the QZ to 5Z levels.

It should be noted that the CCSD(T) results obtained with Gaussian 16 using the aug-cc-pVQZ basis set are not directly comparable to ORCA CCSD(T)/aug-cc-pVQZ values, as in the Gaussian calculations were performed with the default frozen-core approximation, while all ORCA computations employed NoFrozenCore. The resulting differences range from 0.002 eV for sulfur to 0.093 eV for carbon, depending on the relative contribution of core-valence correlation to the ionization energy of each atom. Moreover, the same frozen-core approximation was also used in all Gaussian calculations considered here.

The CCSDT calculations are considerably more demanding than CCSD(T), when the same computational settings are used. For the fluorine atom with the aug-cc-pV5Z-DK basis set, employing the same setup, the CCSDT calculation took approximately 10 hours, while the corresponding CCSD(T) calculation required only 8 minutes and 27 seconds. This difference shows the significantly higher cost associated with the iterative treatment of triple excitations in CCSDT.

Table 10 summarizes the energies, types and degeneracies of the HOMO and HOMO – 1 orbitals for the C, N, O, F, P, and S atoms, together with the HOMO – HOMO-1 energy gap, calculated at QCISD(T)/aug-cc-pVQZ level of theory. The spin character of HOMO confirms the behavior in the IP–EOM–CCSD calculations, indicating that for C, N, and P, the lowest ionization energy is associated with the α electron removal, while for O, F, and S, it comes from the β electron removal. The smallest HOMO – HOMO-1 gap is obtained for sulfur (0.98 eV), while the largest gap corresponds to the phosphorus (4.47 eV).

The overall performance of the computational methods is summarized in Table 11, which reports the mean absolute error (MAE) and mean signed error (MSE) for each method and basis set combination, averaged over the C, N, O, F, P, and S atoms.

Among all combinations of method and basis set, the lowest MAE values are obtained with the CCSDT CBS extrapolation using the aug-cc-pVXZ-DK basis sets (MAE = 0.0168 eV, MSE = 0.0050 eV), closely followed by CCSD(T) CBS extrapolated correlation (Q,5) with the same basis sets (MAE = 0.0174 eV, MSE = 0.0035 eV). The nearly identical MAE values confirm that the perturbative triples treatment in CCSD(T) is enough for all atoms studied. The small positive MSE indicates a slight overestimation at the CBS limit.

The aug-cc-pVXZ CBS results show slightly larger MAE values of 0.0194 eV for CCSD(T) and 0.0188 eV for CCSDT, with MSE values of 0.0132 eV and 0.0147 eV, respectively. The increase in both MAE and MSE relative to the DKH2 results reflects the absence of scalar relativistic corrections, which systematically raise the computed ionization energies. Nevertheless, the CBS results obtained with these basis sets remain within chemical accuracy on average.

For IP-EOM-CCSD, the CBS extrapolation yields MAE = 0.0575 eV (MSE= 0.0575 eV) with DKH2 and MAE = 0.0672 (MSE = 0.0672 eV) without DKH2. In both cases, MAE equals MSE, indicating that the method systematically overestimates the ionization energy for all atoms at the CBS limit. At the finite basis set level, the performance of IP-EOM-CCSD is comparable to CCSD(T), with MAE values of 0.0463 eV (QZ-DK) and 0.0390 eV (5Z-DK), suggesting that the basis set incompleteness error partially compensates for the missing triple excitations.

The ANO-RCC basis sets have larger MAE values, ranging from 0.1433 eV (TZP) to 0.0566 eV (Full) for CCSD(T). The MSE values are consistently negative, indicating an underestimation due to the absence of diffuse functions. The CBS(T,Q) extrapolation (MAE = 0.0784 eV) does not improve over the Full result, consistent with the observation that the extrapolation cannot compensate for the missing diffuse functions.

Across the Gaussian 16 methods that were used with the aug-cc-pVQZ basis set, B2PLYP exhibits the largest MAE (0.3091 eV) with a positive MSE (0.2536), mostly influenced by the overestimation for fluorine. The composite methods G3 (MAE = 0.0491 eV, MSE = -0.0491 eV) and CBS-QB3 (MAE = 0.0578 eV, MSE = -0.0407 eV) provide the best performance among these methods, although none of them achieves chemical accuracy on average. For G3, the equality between MAE and MSE indicates a systematic underestimation for all six atoms. QCISD(T) (MAE = 0.0637 eV), MP2 (MAE = 0.0826 eV), and CCSD(T) (MAE = 0.0756 eV) with the aug-cc-pVQZ basis set perform somewhat worse. For comparison, the corresponding ORCA CCSD(T)/aug-cc-pVQZ result gives a lower MAE of 0.0490 eV (MSE = -0.0474 eV). The difference between the Gaussian and ORCA CCSD(T) values can be attributed mainly to the frozen-core approximation employed in Gaussian calculations, while the ORCA results were obtained with NoFrozenCore.

Overall, the statistical performance is summarized graphically in Figure 11, which reports the MAE values averaged over the C, N, O, F, P, and S atoms. The CCSD(T) and CCSDT CBS extrapolations with aug-cc-pVXZ-DK have the lowest MAE values, while the ANO-RCC and IP-EOM-CCSD results exhibit consistently higher error. Among the methods with the frozen-core approximation, G3 and CBS-QB3 provide the best performance, however none of them achieve chemical accuracy on average.

4. Conclusions

In this work, the first ionization energies of the H, C, N, O, F, P, and S atoms were investigated using several quantum chemical methods, with the aim of identifying the computational protocols able to reproduce the experimental values within chemical accuracy. The methods included the electron propagator approaches OVGF and P3+, as implemented in Gaussian 16, the coupled-cluster methods CCSD(T), CCSDT, and IP-EOM-CCSD, as implemented in ORCA 6.1.0, and the composite methods G3 and CBS-QB3. Several post-Hartree-Fock and DFT methods were also tested using the energy-difference approach.

The OVGF and P3+ methods did not reach chemical accuracy on average for the six atoms (C, N, O, F, P, S). For all basis sets, OVGF performed better than P3+, while P3+ constantly underestimated the experimental ionization energies, as shown by the negative MSE. Among the tested basis sets, aug-cc-pVQZ gave the best results for both methods, followed by cc-pVQZ and 6-311+G(2df,p). Chemical accuracy was achieved only in some cases, specifically with OVGF for phosphorus for all three basis sets, and for oxygen and fluorine when the aug-cc-pVQZ basis set was used. The best agreement obtained with the propagator methods was found for P atom, which has a large HOMO – HOMO-1 gap.

The CCSD(T) and CCSDT methods, used together with the aug-cc-pVXZ-DK basis sets and the CBS extrapolation scheme provided the most accurate result. The corresponding MAE values were 0.0174 eV for CCSD(T) and 0.0168 eV for CCSDT. The difference between the two methods at the CBS limit was small for each atom, showing that the perturbative triples correction in CCSD(T) is enough for an accurate description of these ionization energies. Therefore, the much more computationally expensive iterative treatment of triple excitations in CCSDT did not bring significant improvement for the ionization energies of the atoms that were studied. Moreover, considering the much higher computational cost of CCSDT, CCSD(T) combined with CBS extrapolation represents the best compromise between the accuracy and the hardware resources.

The scalar relativistic contribution introduced through the DKH2 Hamiltonian was small, but consistent. Its size increased with the atomic number, from about 0.0002 eV for hydrogen to approximately 0.0142 eV for sulfur. For hydrogen, which is a one-electron system, the ionization energy is described at the Hartree-Fock level. All HF values obtained with the tested basis set families remained within chemical accuracy, and the CBS-QB3 gave the closest agreement with the experimental value.

The IP-EOM-CCSD method showed a different behavior. Since this method calculates the ionization energies directly from the ionized roots of the neutral reference state, it provides an independent comparison with the energy difference results. At the finite basis sets, IP-EOM-CCSD often gave values close to those obtained with CCSD(T), due to the basis set incompleteness error, which compensated for the missing triple excitations. At the CBS limit, however, this effect was reduced, and the method began to overestimate the experimental values.

The ANO-RCC basis sets showed a less uniform convergence than the aug-cc-pVXZ-DK basis sets. The CBS(T,Q) extrapolation based on the TZP and QZP contractions did not always improve the results, and in some cases, was worse than the Full result. This suggests that ANO-RCC basis sets may not be ideal for the calculations of these atomic ionization energies.

The Gaussian 16 results obtained with the default frozen-core approximation also showed that core-valence correlation is important for achieving chemical accuracy. Among these methods, G3 and CBS-QB3 gave the best overall agreement with the experiment, but none of them reaches chemical accuracy on average over the six atoms. The B2PLYP functional shows a particularly large overestimation for fluorine, which influenced its overall MAE. The comparison between Gaussian CCSD(T)/aug-cc-pVQZ and the corresponding all-electron ORCA result showed that the differences can be significant due to the frozen-core approximation. This confirms that the treatment of core-valence correlation must be considered carefully when ionization energies are compared at this level of accuracy.

Overall, the results show that the most robust approach for the atomic ionization energies studied here is CCSD(T) with the aug-cc-pVXZ-DK basis sets, all-electron correlation, DKH2 scalar relativistic corrections, and the CBS(Q,5) extrapolation scheme. CCSDT confirms the accuracy of the perturbative triples treatment in CCSD(T), while IP-EOM-CCSD provides a useful direct comparison. Thus, the consistency of the results obtained in this study supports the use of this protocol as a reference for the ionization energies of the H, C, N, O, F, P, and S atoms found in proteins and amino acids, as well as in the related benchmark studies.