Self-interaction correction schemes for non-collinear spin-density-functional theory

We extend some of the well established self-interaction correction (SIC) schemes of density-functional theory to the case of systems with noncollinear magnetism. Our proposed SIC schemes are tested on a set of molecules and metallic clusters in combination with the widely used local spin-density approximation. As expected from the collinear SIC, we show that the averaged-density SIC works well for improving ionization energies but fails to improve more subtle quantities like the dipole moments of polar molecules. We investigate the exchange-correlation magnetic field produced by our extension of the Perdew-Zunger SIC, showing that it is not aligned with the local total magnetization, thus producing an exchange-correlation torque.


I. INTRODUCTION
In practical (spin) density-functional theory (DFT) calculations, one needs to select an approximate functional of the density to compute the exchange-correlation energy and the corresponding potential. 1Most of the commonly employed approximations are known to suffer from the so-called self-interaction error, 2 an error that implies that the electron can interact with itself via the total electronic density.The self-interaction error can lead to problems in the prediction of the electronic properties of molecules and materials.For example, it can cause an underestimation of the bandgap of insulators and semiconductors, and an underestimation of the ionization potential and electron affinity of molecules.Thus, correcting for the self-interaction error is important for obtaining reliable DFT predictions of the electronic properties of molecules and materials. 3he search for schemes correcting the self-interaction error, known as self-interaction correction (SIC), has been pioneered by Perdew and Zunger. 2 Their proposed method, now referred to as the Perdew-Zunger selfinteraction correction (PZ-SIC), leads to an exchangecorrelation energy functional that is an explicit functional of the orbitals and, hence, an implicit density functional.Implementations of the PZ-SIC approach are often done in a generalized Kohn-Sham sense, 4 where the exchangecorrelation potential depends on the orbital on which it acts.Alternatively, and in the spirit of the original Kohn-Sham DFT, a local multiplicative exchange-correlation potential can be constructed from PZ-SIC using the optimized effective potential (OEP) technique.5 The so constructed exchange-correlation potentials have the correct asymptotic behavior and exhibit discontinuities as a function of particle number.3,6 It is possible to solve the OEP equations exactly, 7,8 but this is known to be numerically challenging, and one often resorts to the scheme introduced by Krieger, Li, and a) Electronic mail: nicolas.tancogne-dejean@mpsd.mpg.deIafrate (KLI) to approximate the full solution of the OEP equations. 9A further simplification of the KLI approach is the Slater approximation, which neglects the orbitaldependent part in the OEP equations and replaces it by an orbital-averaged term.9 The so-called globally averaged method (GAM) is defined in a similar spirit.10,11 An even more drastic approximation for the SIC consists in replacing in the PZ-SIC the orbital-dependent part directly by an averaged value for all orbitals, leading to the average-density SIC (AD-SIC).12 More recently, Perdew and coworkers proposed new schemes like the local-scaling SIC 13 which are intended to fix some of the known deficiencies of the original PZ-SIC.
To our knowledge, all of these methods have so far only be proposed and employed in the context of collinear spin DFT (SDFT).However, there exist many electronic systems in which noncollinear magnetism, spin-orbit coupling (SOC) and other relativistic effects are relevant, and often the DFT practitioners are left with no other option than to use the local spin density approximation (LSDA), which suffers from self-interaction error.5][16] This allows one to include effects stemming from the noncollinear magnetism and at the same time improve upon the LSDA.
Extending the existing SIC schemes to treat noncollinear magnetism requires care: important fundamental conditions are the local SU(2) gauge invariance of the exchange-correlation energy, and the requirement that the method properly reduces to the collinear limit.Moreover, an important question is whether the exchangecorrelation magnetic field produced by a noncollinear SIC can exert a local torque on the magnetization. 17,18If such a torque exists, it must satisfy the condition that the system cannot exert a global torque on itself (this is known as the zero-torque theorem of SDFT). 19It is the goal of this work to discuss these points.
The paper is organized as follows.In Sec.II, we present the motivation underlying our proposed SIC and extend the collinear formulation of PZ-SIC and AD-SIC to the noncollinear case.Then, in Sec.III we report numerical results obtained for several isolated systems, for which we analyze the effect of the SIC on the electronic and magnetic properties of atoms, small molecules, and clusters.We also discuss its effect on the local texture of the exchange-correction torque.Finally we draw our conclusions in Sec.IV.

