Self-consistent self-interaction corrected density functional theory calculations for atoms using Fermi-Löwdin orbitals: Optimized Fermi-orbital descriptors for Li–Kr

In conventional semi-local density functional theories (DFTs), the classical Coulomb energy of the electrons, which includes the interaction of each electron with itself, is included exactly in the total energy expression, while the exchange-correlation energy contribution, including the self-exchange interaction, is approximated. It has been recognized for nearly four decades¹ that this inevitably results in an incomplete cancellation of electron self-interaction and leads to self-interaction errors (SIE) that can cause significant problems in DFT calculations. For example, the DFT potential for an electron in an atom tends to zero exponentially at large distances from the nucleus, instead of exhibiting the physically correct −1/r behavior. As a result, the energy of the highest occupied orbital (HOMO) is always too high in DFT, making the HOMO level a poor predictor of the first ionization energy. Perdew and Zunger introduced a self-interaction correction (SIC) that significantly reduces the SIE.¹ Their PZ-SIC approach removes self-interaction on an orbital-by-orbital basis, leading to an orbital-dependent total energy and a corresponding set of orbital-dependent Schrödinger-like equations instead of the single Kohn-Sham equation in DFT. Pederson and co-authors subsequently showed that, for a fully variational solution that minimizes the DFT-SIC total energy, the orbitals must satisfy additional conditions called the localization equations (LEs),² so-named because the optimal orbitals tend to be localized to give the most negative correction. Conceptually, the difference between the uncorrected and corrected theories can be stated as follows: in DFT, one needs the correct electron density to obtain the minimum energy, but in DFT-SIC, one needs not only the correct density but also the correct representation of that density in terms of a specific set of orbitals. Determining these orbitals by explicitly satisfying the LE presents a computational bottleneck that is sufficiently severe that it has kept the PZ-SIC method from being widely used,³ though important recent work has been done.⁴ 

Recently, Pederson, Perdew, and Ruzsinszky introduced an alternate implementation of the PZ-SIC method that sidesteps the LE.⁵ In the Fermi-Löwdin-orbital self-interaction correction (FLO-SIC), local orbitals are obtained through a well-defined unitary transformation of the occupied space that generates the Fermi orbitals.⁶ The ingredients required to define this transformation consist of a set of M Fermi-orbital descriptors (FODs), where M is the total number of electrons. Each FOD is a position in 3-dimensional space that can be thought of loosely as the classical position for the corresponding electron. Because different choices of the FODs result in different local orbitals, the DFT-SIC total energy depends on the FOD positions. Pederson derived the derivatives of the DFT-SIC total energy with respect to the FODs.⁷ These “FOD forces” can be used to determine the lowest energy arrangement of the FODs, in a procedure that is precisely analogous to optimizing the positions of the atoms in a molecule. The FLO-SIC method requires optimizing 3M parameters to obtain optimal local orbitals, while explicitly satisfying the LE by finding an optimal unitary transformation on the occupied orbital space involves O(M²) parameters. Formally, at least, this represents a significant computational advantage for the FLO-SIC method compared to methods that impose the LE constraint. In practice, this advantage may be enhanced by the apparent transferability of FOD positions. Results obtained in early, non-self-consistent applications of FLO-SIC suggest that the optimal FOD positions for atoms in similar bonding arrangements are similar.^8–11 Since locating optimal FOD positions is a critical step in the FLO-SIC method and because of the evidence that FODs are transferable, it is useful to determine optimized FOD positions for isolated atoms that could be used as starting points for molecular calculations. That is the subject of this paper, as we present the results of using an unbiased search method and fully self-consistent FLO-SIC calculations¹² to obtain optimal FOD positions for the atoms from Li to Kr.

