Noniterative Fermi–Löwdin Orbitals for Self-Interaction Correction

TL;DR

This paper introduces a faster method for correcting self-interaction errors in quantum chemistry calculations, improving efficiency without losing accuracy.

Contribution

The novel noniterative Fermi–Löwdin orbital self-interaction correction (NIFLOSIC) eliminates iterative relaxation of FODs.

Findings

NIFLOSIC generates localized orbitals slightly more compact than grid-based methods.

NIFLOSIC improves frontier molecular orbital energies and dipole moments significantly.

NIFLOSIC reduces computational cost while reproducing results from fully self-consistent FLOSIC.

Abstract

We introduce the
noniterative Fermi–Löwdin orbital
self-interaction correction (NIFLOSIC) method as a computationally
efficient alternative to traditional Fermi–Löwdin orbital
self-interaction correction (FLOSIC) by eliminating the need for iterative
relaxation of Fermi orbital descriptors (FODs). This is accomplished
using the selected columns of the density matrix localization scheme
[J. Chem. Theory Comput.
2023,
19, 8572] and by exploiting the relationship between the
electron localization function and FODs [J. Chem. Phys.
2025,
162, 144105]. The approach produces
localized orbitals that are slightly more compact than grid-based
selected columns of the density matrix orbitals and generates FODs
in a single, noniterative self-starting step, following density functional
theory calculations. Within a generalized Kohn–Sham framework,
full relaxation of the density minimizes the Perdew–Zunger
energy functional, yielding self-interaction corrected densities and
orbitals. NIFLOSIC reproduces results from fully self-consistent FLOSIC
calculations, while significantly reducing computational cost. Although
the total electronic energy is not suitable for thermochemistry, benchmark
tests across diverse molecular systems demonstrate that NIFLOSIC significantly
improves frontier molecular orbital energies and dipole moments, establishing
a practical and scalable approach for large-scale electronic structure
applications where self-interaction correction is needed.

Figures

Total localized orbital variance, σ, for a set of 676 molecules. The abscissa corresponds to Foster–Boys (FB) orbitals (minimum possible variance) while the ordinate shows the noniterative Fermi–Löwdin orbitals (NIFLO) and selected columns of the density matrix, grid version (SCDM-g) orbitals. The inset shows a zoom in the central region of the plot. The red and blue lines show the best linear fit for each class of localized orbitals. The 1:1 line is shown in dashes.

HOMO (left panel) and LUMO (right panel) energies in eV for each method. The set of molecules is composed of CH4, H2O, SiH4, CO, CO2, SO2, N2, P2, NaCl, C2H6, CF4, H2NNH2, HOOH, H2CO, C2H2, and C5H5N5O. Dashed lines show the 1:1 agreement.

Tables

1. HOMO and LUMO Mean Absolute Deviation (MAD) and Mean Deviation (MD) with Respect to FLOSIC Values for the Data Shown in Figure

	HOMO		LUMO
molecule	DFT	NIFLOSIC	DFT	NIFLOSIC
MAD	5.85	0.32	0.24	0.24
MD	5.85	0.15	–0.18	0.12

2. Dipole Moments (in Debye) Calculated Using DFT, NIFLOSIC, and FLOSIC Methods

molecule	DFT	NIFLOSIC	FLOSIC
H2O	1.80	1.87	1.93
CO	0.29	0.89	0.69
SO2	1.48	1.89	1.79
NaCl	8.53	8.99	9.00
H2NNH2	2.63	2.79	2.83
HOOH	1.52	1.57	1.61
H2CO	2.23	2.96	2.61
C5H5N5O	6.55	6.88	6.73

3. Summary of FLOSIC and NIFLOSIC Calculation Details

	NIFLOSIC	FLOSIC
t/t DFT	13	1326
|g| (mE h /Bohr)	3.84	0.07
ΔE (mE h )	-	–26.7
ΔR (Bohr)	-	0.17

Full text

Introduction

Density functional theory (DFT) ?,? has become a cornerstone in electronic structure calculations, largely due to its effective balance between computational cost and accuracy. ?,? However, a significant limitation of approximate exchange–correlation (XC) functionals is their failure to eliminate the self-interaction (SI) of electrons, leading to the well-known self-interaction error (SIE). This error underlies several recognized shortcomings of lower-rung density-functional approximations, such as producing incorrect dissociation curves for asymmetric molecules, ?,? or predicting overly high energies for the highest occupied molecular orbitals (HOMO) in negatively charged species, which results in an artificially positive charge.?