II. THEORY
We begin by defining the concept of self-interaction for the general case of noncollinear spin systems.Selfinteraction is usually introduced separately for exchange and correlation.Thus, let us first consider the exact exchange energy of a system of N electrons, 20 where Tr is the trace over spin indices of the one-particle spin density matrix γ στ (r, r ′ ) = N j ψ jσ (r)ψ * jτ (r ′ ), constructed from two-component spinor Kohn-Sham orbitals, where σ =↑, ↓ and likewise for τ .Here, the double underline in γ represents a 2 × 2 matrix in spin space. 16he Hartree energy is given by where n(r) = Tr[γ(r, r)] is the total charge density of the system.
From the above definitions of E x and E H , it is straightforward to show that in the one-electron case we have where n i and m i are the single orbital charge and magnetization densities.This is the generalization of the result shown in Ref. 2 for the collinear case, and forms the basis of the self-interaction corrections that we are proposing below.
More generally, for a single orbital there is no correlation energy, so we can write that the exchange-correlation energy should fulfill the constraint Importantly, we remark here that both the exchange energy, Eq. ( 1), and the Hartree energy, Eq. ( 2), are invariant under local SU(2) rotations of the spin.We thus obtain from Eq. ( 4) that the property remains true if we rotate the orbitals such that their magnetization aligns with the z direction: where Rz m i is a symbolic operator notation for performing a rotation on the spin parts of all orbitals such that they are reckoned with respect to a given global z-axis, and then constructing the resulting orbital magnetizations.
This allows us to make the link with the collinear result, see Eq. (30) of Ref. 2. Of course, when starting from the noncollinear formulation of SDFT, one needs to break some symmetries to reduce the four-component noncollinear theory based on the variables (n, m) into a two-component collinear theory based on the variables (n, m z ).This can be achieved for instance using a uniform magnetic field of small magnitude along the z-axis, which causes the orbitals to align their magnetization along this direction.In other words, the system needs to be told to choose the z-axis as its spin quantization axis.
From this, we obtain a set of necessary conditions to be able to employ Eq. (3) to build a SIC.The first condition is that the approximate exchange-correlation functional must be locally SU(2) gauge invariant, i.e., it produces the same exchange-correlation energy independently of the orientation of the orbital magnetization.
The second condition is that the noncollinear and collinear functionals should produce the same energy for the same density, for a magnetization along the z direction.In other words, xc [n iσ , 0] appears in the definition of PZ-SIC, see below.
These conditions are naturally fulfilled by the LSDA when using the method proposed originally by Kübler et al. 21The first condition is also fulfilled by the more recently proposed noncollinear exchange meta-GGA, 20,22 which also recovers properly the result of the Becke-Roussel collinear exchange functional 23 for closed-shell systems.