In Sec. II, we describe the search procedure and other aspects of our computational methodology. We then present the optimal FOD arrangements for the atoms from Li to Kr and describe how the optimal positions evolve from atom-to-atom across the periodic table. We show that the FLO-SIC HOMO energies for the atoms are dramatically improved over corresponding DFT HOMO energies as predictors of the observed ionization energies. Finally, we also show that the optimized DFT-SIC total energies obtained with the FLO-SIC method are in excellent agreement with DFT-SIC total energies obtained previously by explicitly satisfying the LE constraints. This suggests that the FLO-SIC method yields the same solutions for the atoms as the earlier methods, while avoiding the computationally expensive constraints. We conclude the paper with a summary and an outlook regarding future FLO-SIC calculations.

In the PZ-SIC theory, the total energy of an M-electron system is written as

where $USIC[ρi]=−12∫∫ρi(r→)ρi(r→′)|r→−r→′|dr→dr→′−Uxc[ρi,0]$⁠. Here U^xc is the exchange-correlation energy functional and $ρi=|ϕi|2$ is the density associated with the ith electron orbital, $ϕi$⁠. The total electron density is $ρ=∑iρi$⁠. By varying $EDFT–SIC$ with respect to the orbitals in the standard way, the orbital-dependent, Schrödinger-like PZ-SIC equations for an M-electron system are obtained as

where the SIC Hamiltonian for each occupied orbital is

Here H⁰ is the DFT Hamiltonian and $VSIC[ρi]$ is the SIC potential for the ith orbital. $λiji=⟨ϕj|Hi|ϕi⟩$ are Lagrange multipliers introduced to insure the orthonormality of the orbitals. The specific form of H⁰ and V^SIC depends on the exchange-correlation functional used. In the calculations presented here, we use the local spin density approximation (LSDA)¹³ and an exchange-only version of that functional.

The local orbitals used for Eq. (2) are obtained in a two-step process beginning with the construction of the Fermi orbitals according to

where $a→i$ is the FOD for orbital i. It is easily shown that F_i are normalized but not mutually orthogonal.⁵ In the second step, a symmetric Löwdin orthogonalization is performed on F_i to obtain the orthonormal local orbitals $ϕi$⁠. Since the procedure of Eq. (4) uses the density matrix and the total electron density $ρ$⁠, it is unitary invariant. Any choice of orbitals spanning the same occupied space will result in the same local orbitals for a given set of $a→i$⁠.

Minimizing the self-interaction-corrected total energy, E^DFT–SIC, requires the DFT-SIC equations to be solved self-consistently¹² and the FODs to be arranged such that the FOD forces are zero. The first process optimizes the occupied orbital space and the second insures that the FLO-SIC local orbitals result in the most negative SIC correction for the given occupied space. To solve Eq. (2) self-consistently, we use the procedure described in Ref. 12. To minimize the energy with respect to FOD positions, we compute derivatives of E^DFT-SIC with respect to the components of $a→i$⁠,⁷ which can be thought of as forces on the FODs. These are then used in a gradient-based algorithm to find the local minimum of the energy with respect to the FOD positions. How these two processes are combined to obtain the minimum energy is discussed further below.

We used the following unbiased search process to find the minimum energy FOD positions. First, we created a large number (∼1000) of starting sets of FODs for each atom. For Li to Ar, the starting sets were made by placing all FODs, except for the one corresponding to the 1s orbital, at random positions inside a sphere of radius 2.0 Bohrs centered on the position of the nucleus. The 1s FOD was always placed at the center of the sphere. For K to Kr, the random starting points were created by first placing the FODs for the 1s, 2s, and 2p orbitals at the optimal positions they have in the Ne atom (see below). The remaining FODs were then placed randomly inside a spherical shell centered on 1s FOD. The inner radius of the shell was 0.7 Bohr and the outer radius was 2.7 Bohrs.

From each random starting point, we used a conjugate gradient (CG) algorithm to relax the FODs to a local minimum energy configuration. At each step in the minimization procedure, the current FOD positions are used to determine the local orbitals and compute the DFT-SIC total energy and the FOD forces. The latter are fed into the CG algorithm and the FOD positions are updated. The new positions are then used to repeat the process. The steps are repeated until the largest FOD force drops below a specified tolerance of the order of 0.001 Hartree/Bohr. Note that all FOD positions were allowed to move during the relaxation, including those that were placed at non-random positions in the starting configurations.