In the early 1970s, Lindgren? proposed correcting HOMO levels by modifying the local density approximation (LDA) Hamiltonian with a self-Coulomb interaction. In 1981, Perdew and Zunger introduced a method to eliminate the one-electron SIE on an orbital-by-orbital basis.? Commonly referred to as the Perdew–Zunger (PZ) self-interaction correction (SIC), this approach modifies the energy functional to depend explicitly on individual orbitals

$$[eqn]$$

In eq, n _ i _ ^σ^ denotes the spin density of a single orbital (σ = ↑, ↓), while E XC and E H are the exchange–correlation and Hartree energies, respectively. Minimizing this functional results in localized orbitals. Despite its conceptual appeal, the PZ-SIC method has not seen widespread use in routine calculations, in part due to the computational expense of minimizing the energy, which requires variations of the individual orbital densities in E PZ‑SIC. Moreover, applying PZ-SIC alongside conventional approximate functionals often degrades the accuracy of various properties. This challenge has motivated the development of DFT approximations compatible with SIC. ?−? ? ? ? ? ?

Several strategies have been developed to minimize E PZ‑SIC, many of which involve fully relaxing the (localized) orbitals using diverse techniques. These approaches differ primarily in how they handle the additional variational parameters required to describe localized orbitals. Some methods enforce localization conditions directly through localization equations, ?−? ? ? ? while others adopt alternative strategies such as the constrained gradient search proposed by Vydrov and Scuseria,? the two-step unitary optimization method introduced by Lehtola and Jónsson, ?,? or the Koopmans-compliant functional framework developed by Ferretti et al.? Additionally, several implementations utilize variants of the optimized effective potential (OEP) method. ?−? ?

An alternative method for minimizing E PZ‑SIC involves the use of Fermi–Löwdin orbitals (FLOs). ?,? This technique, known as FLOSIC, is based on the construction of localized orbitals by parametrizing them as Fermi orbitals?

$$[eqn]$$

where a denotes spatial points referred to as Fermi orbital descriptors (FODs), and n(r′, r) is the density matrix (with spin indices omitted for clarity). Originally proposed by Luken et al. ?,? to build nonorthogonal localized orbitals, this approach takes advantage of the real-space decay of the density matrix in systems with a finite energy gap, ?,? and identifies the rows of a scaled density matrix as localized orbitals. The nonorthogonal Fermi orbitals f _ a _(r) are orthogonalized using the Löwdin scheme,? resulting in Fermi–Löwdin orthonormal orbitals. The denominator in $[eqn]$ in eq ensures a balanced weighting during symmetric orthogonalization by normalizing the Fermi orbital. Since FLOSIC introduces FODs as additional variational parameters, minimizing E PZ‑SIC requires full relaxation of both the electron density and all FODs. Early implementations of FLOSIC employed Jacobi (or Givens?)-type rotations to eliminate the overlap between occupied and virtual orbitals at each self-consistent iteration.? More recently, a mean-field approach for density relaxation within FLOSIC has been proposed.? Other implementations based on unified Hamiltonian frameworks and effective potentials have also emerged. ?−? ? ? ? In these methods, FODs are relaxed between self-consistent field cycles in a double-loop fashion, keeping the density fixed during each FOD update. Alternatively, one can reformulate the problem as a FOD minimization where each evaluation of the energy and its gradients is performed after self-consistency. Importantly, it has been shown that given a set of FODs, and due to the localized nature of the FLOs, the cost of a FLOSIC calculation can scale similarly to that of DFT in the large system regime.? Thus, methods that allow the derivation of FODs directly from inexpensive evaluations of the Kohn–Sham (KS) orbitals, or alternatively density matrices, would be highly desirable.