A. Noncollinear Perdew-Zunger SIC
Based on Eq. ( 3), we can propose a generalization of the PZ-SIC to the noncollinear case.Let us first start by reviewing briefly the collinear case.The idea behind the PZ-SIC consists in removing all the single-electron self-interaction errors for a given density functional approximation.This leads to the energy functional In this expression, n ↑ and n ↓ refer respectively to the up and down channels of the total electronic density, and the f i,σ are occupation numbers.For each Kohn-Sham orbital φ i one needs to compute the corresponding Hartree and exchange-correlation energy from its individual density n iσ and subtract it from the energy computed from the total density.This above expression is intrinsically limited to the collinear case, but can be easily generalized to the noncollinear case.Indeed, in the latter case the exchangecorrelation functional is not a functional of the density in the two spin channels (E xc [n ↑ , n ↓ ]) but a functional of the total density n and the local magnetization m.This immediately suggests generalizing Eq. ( 6) to the noncollinear case as (7) This correction removes the self-interaction of each orbital φ i as in the collinear case.
In practice, the noncollinear PZ-SIC scheme can be challenging to implement.First of all, it requires finding the local effective potential originating from this orbital-dependent scheme, unless one wants to resort to using a generalized Kohn-Sham scheme that allows for orbital-dependent potentials. 4Finding this local multiplicative potential is usually achieved by solving the OEP equation, 5,7 or some simplified version of it like the KLI approximation. 9 more subtle complexity comes from the fact that different orbitals can produce the same density.For a typical density functional approximation like LSDA, this is not a problem.However, this becomes a well-known problem with PZ-SIC, whose results depend on the orbitals and hence vary under a unitary transformation of the orbitals, [24][25][26][27][28] unless one minimizes explicitly over all possible unitary transformations, 7,29 or use specific orbitals that make the SIC a true density functional. 30We will briefly discuss this point below with numerical examples.
Finally, let us comment on an important difference between the collinear case and the noncollinear case, which concerns the practical solution of the KLI equations to get to an approximate solution to the full OEP equation.When solving these equations, the potential is defined up to a constant, which is fixed by imposing for isolated systems that v xc,σ → 0 for r going to infinity. 9This leads to a different constant for the up and down potentials in the collinear case.However, in the noncollinear case, we end up with a single constant, as we have a 2 × 2 matrix in spin space for the potential.As a consequence, for an open-shell system without SOC, for which we can compare directly the collinear and noncollinear results, the potentials for the majority spin are very similar, but in the minority spin channel they may be different.

B. Noncollinear averaged density SIC
While the PZ-SIC is known to produce very good results, it is also known to be numerically expensive to evaluate, as one needs to solve one Poisson equation and compute the exchange-correlation energy for each occupied Kohn-Sham state, and one further needs to solve the OEP equations to obtain a local multiplicative potential needed to perform Kohn-Sham SDFT calculations.This is why several simplified methods have been proposed.Among them, the most effective method is probably the AD-SIC, which, a bit surprisingly given its simplicity, can produce excellent results for atoms compared to PZ-SIC.The motivation of this method is that if all orbitals have a similar localization, we can replace their density in Eq. ( 6) by their averaged density. 7This is particularly suited for calculations with identical atoms and pseudopotential-based simulations as orbitals are similar in these cases.However, AD-SIC suffers from a sizeconsistency problem as it is explicitly based on the number of electrons, 7 which makes it unsuitable for extended systems.In this section, we show how to generalize the AD-SIC to the noncollinear case.
In the collinear case, the AD-SIC is obtained by replacing in Eq. ( 6) the orbital and spin-resolved density n iσ by the average spin-resolved density n σ /N σ , where N σ = drn σ (r) is the number of electrons in the spin channel σ.This directly leads to the collinear AD-SIC energy functional: Following this logic, one could be tempted to average not the up and down densities of collinear SDFT, but the full spin-density matrix of non-collinear SDFT, or equivalently the local charge and magnetization densities.Inserting this into Eq.( 7), one would directly obtain However, this choice does not produce the correct collinear limit.In order to illustrate this, let us consider a Li atom in a uniform magnetic field aligned along the z direction.In this case, the system has three electrons, two residing in the 1s level, and one in the 2s level.It is straightforward to see that the one electron in the 1s (spin-channel α) and the one in the 2s level have their orbital magnetization antialigned with the external magnetic field, while the second 1s electron (spin-channel β) has its orbital magnetization aligned with the external magnetization.The densities corresponding to these states are denoted n 1s,α , n 2s,α and n 1s,β .Assuming that the approximate functional which we want corrected with AD-SIC fulfills the requirements mentioned in the introduction [SU(2) gauge invariance, and the same energy for a single orbital density with m z > 0 in the noncollinear case and for the same density in the up channel for the collinear functional] we can treat the same Li atom as a collinear electronic system with a static magnetic field along the z axis.
Let us now compute the collinear and noncollinear AD-SIC corrections for this Li atom.The AD-SIC for the collinear-spin case, Eq. ( 8), is If we use the proposed averaged density SIC as in Eq. ( 9), we get Clearly, this expression will not lead to the desired collinear limit, as seen directly from the Hartree term.However, it is possible to recover the collinear limit using the same logic as proposed originally by Kübler et al. 21for treating LSDA with noncollinear spin.By diagonalizing first the spin-density matrix, we obtain two densities, n ↑ and n ↓ , which we can average by normalizing them by their integrals (thus defining the number of "up" and "down" electrons in the frame defined by the local magnetization).Similarly to the LSDA case, the potential is computed in the local frame and independently for the up and down channels, and then rotated back to the global frame using the total magnetisation.This procedure will produce the collinear limit expected in the above Li atom example.
The direct consequence of this procedure is that both the LSDA energy/potential and the AD-SIC corrections are evaluated in the same frame, which makes this approach consistent and also invariant under local and global SU(2) rotations.However, the price to pay is that the exchange-correlation magnetic field originating from the AD-SIC correction term is aligned with the local magnetisation, meaning that no exchange-correlation torque is produced by the correction scheme.

