Exchange-Correlation Functional Comparison of Electronic Energies in Atoms Using a Grid Basis ()
1. Introduction
Density functional theory (DFT) is a quantum mechanical method widely used in chemistry [1] and materials science [2] to calculate system properties from first principles within its one significant approximation: the exchange-correlation functional. DFT calculations of atomic total electronic energies are important for calibrating advancements in exchange-correlation functionals among codes [3] and producing reliable pseudopotentials [4]. Recent advances in exchange-correlation functionals continue to show reduced error in molecular [5] [6] and solid-state [6] [7] test sets compared and call for an investigation of the functionals' effect on energies calculated in isolated atoms.
Non-empirical pseudopotentials constructed for use with plane waves in DFT simulations of systems with periodic boundary conditions rely on atomic calculations in which the electrons are represented through Kohn-Sham states constituting the electron density on a radial grid [8]. Existing research has explored the effects of exchange-correlation functional at the local-density [8] and generalized-gradient [9] [10] levels of approximation, including the effects of incorporating spin-polarization and relativistic effects [8] on atomic density functional total energies and ionization energies.
In this work, we compare the published values of these total and (first) ionization energies of all-electron atoms to two existing codes in order to establish a baseline of comparison. We then extend these studies to the meta-generalized gradient and exact-exchange and hybrid levels of exchange-correlation approximation, investigating the effects of these additional levels of complexity on the total energies and ionization energies of atoms calculated using a grid basis, comparing ionization energy to experiment. Finally, we also investigate the effect of combining relativistic effects in both scalar-relativistic and fully relativistic calculations to these energies using the same levels of exchange-correlation approximation.
2. Method
Quantum ESPRESSO [11] [12] is “an integrated suite of Open-Source computer codes for electronic-structure calculations and materials modeling at the nanoscale, based on density-functional theory, plane waves, and pseudopotentials” that includes an atomic density-functional code atomic OPIUM [13] is a stand- alone pseudopotential generation code that implements both density functional and Hartree-Fock theory to calculate electronic total energies. Both atomic and OPIUM use a radial grid to represent the single-particle Kohn-Sham states that is of the form:
(1)
where the parameters a and b determine the relationship of the N points that are indexed from
to
. In atomic, a by default is set to 1/Z where Z is the atomic number of the atom being simulated, and, in OPIUM, a defaults to
. In both codes, b is an exponential grid spacing parameter.
As of Quantum ESPRESSO version 6.3, atomic includes support for a variety of exchange-correlation functionals at the local-density approximation (LDA) and generalized-gradient approximation (GGA) levels while, at version 4.1, OPIUM implements exchange-correlation functionals at the LDA and GGA levels as well as the Hartree-Fock [14] (HF) and hybrid (PBE0) [15] functional levels. Yao and Kanai [16] implemented the TPSS [17] and SCAN [18] meta-GGA exchange-correlation functionals in a modified version of atomic. In this investigation, electronic orbital occupation of the Kohn-Sham states was chosen to match the ground-state electron configurations of the atoms established by experiment [19].
3. Results
3.1. Total Energy
To corroborate the existing implementations of exchange-correlation in atomic and OPIUM against the published values, Table 1 presents the mean absolute relative error of the DFT total energies produced by the atomic and OPIUM codes from the published values for LDA/LSD (local spin density) [8] [20] and LDA & GGA [22] exchange-correlation functionals. The atomic data use the internal exchange-correlation functional implementation in Quantum ESPRESSO, and tests using the libxc [27] implementations of exchange-correlation in Quantum ESPRESSO produce statistically indistinguishable differences. Table 2 shows the mean absolute relative error of the HF total energies produced by the OPIUM code from prior published values [28] [29].
The NIST data are reported to 10−6 Ha so a difference of 9 × 10−7 Ha for the smallest total energy (E(H) [LDA] = −0.445671 Ha) would produce the largest possible absolute relative error of 2 × 10−6. The calculated LDA MAREs for both atomic and OPIUM are below this threshold, and the atomic LSD MARE also
Table 1. Mean absolute relative error (MARE) for the data in each reference for atomic code in Quantum ESPRESSO [11] [12] and the OPIUM code [13], for varying exchange-correlation functionals and spin-polarization treatment. The NIST [8] [19] [20] data use the Vosko-Wilk-Nusair [21] (VWN) parameterization of the local density approximation (LDA) while the Lee and Martin [22] data use the Perdew-Zunger [23] parameterization for LDA. Lee and Martin compare two generalized-gradient approximations (GGAs), Perdew-Wang 1991 [24] (PW91) and Perdew-Burke-Ernzerhof [25] [26] (PBE).
Table 2. Mean error (ME), mean absolute error (MAE), mean relative error (MRE), and mean absolute relative error (MARE) for the total energy using the Hartree-Fock method in the OPIUM code with respect to the values reported in Davidson et al. [28] for elements Z = 3 - 10.
lies below this threshold (the spin-polarized VWN functional is not implemented in OPIUM).
The data of Lee and Martin are reported to 10−3 Ry so a difference of 9 × 10−4 Ry for the total energy of hydrogen would produce the largest possible absolute relative error of 1 × 10−3. The calculated LDA, PW91, and PBE MAREs for both atomic and OPIUM lie below this threshold.
Davidson et al. [2] report Hartree-Fock energies to 10−6 Ha calculated using Slater-type orbitals with an 11s, 10p, 9d, 8f, 6g, 4h, 2i basis set. A comparison with Hartree-Fock calculations using the exponential radial grid in OPIUM shows mean relative and mean absolute relative errors of 0.04%, greater than the 0.0001% difference that would be produced by 9 × 10−7 Ha difference on the smallest total energy in this set (E(Li) [HF] = -7.432727 Ha). The values for Z = 6 - 8 deviate in the 0.05% 0.2% range while the values for Z = 3 - 5, 9 - 10 deviate on the order of 1 × 10−6 %, implying that the only numerically significant difference lies in the total energies for C, N, and O.
For completeness, we present the total energies for elements Z = 1 - 92 in Tables A1-A7 of the Appendix calculated using the DFT exchange-correlation functionals LDA-PW, GGA-PBE, GGA-PBEsol, metaGGA-TPSS, metaGGA-SCAN, and hybrid-PBE0 together with Hartree-Fock.
Following Kotochigova et al., we plot −Z7/3Etotal against Z to contrast the DFT exchange-correlation functionals with Hartree-Fock (Figure 1). At the global level, all of these methods follow the same trend. Zooming in on the Z = 2 - 10 region, we see the LDA energies lowest followed by the GGA-PBEsol energies, and the GGA-PBE, metaGGA-TPSS, metaGGA-SCAN, and hybrid PBE0 energies converging at this scale for Z > 5.
In addition to the atomic number-scaled total energies, we plot the relative difference from the corrected Thomas-Fermi energy (Figure 2). For this quantity, the differences between the functionals appear most starkly in low Z elements. In Period 3, the metaGGA functions group with the largest relative difference, followed by GGA-PBE and hybrid-PBE0. The reduced gradient expansion coefficient in both the exchange and correlation components of GGA-PBEsol
Figure 1. Scaled total energy −Z−7/3Etotal with atomic number Z for all methods for all atoms (left) and Period 2 atoms (right).
compared to GGA-PBE makes its relative differences from the corrected Thomas-Fermi energy closer to the LDA-VWN values.
3.2. Ionization Energy
The (first) ionization energy of an atom is calculated from the difference between the total energies of the all-electron atom and the singly ionized cation [30]. The density functional theory values for ionization energy reproduce the qualitative trends of the experimentally determined values [8], and including gradient and higher-order corrections does not significantly alter these trends (see Figures 3-5). Ionization energies rise across each period of elements, starting low in Group 1 and rising to a peak in Group 18, before dropping for the Group 1 element in the subsequent period. For increasing Z across Groups 4 - 12, ionization energy shows fluctuations of up to nearly 5 eV reflecting the complexities of the d- and f-shell orbital energies in both the neutral atom and the singly ionized cation.
Figure 2. Relative difference of total energy from corrected Thomas-Fermi energy ΔE/ECTF with atomic number Z for all methods.
Figure 3. Ionization energy Ei with atomic number Z for two GGA functionals: PBE and PBEsol.
Figure 4. Ionization energy Ei with atomic number Z for two meta-GGA functionals: TPSS and SCAN.
Figure 5. Ionization energy Ei with atomic number Z for Hartree-Fock and the hybrid PBE0 functional.
3.2.1. Comparison to Experiment
For more detailed consideration of the calculated ionization energies, we investigate the difference of the calculated ionization energy from the experimental values with increasing atomic number (Figures 6-8) as reported in the NIST Computational Chemistry Comparison and Benchmark Database [31]. We collate into Figure 9 the comparison of the differences using all of these levels of theory. The mean-field approximation of density functional theory shows up with peak differences occurring at Z = 2, 27, 28, 78, 79 and Group 16 & 17 elements, the difference from experiment trending downward with increasing Z without any attention to relativistic effects in the calculation (we consider these effects in Section 3.2.2).
3.2.2. Relativistic Effects
The effects of special relativity on the electrons can be accounted for by transforming the Kohn-Sham equations into Dirac-like equations with a spinor containing two components, each indexed by the quantum numbers n,
, and j where j is the total angular momentum. These equations are coupled first-order equations of the spinor components containing the Dirac quantum number
, and the charge density is constructed by summing the sum of the squares of the
Figure 6. Difference in ionization energy from experiment ΔEi versus atomic number Z for two GGA functional.
Figure 7. Difference in ionization energy from experiment ΔEi versus atomic number Z for two meta-GGA functional.
Figure 8. Difference in ionization energy from experiment ΔEi versus atomic number Z for HF and the PBE0 functional.
Figure 9. Difference in ionization energy from experiment ΔEi versus atomic number Z for HF and the PBE, PBEsol, TPSS, SCAN, and PBE0 functionals.
Table 3. Mean absolute relative error (MARE) for the ionization energy with respect to the experimental values for elements Z = 1 - 36, varying relativistic treatment and exchange-correlation functional.
Table 4. Mean absolute relative error (MARE) for the ionization energy with respect to the experimental values for elements Z = 37 - 92, varying relativistic treatment and exchange-correlation functional.
spinor components, multiplied by the occupancies and divided by the spherical surface area, over the quantum numbers n,
, and j. This is the so-called fully relativistic treatment. A simplification reduces the coupled first-order equations into a single second-order equation for the large spinor component and averages over the spin-orbit components—this is the scalar relativistic case.
In Figures 10-12, we show the effects of the two models of relativity together with the non-relativistic calculations on the difference of calculated and experimental ionization energies for elements in periods four, five, and six (Z = 19 - 36, 37 - 54, and 55 - 86, respectively) using the GGA-PBE functional. The qualitative ordering of the three calculations on each particular element remain the same in LDA-PW and GGA-PBEsol as in GGA-PBE.
Figure 10. Difference in ionization energy ΔEi from experiment with atomic number Z for Period 4 elements (Z = 19 - 36) using the PBE exchange-correlation functional varying relativistic treatment.
Figure 11. Difference in ionization energy ΔEi from experiment with atomic number Z for Period 5 elements (Z = 37 - 54) using the PBE exchange-correlation functional varying relativistic treatment.
Figure 12. Difference in ionization energy ΔEi from experiment with atomic number Z for Period 6 elements (Z = 55 - 86) using the PBE exchange-correlation functional varying relativistic treatment.
4. Conclusion
We present density functional atomic total energies for an exponential radial grid basis using exchange-correlation approximations ranging from local density through generalized gradient and meta-generalized gradient to exact-exchange hybrid levels of complexity in comparison with the same energies calculated using the Hartree-Fock method. We demonstrate that, while these increased levels of complexity change the absolute energies and first ionization energies of the atoms, the relative trends with increasing atomic number and difference from experiment remain largely the same. Adding scalar and spinor relativistic corrections to the density-functional Hamiltonian can alter significantly the difference of the first ionization energy from experiment, changing the sign of the difference, but the mean absolute relative error in the first ionization energy for elements Z = 37 - 92 reduces only from 7% to 5%.
Appendix—Total Energies
Table A1. Total spin-polarized energies (in Ha) for atoms in the first two periods of the periodic table (Z = 1 - 10) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction. Values are depicted to 10−4, indicative of the digit where no grid basis error is evident.
Table A2. Total spin-polarized energies (in Ha) for atoms in the third period of the periodic table (Z = 11 - 18) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction. Values are depicted to 10−4, indicative of the digit where no grid basis error is evident.
Table A3. Total energies for atoms (in Ha) in the fourth period of the periodic table (Z = 19-36) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction.
Table A4. Total energies (in Ha) for atoms in the fifth period of the periodic table (Z = 37 - 54) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction.
Table A5. Total energies (in Ha) for atoms in the sixth period of the periodic table without the lanthanoids (Z = 55 - 56 & 72 - 86) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction.
Table A6. Total energies for lanthanoid atoms in the sixth period of the periodic table (Z = 57-71) with varying exchange-correlation functionals (as well as DFT to HF) and no relativistic correction.
Table A7. Total energies (in Ha) for atoms in the seventh period of the periodic table up to U (Z = 87-92) with varying exchange-correlation functionals and no relativistic correction.