The FLO-SIC methodology is implemented in a development version of the NRLMOL code.^14,15 NRLMOL features extensive and highly optimized Gaussian orbital basis sets¹⁶ and an efficient variational integration scheme¹⁴ that can be adjusted to yield highly accurate numerical integrals. An extremely fine integration grid is needed to accurately compute integrals involving the orbital dependent $ViSIC[ρi]$⁠. With the combination of an extensive basis set and accurate integration, we estimate the absolute accuracy of the DFT-SIC total energies presented below to be of the order of 1 × 10⁻⁴ Ha.

To illustrate the search method, we consider the case of Cr, which has M = 24 electrons. The minimum energy spin state of the atom in the LSDA has 6 unpaired electrons or 15 electrons in the majority spin channel and 9 in the minority channel. Thus, the starting FOD configurations for Cr were created by placing the 10 non-core FODs in the majority channel and the 4 non-core FODs in the minority channel at random positions within a spherical shell around the 1s and 2s2p³ core FODs in each channel.

A typical plot of total energy vs CG step during an FOD optimization is shown in Fig. 1. The random starting arrangement of the FODs is also shown in the upper panel. (Only majority spin FODS are shown for clarity.) After approximately 900 CG updates, the energy is converged to approximately 1 mHa and the maximum FOD force is reduced to 0.005 Hartree/Bohr. The spikes seen in the graph correspond to CG steps in which the position of the 1s FOD is displaced significantly. The arrangement of the relaxed FODs is shown in the lower panel of Fig. 1. The CG optimization has clearly led to a significant change in the starting FOD positions. The optimal majority spin arrangement can be described qualitatively as a central 1s FOD, surrounded by a 2s2p³ tetrahedron, and then further bracketed by FODs forming a nearly square pattern on one side and a bent rhombus pattern on the other, plus two additional FODs that are well outside of the rest. The square and rhombus are staggered, and the vertices of the rhombus are aligned with the inner tetrahedron. The relaxed minority spin arrangement features the 1s plus 2s2p³ tetrahedron surrounded by a second 3s3p³ tetrahedron that is inverted relative to the first.

The DFT-SIC total energy is minimized in the FLO-SIC method when the PZ-SIC equations are solved self-consistently¹² (updating the occupied orbital space) and the FODs are arranged so that the FOD forces are zero (finding the optimal local orbitals spanning the occupied space). To make the searches for optimal FODs as efficient as possible, we tried two distinct approaches that differed in how the occupied orbital space is updated relative to updates of the FODs. The first is to update the occupied orbital space and the FODs sequentially by solving Eq. (2) self-consistently at each CG step, using the self-consistent electron density to compute the DFT-SIC energy and FOD forces, and feeding these into the CG algorithm to update the FOD positions. We refer to this as the SCF-CG method.

An alternative approach is to freeze the occupied orbital space while carrying out the CG optimization of the FODs, computing the FOD forces using the same total density at each step in the minimization. We tried two variants of this frozen occupied space (“frozen-CG”) approach: first by simply using the self-consistent solution of a DFT calculation to define the occupied orbital space and second by first solving the PZ-SIC equations self-consistently for the starting set of FODs and using the resulting orbitals to define the frozen occupied orbital space. [It was found in earlier calculations⁷ that the orbitals from a generalized gradient approximation (GGA) calculation are better starting points for DFT-SIC calculations that use the LSDA approximation. We continue that practice here.] Both choices of the frozen occupied space required similar numbers of CG steps to reach optimized FOD positions, so we generally adopted the simpler approach of using the GGA orbitals to define the occupied orbital space. Once the frozen-CG approach results in optimal FOD positions, we switch to the SCF-CG method to obtain a fully optimized DFT-SIC total energy. This typically requires only a small number of additional steps from the end of the frozen-CG search.