A class of orbital localization methods that explicitly avoid demanding iterative minimization schemes has been proposed. ?−? ? In particular, the selected columns of the density matrix grid-based (SCDM-g) method provides a noniterative approach to orbital localization by selecting columns of the density matrix evaluated on a real-space grid,? much like the original idea of Luken and Beratan? The SCDM-g localization technique, which constructs nonorthogonal localized molecular orbitals using a QR decomposition, is independent of predefined atomic orbitals, and its performance is notably robust with respect to grid resolution and basis set contraction schemes. The nonorthogonal localized orbitals can be seen as the density matrix columns that minimize the overlap matrix condition number. These orbitals are then symmetrically orthogonalized to give the set of orthonormal localized molecular orbitals (LMOs).

Using the locality of the density matrix in a local basis set representation, we have introduced a method called DOCSIC (for density matrix as orbital coefficients SIC),? as a simplification of the FLOSIC method. DOCSIC utilizes columns of the orbital representation of the density matrix to construct localized orbitals in PZ-SIC calculations, avoiding the need to incorporate additional variational parameters. This approach enabled a self-consistent generalized KS implementation of SIC and offers a fast alternative to traditional self-interaction correction schemes for cases where SIE is dominant, such as those in the SIE 4 × 4 test set of Grimme et al.? for binding energies. However, DOCSIC provides only moderately localized orbitals and corrects HOMO and LUMO energies only slightly with respect to DFT.

In this work we introduce a new approach to FLOSIC calculations that utilizes FODs and constructs FLOs based on a variation of the SCDM-g method. The method can be used to obtain localized orbitals, determine FODs for FLOSIC calculations, and allows for relaxing the density matrix to minimize the PZ energy functional.

Theory and Implementation

The FLOSIC approach generates localized orbitals by utilizing Fermi orbitals, as defined in eq. These are subsequently symmetrically orthonormalized, resulting in FLOs. Within the PZ-SIC framework, these orbitals are used to remove SIE on an orbital-by-orbital basis. As mentioned in the Introduction, a key ingredient in the FLOSIC approach is the determination of the set of FODs needed to construct these Fermi orbitals. In practice, for a given molecular system, an initial set of FODs is determined (either empirically or from existing data) and then iteratively relaxed to minimize the PZ energy functional. In a recent work, some of us have shown a formal link between FODs and the electron localization function (ELF).? Specifically, ELF values equal to one necessarily imply stationary FODs, while ELF critical points empirically resemble the FOD structure obtained from FLOSIC calculations.

Through a rank-revealing QR factorization, the original SCDM-g method constructs the most linearly independent orbitals, favoring points located where the electron density is significant (or equivalently, where the MOs have a significant value).? Using the connection between the ELF and FODs, here we introduce a methodology based on a modified SCDM-g localization scheme ?,? to determine a set of FODs for FLOSIC calculations in a self-starting, noniterative fashion. To this end, we normalize the MO values and multiply by the ELF at each grid point, aiming to favor the selection of points where the ELF is maximal, which is a desirable feature for the FLOSIC method as discussed before.? This scheme balances the ability of the rank-revealing QR factorization to find most linearly independent Fermi orbitals and maximal ELF points at the same time. The process is as follows:

• 1.From a DFT calculation, evaluate all the N molecular orbitals (MO) {ψ_ b _} (N occ occupied and N virt virtual), the density n(r), and the ELF, F ELF(r) on all the K points of a quadrature grid {r _ j _} (the standard integration or a coarser grid) and evaluate the N × K matrix $[eqn]$ . Normalize this quantity to obtain $[eqn]$ . Evaluate $[eqn]$ on the grid.
• 2.Perform a rank-revealing QR factorization with pivoting on the rectangular N × K matrix A of elements A _ bj _. This QR factorization with pivoting provides a solution to AΠ = QR, where Q and R are orthonormal and upper-triangular matrices, respectively, and Π is a permutation matrix. Π is chosen so that the diagonal elements of R are nonincreasing, i.e. |R _ ii | ≥ |R _ i+1 i+1|. Keep only the first N occ elements of the permutation matrix, Π occ.
• 3.The grid points {a | r _ j _ ∈ Π occ} are the FODs to be used to construct Fermi orbitals,

$$[eqn]$$

The FLO coefficient matrix X in the atomic orbital (AO) basis {ϕ_μ_} can be evaluated as X = GO ^1/2^, where G = PY, O = Y ^†^ PY is the Fermi orbitals overlap matrix, and the matrix Y has elements $[eqn]$ . The single-orbital density matrices P _ a _ can be evaluated from the elements of X as $[eqn]$ .