III. NUMERICAL RESULTS
We have implemented the above equations in the realspace code Octopus 31 in order to perform tests.For the case of PZ-SIC, we only computed the solution of the OEP equations at the KLI level, using the explicit solution for noncollinear spin proposed in our recent work (see supplementary information of Ref. 20).

A. Isolated Xe atom
In order to investigate the interplay between SIC and SOC, as well as numerical and theoretical problems related to the various schemes, we first consider the case of an isolated Xe atom.We use a grid spacing of 0.30 Bohr, employing norm-conserving fully relativistic Hartwigsen-Goedecker-Hutter (HGH) pseudo-potentials. 32The simulation box is taken as a sphere of radius 12 Bohr centered at the atomic center.In Fig. 1 we show the splitting of the 5p electronic levels into 5p 1/2 and 5p 3/2 levels for LSDA, LSDA+AD-SIC, and LSDA+PZ-SIC.In all cases the collinear limit is correctly recovered for PZ-SIC and AD-SIC.We found that the inclusion of the SIC does not change how SOC affects the energy levels, and the Splitting of 5p levels of Xe due to SOC versus the spin-orbit strength computed for LDA (blue curves), LDA+AD-SIC (orange curves), and LDA+PZ-SIC (red curves).The symbols (square, circle, and triangle) indicate the results obtained for the corresponding spinunpolarized calculations.
degeneracy of the energy levels is properly described by our corrections.As visible from the figure, we nicely recover the collinear limit, indicated by the symbols in Fig. 1.We also checked that in the case of vanishing SOC strength, using a small magnetic field along x, y, or z directions produces identical results, as expected from the SU(2) invariance of our proposed formulation.However, we note that for more complicated molecules, the collinear limit is not always recovered, see below.
Let us now comment on the dependence on a unitary transformation of the orbitals used in the evaluation of Eq. ( 6) and Eq.(7).In order to reveal this, we define a new set of orbitals where U is a unitary matrix.The two sets of orbitals, {φ j } and { φi }, have the same density, but their contribution to their PZ-SIC energy is different.To illustrate this we consider here three cases: i) the "minimizing" orbitals obtained directly from the solution of the Kohn-Sham equations, ii) the result of the so-called subspace diagonalization procedure in which the unitary matrix is found by diagonalizing the Hamiltonian matrix in the subspace of the "minimizing" orbitals, iii) the localization method known as the SCDM method 33 that produces Wannier functions.In Table I we report the total energy and ionization potential of Xe for the first and the last approach for both the collinear and the noncollinear case.We find no difference between the "minimizing orbitals" and the ones obtained by subspace diagonalization.As expected, more localized orbitals produce a lower total energy and a higher ionization potential.Overall, it is apparent from these results that our noncollinear functional suffers from the same problems as the collinear formulation.One solution would be to implement a minimization of the PZ-SIC energy correction with respect to the unitary transformation U , which we defer to some future work.In the following, unless specified explicitly, orbitals from the subspace diagonalization are always employed.