To better compare these approaches, consider how the DFT-SIC energy changes at various stages of the calculation. For the random starting FOD arrangement shown in Fig. 1, the non-self-consistent DFT-SIC total energy obtained using the GGA frozen orbital space is −1045.867 Ha. (The GGA total energy without SIC is −1042.208 Ha.) Solving the PZ-SIC equations self-consistently with this set of FODs reduces the DFT-SIC energy to −1045.902 Ha, a decrease of 0.035 Ha or approximately 0.95 eV. Optimizing the FODs using the SCF-CG method ultimately results in an energy of −1046.121 Ha, a further decrease of about 0.219 Ha or nearly 6 eV. Following the frozen-CG approach from the same random starting point, the CG optimization of the FODs lowers the initial non-self-consistent DFT-SIC energy by 0.216 Ha, to −1046.085 Ha. This is nearly the same decrease as obtained through FOD optimization in the full-SCF method. From the optimized FOD positions, a single SCF calculation lowers the energy by 0.033 Ha, nearly the same shift as obtained in the first self-consistent calculation in the full-SCF method. The additional optimization with the full-SCF method lowers the energy by only another 2 mHa, to −1046.121 Ha, the same final relaxed energy found earlier. This analysis indicates that the frozen-CG method succeeds in getting the FODs very close to the positions they have in the fully optimized (FOD positions and the occupied space) solution. Since it requires much less computer time than the SCF-CG approach, we adopted the frozen-CG method as the standard approach for obtaining the results reported below.

The arrangement of the optimal FODs for an atom is characterized by a shell structure corresponding to the principal quantum number of the electrons. Consider the arrangements for the closed-shell atoms Ne, Ar, and Kr shown in the upper panel of Fig. 2, where the FODs for orbitals of one spin are shown for each atom. For the opposite spin, the pattern is the same, although to minimize the energy, the positions are inverted through the origin of coordinates. This inversion lowers the energy by 0.001 Ha for the Ar atom, for example.

For Ne (1s²2s²2p⁶), the FODs corresponding to the 2s and 2p local orbitals are arranged on the vertices of a tetrahedron centered on the position of the nucleus. As discussed in Ref. 5, the local orbitals generated by FODs at these positions are sp³ hybrids that yield more negative SIC corrections than pure s- or p-orbitals. For Ar (1s²2s²2p⁶3s²3p⁶), the FODs corresponding to the 3s and 3p orbitals form a tetrahedron outside that of the n = 2 shell, with its vertices staggered with respect to those of the inner tetrahedron. For Kr, the FODs corresponding to the 4s and 4p orbitals (green) form yet another tetrahedron, this one well outside the positions of the n = 3 (pink) and n = 2 (light blue) FODs. The 3s, 3p, and 3d orbitals are completely filled in Kr, and the corresponding nine FODs all lie at a roughly equal distance from the nucleus in an anti-prism arrangement.

How does the optimal FOD pattern change for atoms across a row of the periodic table? To illustrate this, the patterns for the FODs in the majority spin channel for B, C, and N are shown in the middle of Fig. 2. For B, there are two n = 2 FODs in the majority channel that lie on either side of the central 1s FOD in a nearly linear arrangement that creates sp hybrid local orbitals. For C, the three n = 2 FODs lie at the vertices of a roughly equilateral triangle, centered on the 1s, creating sp² hybrids. For N, the four n = 2 FODs sit at the vertices of a tetrahedron, making sp³ hybrids. In each of these cases, the FODs in the minority channel are arranged in a similar way: a single 2s FOD outside of the central 1s FOD. Going beyond N to O, F, and Ne, the majority spin FODs retain their tetrahedral arrangement, while the minority spin FODs repeat the patterns for B–N.

A similar evolution of the valence FODs is found for the 3rd row atoms and for Ga to Kr in the 4th row. In these cases, the tetrahedral arrangement of the n = 2 core FODs is maintained.