The resulting noniterative Fermi–Löwdin orbital (NIFLO) single-orbital coefficient matrix X and density matrices P _ a _ can be used to evaluate the orbitals and their corresponding densities, while the grid points a can be used as FODs for FLOSIC calculations. The effective generalized KS Hamiltonian can be obtained as

$$[eqn]$$

where the first term in the r.h.s of eq is the standard KS Hamiltonian, and E SIC refers to the second term on the r.h.s of eq. This derivative can be written as

$$[eqn]$$

We have implemented the NIFLO approach using an in-house modification of the PySCF electronic structure code.? The derivatives in eq were evaluated using the automatic differentiation Autograd from Torch.? We note that, alternatively, one can obtain H _ SIC _ using the procedure outlined in ref ?. The current implementation assumes real-valued density matrices and hence leads to real-valued localized orbitals. With the generalized multiplicative Hamiltonian matrix in eq, the self-consistent procedure can be performed as in any standard DFT calculation. It should be emphasized that this procedure employs the grid points {a} utilized to generate the localization transformation from orbitals obtained from a DFT (not SIC) calculation. However, since this transformation depends also on the density matrix, the transformation itself is updated at every iteration. We refer to this approach as NIFLO self-interaction correction, or NIFLOSIC. A version of the code, including the initial FOD generator, NIFLO localized orbitals, and NIFLOSIC and standard FLOSIC calculations, is publicly available in https://github.com/peraltajuan/NIFLO.

Results

All calculations in this work utilize the density functional approximation of Perdew, Burke, and Ernzerhof (PBE) ?,? with the def2-TZVPD basis set. ?,? We first assess the locality of the resulting NIFLOs by comparing the total orbital variance $[eqn]$ with the SCDM-g and Foster–Boys (FB) orbitals. The latter minimize the mean extension of the orbitals from their centroids, and thus have the minimum possible orbital variance. In Figure, we show this comparison for a set 676 randomly selected molecules from the GDB-11 database containing from 3 to 11 first-row elements.? The plot shows that the NIFLO method yields localized orbitals that are more compact than their SCDM-g counterparts. Both localized orbitals are significantly more localized than noniterative Cholesky orbitals (not shown in Figure).?

Total localized orbital variance, σ, for a set of 676 molecules. The abscissa corresponds to Foster–Boys (FB) orbitals (minimum possible variance) while the ordinate shows the noniterative Fermi–Löwdin orbitals (NIFLO) and selected columns of the density matrix, grid version (SCDM-g) orbitals. The inset shows a zoom in the central region of the plot. The red and blue lines show the best linear fit for each class of localized orbitals. The 1:1 line is shown in dashes.

SIE leads to an artificial electron delocalization in functionals of LDA and generalized gradient approximation (GGA) families, which systematically raises the energy of the HOMO and lowers the energy of the lowest unoccupied molecular orbital (LUMO). This results in a well-known unphysical underestimation of the HOMO–LUMO gap in molecules and band gap in solids, often associated with the delocalization error in DFT. ?−? ? ? ? Thus, it makes sense to assess the ability of the NIFLOSIC method to reproduce frontier orbital energies as calculated by FLOSIC. To this end, we compare the HOMO and LUMO energies for a diverse set of molecules using standard DFT, the NIFLOSIC, and fully self-consistent FLOSIC approaches. Figure summarizes the results for 17 representative molecules (structures taken from ref ?), ranging from simple hydrides to more complex systems. The HOMO energies predicted by conventional DFT are significantly underestimated with respect to FLOSIC, as it is well-known, with a mean absolute deviation (MAD) of 5.90 eV with respect to the FLOSIC energies. NIFLOSIC HOMO energies, in contrast, align much better with FLOSIC HOMO energies, with a MAD of 0.30 eV. A statistical summary of these results is shown in Table. The effect of removing SIE in LUMO energies is less pronounced, as shown in Figure and this is reflected on both the NIFLOSIC and FLOSIC methods.