B. Diatomic closed-shell systems
We continue analyzing the effect of our proposed functional on small closed-shell molecules for which SOC is known to be important for their electronic structure.It is known that SOC plays an important role on the bond length of closed-shell dimers, as well as their harmonic frequency and their dissociation energy. 34However, the choice of the functional is also important for these properties, 34 and we expect the SIC to be relevant for improving the theoretical modelling of these molecules.
We start by considering the Bi 2 molecule.We performed calculations at the experimental bond length 35 of 2.661 Å for LSDA, LSDA+AD-SIC, and LSDA+PZ-SIC.We used a grid spacing of 0.30 Bohr, employing norm-conserving fully relativistic Hartwigsen-Goedecker-Hutter (HGH) pseudo-potentials. 32The simulation box was obtained from the union of two spheres of radii 12 Bohr centered on each atoms.As shown in Fig. 2, the inclusion of the SIC does not change how SOC affects the energy levels of the molecules, and the degeneracy of the energy levels is properly described by our corrections.
As in the case of Xe, the AD-SIC properly recovers the collinear limit, while we found that the PZ-SIC becomes unstable when the SOC strength is set to zero.Indeed, in this case Bi 2 is non-magnetic, and hence any local SU(2) rotation of the spins associated with a given orbital leaves the energy unchanged but changes the potential.In order to get a converged ground state in absence of SOC, we apply a tiny magnetic field.Unlike the case of Xe, we found here two possible solutions.Aligning the external magnetic field along the molecular axis, we get the limit of vanishing SOC strength.Aligning the magnetic field perpendicular to the molecular axis, we get the same eigenvalues as in the collinear calculation.
We also performed similar simulations for other diatomic molecules using their experimental geometry, see Table II.For all molecules we employ a grid spacing of 0.3 Bohr, a radius for atom-center spheres of 12 Bohr, except for Au for which we included semi-core states and used a grid spacing of 0.25 Bohr.Overall, we find that the inclusion of the SIC drastically improves the agreement with respect to the experiment for the ionization potential, as expected from the vast literature on collinear SIC.
We also investigated the polar diatomic molecules HI, IF, PbO, and TlF at their experimental geometry and compared the dipole moments for different level of description with the experimental values.Consistent with the collinear case, 43 we found that the dipole moment on average deviates more from the experimental value when using SIC than simply using noncollinear LSDA.Importantly, the limitations of the approximation of an averaged density used to get to AD-SIC appears more clearly on the dipole moments than on the ionization energy.We also performed geometry relaxation.As found in the collinear case, 24,44 we obtain that including SIC shortens the bonds, resulting here in underestimated bond lengths compared to the LSDA, the latter being in better agreement with experimental values.