To understand the evolution of FOD positions in the transition metal (TM) series, it is convenient to start with the arrangement for Sc ([Ar]4s²3d) shown in the bottom panel of Fig. 2. In the majority channel, the 4s FOD lies well outside the others, and the 3d FOD joins the n = 3 shell. In the minority channel, the 4s FOD also lies well outside the n = 3 and n = 2 tetrahedra. As additional 3d electron states are occupied, the corresponding FODs join the n = 3 shell. At Mn ([Ar]4s²3d⁵), the n = 3 shell is completed in the majority channel, revealing the anti-prism pattern, with the 4s FOD remaining well outside. A similar evolution repeats itself in the minority spin channel as more 3d electrons are added, culminating in Zn ([Ar]4s²3d¹⁰), which has the complete n = 3 anti-prism and a more distant 4s FOD in both spin channels.

The results shown in Fig. 2 suggest that the arrangements of the core FODs follow simple heuristic rules based on viewing the FODs qualitatively as average electronic positions. The “electrons” simply position themselves as near to the positively charged core as possible, while maintaining a maximum separation from each other to minimize the Coulomb repulsion. This implies that, for example, the four FODs of the 2s2p³ shell should form a tetrahedron. The same rule implies that a full 3s3p³ tetrahedron should be staggered with respect to the inner n = 2 tetrahedron, as seen for Ar.

By using such rules, the problem of optimizing 3M positions in a molecule may be greatly simplified. A complete 2s2p³-shell can, in principle, be described by a single FOD position and the remaining FODs can then be constructed using local tetrahedral symmetry. The chemical environment of the molecule may impose a preferred orientation on the core FOD shells; however, by rigidly rotating the shells, one could find the preferred orientations on an atom-by-atom basis.

In Fig. 3 the average distance from the nucleus to the FODs in the n = 2, 3, and 4 shells (d) is shown versus the atomic number (Z). (Only average values for the majority spin channels are shown.) Generally, the average distance to the FODs in a shell decreases with increasing atomic number. This parallels the decrease in the radial expectation value of orbitals with the same principal quantum number as the atomic number increases. It is expected that the average distances for core FOD shells will be approximately the same in molecules as in the isolated atoms. A table containing the average value of d for the shells in both spin channels is included in the supplementary material.

In Fig. 4 we show the distance d from the nucleus to the FODs for Ne-like (O²⁻, F⁻, Ne, Na⁺, and Mg²⁺) and Ar-like (Cl⁻, Ar, K⁺, and Ca²⁺) ions. These ions are the most common charge states for these atoms in molecular settings. In all cases, we created starting arrangements with the n = 2 and n = 3 FODs (for the Ar-like ions) at the vertices of tetrahedra and relaxed to optimal positions using CG. d decreases approximately linearly with an increasingly positive charge state.

It is interesting to compare the optimized DFT-SIC total energies obtained using FLO-SIC to those found previously using different approaches. Vydrov and Scuseria¹⁷ studied a number of atoms and small molecules using a direct minimization approach to solve the DFT-SIC equations using the gradient of the total energy with respect to the local orbitals. The gradient expression includes terms that stem from the orthogonality constraint on the local orbitals, effectively incorporating the effects of the localization equations so that at self-consistency, the results satisfy the LE. They presented total energies for the atoms from Li to Ar for a variety of different exchange-correlation functionals.

In Table I we compare our FLO-SIC results for the PW91 LSDA exchange-correlation functional¹³ to Vydrov’s and Scuseria’s LSDA results using the SVWN5 functional.¹⁸ [These functionals differ only slightly in their treatment of correlation, and atomic total energies computed using them are generally negligible (i.e., typically less than 1 mHa).¹⁹] Experimental total energies are also shown for comparison.²⁰ The FLO-SIC results are extremely close to the earlier results. In twelve of the sixteen cases, the total energy difference is less than 3 mHa, with the FLO-SIC energy nearly always lower in energy. For Ne and Na, the FLO-SIC result is lower by 16 and 7 mHa, respectively. These small differences may be due to the use of different basis sets. The Gaussian orbital basis sets used for the NRLMOL calculations are highly optimized¹⁶ and roughly equivalent to quadruple zeta quality. Vydrov and Scuseria¹⁷ used the 6-311 + G(3df,2p) basis set. Some of the differences might also be due to the use of different numerical integration schemes and the challenge of accurately integrating over the SIC potentials.