HOMO (left panel) and LUMO (right panel) energies in eV for each method. The set of molecules is composed of CH4, H2O, SiH4, CO, CO2, SO2, N2, P2, NaCl, C2H6, CF4, H2NNH2, HOOH, H2CO, C2H2, and C5H5N5O. Dashed lines show the 1:1 agreement.

1: HOMO and LUMO Mean Absolute Deviation (MAD) and Mean Deviation (MD) with Respect to FLOSIC Values for the Data Shown in Figure

SIE causes spurious charge transfer within a molecule and results in inaccurate predictions of the molecular dipole moment. Correcting for SIE generally increases the degree of charge transfer, thus increasing the size of calculated dipole moments.? In Table we show the molecular dipole moments calculated for the 8 nonzero dipole cases in Table. NIFLOSIC corrects the DFT dipole moment toward the FLOSIC value in all cases, indicating that the resulting density from the self-consistent NIFLOSIC calculation effectively captures the self-interaction correction.

2: Dipole Moments (in Debye) Calculated Using DFT, NIFLOSIC, and FLOSIC Methods

The results in Table showcase the potential of the NIFLOSIC method for performing self-interaction corrected band structure calculations in solids at a fraction of the cost of full FLOSIC calculations. As a reference, for the set of molecules analyzed in this work, NIFLOSIC calculations take on average about 13 times more than DFT calculations, while the full FLOSIC calculations take about 1300 times that (as shown in Table). This large difference is mainly due to the additional time required to fully relax the FODs in FLOSIC calculations, while in NIFLOSIC FODs are determined noniteratively only once.

3: Summary of FLOSIC and NIFLOSIC Calculation Details

In this work, the grid points {a} determined for the construction of the NIFLOs are subsequently employed as initial positions for the FODs, while the converged NIFLOSIC density serves as the initial guess for the FLOSIC calculations, underscoring a practical applicability of the method. This allows us also to inspect the behavior of key quantities upon relaxing FODs. Table shows the average energy decrease for the set of 17 molecules in Table. One immediate observation from this is that the NIFLOSIC total energies are close to the fully relaxed FLOSIC energies, but not sufficiently close to perform general thermochemistry (the gain in energy is on average ∼17 kcal/mol or 0.73 eV), except in cases where SIE is dominant? or potentially with SIC schemes where the total energy is relatively insensitive to the FOD positions. The average displacement per FOD in the self-consistent process is ∼0.17 Bohr and the FOD gradient norm for each FOD decreases in average from 3.84 to 0.07 mE_ h /Bohr. The NIFLOSIC method provides a fast nonempirical self-starting set of orbitals and FODs for FLOSIC calculations, at a cost that is on average about 12 times the cost of a DFT calculation for the small molecules considered here. In contrast, the full FLOSIC optimization computational cost using a tight tolerance of 10^–5^ mE h _/Bohr for the projected FOD gradient is about 100 times the cost of a NIFLOSIC calculation.

Summary and Outlook

The NIFLOSIC method provides a computationally efficient alternative to traditional FLOSIC by eliminating the need for iterative relaxation of Fermi orbital descriptors. The method is shown to produce localized orbitals that are slightly more localized than SCDM-g orbitals, and also to generate FODs in a single, noniterative self-starting step after DFT calculations. Fully relaxing the density to minimize the PZ energy functional in a generalized KS scheme provides self-interaction corrected density and orbitals. NIFLOSIC yields results that closely reproduce fully self-consistent FLOSIC calculations while remarkably reducing the computational cost. The method can also be utilized to produce initial FODs and densities for FLOSIC calculations. One weakness of the method is that the total electronic energy may be unsuitable for thermochemistry evaluation. However, benchmark tests across diverse molecular systems demonstrate that NIFLOSIC significantly improves HOMO energies and dipole moments, positioning it as a practical and scalable approach for large-scale electronic structure applications where self-interaction correction is essential.

This work offers a method to derive FODs and FLOs directly from KS orbitals, eliminating the FOD optimization and potentially restoring DFT-like computational scaling. Although our results discourage the use of this method with current density functionals for thermochemistry purposes, looking into the future, the resulting computational efficiency encourages the derivation of new explicitly self-interaction-corrected density functionals that are insensitive to the localization scheme, rather than appending new terms onto density-functionals that were developed to be used without self-interaction corrections.