C. Magnetic cluster
We now investigate the effect of SIC on the properties of small magnetic clusters by specifically considering the iron dimer, Fe 2 . 45Clusters of this type have been widely studied by means of LSDA, including SOC (see for instance Ref. 46 and references therein).Unless stated differently, SOC is included throughout.In all calculations we employ a grid spacing of 0.15 Bohr, a radius for atomcenter spheres of 12 Bohr, and we included the semi-core states for Fe atoms.A small Fermi-Dirac smearing of 10 meV for the occupations was also used.The Fe-Fe distance was taken for the Fe dimer to be the experimental one of 2.02 Å. 47 In all cases that included both LSDA exchange and correlation energy, we found a total magnetic moment of 6µ B , in agreement with in prior works.We note that our LSDA value matches well the atomic magnetic moment reported in the pioneering work of Oda et al. 48The fact that the atomic magnetic moments computed on a sphere around the atoms decrease indicates that for Fe 2 , the SIC tends to push away the magnetization from the atomic center, while the increase of the ionization potential is consistent with an increased localization of the orbitals.This points toward a non-negligible contribution of itinerant electrons to the magnetic properties in this cluster.We also computed the values for exchangeonly LSDA, together with SIC corrections.The total magnetic moments are not properly predicted in these cases, demonstrating the key importance of correlations in order to obtain reliable magnetic structures.We finally turn our attention to the exchangecorrelation torque τ (r), defined as where m is the local magnetization density and B xc is the exchange-correlation magnetic field.We computed this quantity using LSDA and LSDAx with PZ-SIC, and also with the Slater potential.As a reference here, we consider the Slater potential, which was shown to give reasonable results compared to the result of exact-exchange potential computed at the level of KLI. 20From our results (see Figs. 3a and d) the Slater potential produces a small exchange-correlation torque around the atoms, where the symmetries of the system are clearly apparent.Our results for PZ-SIC (Figs. 3b, c, e, and f) show that PZ-SIC also produces a non-vanishing torque around the atoms.
While it shows, as required by the zero-torque theorem, alternating positive and negative patterns that are also in accordance with the symmetries of the system, the overall shape and magnitude strongly differs from what is obtained from Slater potential.Importantly, we want to stress here that like the energy, the torque obtained from PZ-SIC depends upon the unitary transformation of the orbitals.This quantity therefore needs to be analyzed with great care, and we aim in the future at implementing a minimization over unitary transformations in order to eliminate this ambiguity, similarly to prior efforts. 29

IV. CONCLUSIONS
To summarize, we presented how to extend some of the existing SIC approaches to the case of non-collinear spins.We then analyzed numerically how these non-collinear SIC schemes behave for various closed-shell and magnetic systems.Overall, we found that our noncollinear schemes exhibit similar advantages and deficiencies as the collinear ones.The ionization energies are improved, but bond lengths are found to be worse than those obtained for LSDA.When the localization of individual orbitals is important, the AD-SIC performs poorly for observables that depend on local orbitals, like dipole moments or magnetic moments.
We further demonstrated that PZ-SIC for noncollinear spin can produce a non-negligible exchange-correlation torque around the magnetic atoms, but we found large differences in the magnitude and texture of the exchangecorrelation torque compared to the result of the Slater potential.
Overall, our work opens the door to a better description of the electronic and magnetic properties of systems when noncollinear effects are important, but we note that some further work, including the computation of accurate benchmarks, is needed in order to get reliable results for the collinear and noncollinear PZ-SIC schemes.Once such SIC schemes are fully established we expect them to become a useful tool for the description of materials with noncollinear magnetism.Square modulus of the four highest occupied πu orbitals of Bi2 obtained for a magnetic field aligned with the molecular axis a)-d); or perpendicular to it e)-h).The plane shown in the figures is the plane perpendicular to the molecular axis.

FIG. 3 .
FIG. 3. Exchange torque for Fe2.Top panels: The y component of the local exchange torque τ (r) in the plane y = 0, computed from a) Slater, b) LSDA with PZ-SIC, c) LSDAx with PZ-SIC.Bottom panels: the same as top panels, but showing x component of the torque in the z = 0 plane.

TABLE I .
Total energy Etot and ionization potential Ip, in Hartree, for the collinear and collinear cases using different orbitals, as explained in the main text.

TABLE II .
Ionization potentials, in eV, of diatomic systems using their experimental geometry, including SOC, for different energy functionals.

TABLE III
a Ref.35

TABLE IV .
Electronic and magnetic properties of Fe2 for different energy functionals.Ionization potential (Ip) is given in eV, and the total (M ) and local magnetic moments (|m|) are given in µB and are obtained by integrating the density on a sphere of radius 1.909 Bohr around the atoms.Exchangeonly LSDA (LSDAx) results and also reported.