Pederson and Lin²¹ also carried out PZ-SIC calculations for the atoms from Li to Ar, using a method that explicitly satisfied the LE in the self-consistency scheme. They used the LSDA functional in the exchange-only (X-only) approximation. To compare to these data, we re-optimized the FOD positions for each atom using the same X-only approximation and recomputed the total energies. (Only a few CG steps were required when starting from the optimal FOD positions from the LSDA calculations.)

In Table I we present our FLO-SIC X-only results along with those of Pederson and Lin.²¹ The comparison is close, but the FLO-SIC results are systematically lower, in some cases by as much as 0.02 Ha. Again the differences are likely due to the use of better basis sets in the FLO-SIC calculations. Pederson and Lin created contracted basis sets from single Gaussian orbitals using Gaussian exponents due to Huzinaga²² and Veillard,²³ respectively. We expect the optimized NRLMOL basis sets to be significantly better. We note that the differences increase with the atomic number as would be expected if they resulted from basis set effects.¹⁶ 

As mentioned in the Introduction, HOMO levels calculated within DFT systematically underestimate the ionization energy. PZ-SIC is known to lower HOMO levels, making them a better approximation of the corresponding electron removal energies.¹ In rigorous PZ-SIC theory, the canonical orbitals are those obtained by diagonalizing the Lagrange multiplier matrix, and the corresponding eigenvalues are approximations of electron removal energies. In the FLO-SIC approach, the Lagrange multiplier matrix is not strictly symmetric (see the discussion in Sec. III D). Thus, we follow the approach of Ref. 5 and diagonalize the following, symmetrized Hamiltonian:

where σ is the spin and $Vijiσ=⟨ϕiσ|VSIC[ρi]|ϕjσ⟩$⁠. We find that the HOMO levels obtained from these eigenvalues give improved approximations of the atomic ionization potentials (IPs) compared to the corresponding DFT energies. In Fig. 5 the absolute values of the FLO-SIC HOMO energies are compared with the HOMO from the uncorrected LSDA calculations and with the experimental IPs for all the atoms from Li to Kr. The FLO-SIC HOMO energies track the IPs closely, typically lying within a few tenths of an eV. By contrast, the LSDA HOMO energies underestimate the IPs by roughly half, although the atom-to-atom trends in the LSDA HOMO levels faithfully reproduce the measured trends. It is interesting that dips in the experimental values at O, S, and Se (N = 8, 16, and 34) are better reproduced in the LSDA calculations than in FLO-SIC. It is also interesting that the SIC HOMO levels slightly overestimate the ionization energies for the s-p atoms and slightly underestimate these energies for the 3d transition metal atoms.

For a variational DFT calculation in the X-only approximation, the virial ratio, the ratio of the absolute value of the total energy to the total kinetic energy of the electrons (−E_tot/T), can be shown to equal one for an atom in the limit of an exact solution.²⁴ Because the value of the virial ratio depends in practice on factors such as the degree of self-consistency and the completeness of the basis set, it can be a useful gauge of the quality of a calculation method. For example, Pederson and Lin showed that X-only PZ-SIC calculations for atoms in which the LE are satisfied have virial ratios roughly two orders of magnitude closer to one than calculations ignoring the LE.²¹ In the FLO-SIC method, minimizing the total energy requires both a self-consistent occupied orbital space and optimized FOD positions. Thus, the virial ratio for a FLO-SIC calculation is expected to depend on how well optimized the FOD positions are.

In Table II, the values of the virial ratio are presented for self-consistent X-only FLO-SIC calculations for the atoms from Li to Ar for two cases: (1) a random arrangement of the valence FOD positions and (2) the optimized FOD positions. All calculations utilized an SCF convergence criterion of 1 × 10⁻⁶ Ha, and the same extensive, optimized basis sets were used for the PZ-SIC and LSDA calculations. As shown in Table II, the virial ratios for calculations using optimized FODs are much closer to one than for those using unoptimized FODs. For most atoms, the difference from unity is roughly two orders of magnitude smaller for the calculations using optimized FODs. This is similar to the effect seen by Pederson and Lin due to satisfying the LE.

The orbitals that variationally minimize the total energy in a DFT-SIC calculation can be shown to satisfy the localization equations (LEs).² The LE can be expressed in terms of the Lagrange multipliers as

Equation (6) implies that the matrix of Lagrange multipliers for the correct orbitals is symmetric across the diagonal. In previous PZ-SIC implementations, the LEs were satisfied explicitly as a constraint on the local orbitals. Imposing this constraint was a time-consuming step, as it involved finding the M² elements of a specific unitary transformation on the occupied space that produced the orbitals satisfying the LE. The FLO-SIC method sidesteps the LE constraints and instead optimizes the 3M FOD positions to minimize the energy. Based on this formal comparison, FLO-SIC is expected to be more efficient than previous methods, though we have not directly confirmed this.

Given that the FLO-SIC results for the total energies presented in Table I are essentially the same (within basis set-related differences) as the energies found previously in calculations that explicitly incorporate the LE constraints, it is natural to ask whether the FLO-SIC solutions also satisfy the constraints even though they are not explicitly invoked in the method. In Fig. 6 we show the values of the Lagrange multiplier matrix for the Ar atom for two cases. On the left is the matrix obtained using local orbitals corresponding to one of the random starting FOD arrangements and on the right is the matrix obtained using local orbitals generated by the optimized FOD positions. Different colors are used to indicate the relative size of the matrix elements. It is easy to see that the matrix on the left is not symmetric, demonstrating that the LEs are not satisfied for local orbitals obtained using random FODs. The matrix on the right is highly symmetric, implying that the LEs are satisfied to good approximation by the optimal local orbitals obtained in the FLO-SIC calculation. The RMS difference between the values of the corresponding off-diagonal elements is 3.95 × 10⁻² Ha for the unoptimized case and 1.4 × 10⁻³ Ha for the optimized case. Similar results are found for other atoms.

This result is intriguing. FLO-SIC was introduced as an efficient method for implementing PZ-SIC. Because it avoids the LE constraints, it is not obvious that FLO-SIC should lead to the same optimal solutions as methods that explicitly satisfy those constraints. Yet for the atoms from Li to Ar, this appears to be so.

Another interesting question regarding the FLO-SIC solutions is whether a single optimal arrangement of the FOD positions exists for each atom, or whether multiple local minima can be found, corresponding to different sets of local orbitals. Generally, the energy surface corresponding to FOD positions is complicated. For example, displacing one of the core FODs results in a much larger change in the total energy than the same displacement of a valence FOD. Furthermore, the energy surface is nearly flat for some concerted displacements of FODs. A scissor-like mode that changes the angles in the tetrahedral arrangement of the 2s2p³ FODs for Ne, without changing the distances of the FODs to the nucleus, is one example. Because of such factors, it can be difficult to fully optimize the FOD positions and relaxations from different random starting points can reach seemingly different arrangements that nonetheless correspond to the same total energy.

Figure 7 shows five arrangements for the Cr atom that illustrate the subtlety of the question of distinct minima. Each of these configurations corresponds to a DFT-SIC total energy of −2092.121 ± 0.0012 Ha, which is very close to the accuracy limit in our calculations. The main differences in the arrangements involve the placement of the two outer FODs; the orientation of the remaining inner core FODs is qualitatively similar in all configurations. Thus, for the Cr atom, the placement of two outer valence electron FODs does not change the energy in a significant way. One possible implication of the general structure of the optimal Cr FOD arrangement is that, in molecular environments, the atom will detach the outer electrons to achieve a more symmetric n = 3 shell, converting to a Cr²⁺ state.

An example of clearly distinct local minima involves the 3d and 4s occupations for transition metal atoms. In Fig. 2 the 4s FOD can be seen to lie much further from the nucleus than the n = 3 FODs. We find it possible to use this difference to obtain different minimum energy arrangements of the FODs that differ in the occupation of the 3d and 4s orbitals. For example, for Ni with a net spin of 2 electrons, one can relax to the [Ar]4s²3d⁸ configuration by placing one FOD in the minority channel well outside the n = 3 shell in the starting configuration or to the [Ar]4s3d⁹ configuration by placing all the minority spin valence FODs close to the distance of the n = 3 shell. The former configuration is lower in energy by approximately 1.7 eV. Interestingly, the LMs for both relaxed solutions are approximately symmetric, implying that both solutions satisfy the LE to a similar degree.

The FLO-SIC method⁵ is a recent implementation of the Perdew-Zunger SIC.¹ The orbitals that minimize the DFT-SIC total energy are obtained using a procedure that requires, as ingredients (i) the self-consistent total electron charge density and (ii) the correct Fermi orbital descriptors (FODs), M positions in 3-space for M electrons. In this paper, we described an unbiased method for finding optimal FOD positions for isolated atoms from Li to Kr, based on fully self-consistent FLO-SIC calculations. The optimized arrangements of the FODs display a shell-structure that corresponds to the principal quantum number of the corresponding atomic orbitals. The FOD arrangements for closed shells adopt simple patterns, such as a tetrahedron for the four 2s2p³ orbitals, that can be understood by interpreting the FODs qualitatively as average electron positions. The optimized arrangements for the atoms, especially for the core FODs, are expected to be useful as starting points for future FLO-SIC calculations for molecules, clusters, or condensed matter systems.

An interesting issue that arises in this work is whether PZ-SIC solutions obtained using the FLO-SIC methodology are the same as those obtained using methods that explicitly require the local orbitals to satisfy the localization equations (LE), a necessary condition for the variational minimization of the DFT-SIC total energy. The FLO-SIC method avoids the LE, instead of minimizing the total energy by optimizing the FOD positions. Despite this apparent short-cut, the optimal FLO-SIC total energies for Li to Ar are shown to be equal to (within basis set differences) energies obtained by methods that explicitly impose the LE constraints.^17,21 Thus, for the case of atoms, it appears that FLO-SIC results in the same solution as calculations that apply the LE constraints. Whether this is true in general is a question for further study.

A critique of the PZ-SIC formalism is that its application to functionals more complicated than the LSDA does not improve their performance.²⁵ Properties such as binding energies, for example, are worse with SIC than without it.¹⁷ Analysis of such problems has led to suggestions of schemes for scaling down SIC in situations in which semilocal functionals are expected to perform well, for example, for slowly-varying densities,²⁶ but such schemes have been found to re-introduce SIE.²⁷ The remarkable success of density functionals such as PBE,²⁸ TPSS,²⁹ and the new SCAN functional³⁰ can be attributed to the fact that they were constructed to reproduce known properties of the exact (but unknown) density functional. It is unclear how the introduction of SIC into such functionals ex post facto affects these properties. A possible way forward is to base the construction of higher level functionals on a foundation that is already self-interaction free. An efficient method for implementing SIC such as FLO-SIC could be very useful in that context.

See supplementary material for the optimized FOD positions in the FLO-SIC-LSDA calculations for all atoms from Li–Kr, a table of average FOD distances by principal quantum number, and a table of average FOD distances versus charge state for Ne-like and Ar-like ions.

The authors acknowledge helpful discussions with M. R. Pederson and T. Baruah. The computational work and initial analysis performed at CMU was supported by the U. S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award No. DE-SC0001330; a portion of the analysis was supported under award DE-SC0018331. Financial support from the Deutsche Forschungsgemeinschaft HA5070/3 is gratefully acknowledged. T.H. and J.B. were supported in part by the Office of Naval Research Grant No. N00014-16-1-2464. J.B. gratefully acknowledges support as a graduate fellow from the International Research Support Initiative Program (Pakistan).