PyFLOSIC: Python-based Fermi–Löwdin orbital self-interaction correction

The continuously growing availability of open-source software has given rise to an ongoing paradigm shift in quantum chemical software development. At variance to the traditional model of monolithic programs, where new algorithms need to be re-implemented separately from scratch in each code, present-day computational science is proceeding in leaps and bounds, thanks to the thriving ecosystems of small projects dedicated to solving well-defined subproblems^1–6 that can be easily combined to form a code that is greater than the sum of its parts (see Refs. 7 and 8 and references therein for further discussion). The efficiency of the new development models comes from not having to “re-invent the wheel” within every program; instead, standardized, modular tools designed to solve specific problems can be reused across the board. This kind of cleaner code design and efficient code reuse has enabled fast development and streamlined the adoption of new methodologies, since features added to any one component become straightforwardly available in everything built on top of them.

This work describes pyflosic, an open-source, python-based implementation of the Fermi–Löwdin orbital self-interaction correction (FLO-SIC).^9–12 The core routines of Pyflosic were developed during Lenz Fiedler’s master’s thesis,¹³ and the code is freely available on github (https://github.com/pyflosic/pyflosic). pyflosic builds on the python simulation of chemistry framework (pyscf),^7,14 which is an open-source electronic structure and quantum chemistry code written primarily in python. As a general-purpose quantum chemistry program, pyscf already contains a vast number of methods.^7,14 For example, self-consistent field (SCF) approaches such as Hartree–Fock and density functional theory^15,16 (DFT) are available with various choices for spin treatment, density fitting techniques, and numerical quadrature methods. Second-order orbital optimization is also available in pyscf.¹⁷ 

Since pyscf is highly modular and its parts can be imported like any other python module, it is rather easy to add new functionalities to it or to combine it with existing software,¹⁴ as has already been demonstrated by various authors.^2,5,6 This is also the strategy that was adopted for pyflosic. Before pyflosic, the FLO-SIC method has only been available in the reference implementation within the Naval Research Laboratory Molecular Orbital Library (nrlmol).^18–26 There is also a FLO-SIC implementation in a developer branch of the nwchem program,^27,28 but to the best of our knowledge, this implementation is not currently publicly available.

Similar to nrlmol and nwchem, pyscf—and therefore pyflosic—uses a Gaussian-type orbital (GTO) basis set.^7,14 However, in contrast to the nrlmol program, nwchem and pyscf (and pyflosic) are able to routinely handle basis functions with high angular momentum. The most-commonly used all-electron and effective core potential GTO basis sets are already available within PYSCF, and additional ones can be downloaded from the Basis Set Exchange²⁹ and parsed in from several formats.

The importance of the tractability of calculations within large basis sets is underlined by the recent fully numerical benchmark studies that have shown that the reproduction of atomization energies even within DFT may require several shells of polarization functions.^30–33 pyscf contains fast routines for the evaluation of molecular integrals via the libcint library,¹ which are also used in pyflosic. Moreover, pyflosic inherits OpenMP parallelization from pyscf, enabling the efficient use of multicore computation architectures.⁷ By enabling the use of large basis sets in FLO-SIC calculations, thanks to the fast elementary routines in pyscf, pyflosic enables reliable computational studies even of molecules for which very large and flexible basis sets are required, such as SO₂ and SF₆,^30,31,33 as will also be demonstrated later in this work. Recently developed powerful approaches to handle significant linear dependencies in the underlying molecular basis set^34,35 are also available in pyflosic through pyscf, enabling accurate FLO-SIC calculations even in pathologically over-complete basis sets.^33–35 

Equally importantly, again at variance to the previous implementations of FLO-SIC that feature in-house implementations of only a handful of density functionals, pyflosic contains (again via pyscf) interfaces to the libxc³ and xcfun³⁶ libraries of exchange-correlation functionals, which grant access to a wide range of hundreds of local density approximations (LDAs), generalized-gradient approximations (GGAs), meta-GGAs, hybrid functionals, and range-separated hybrids.⁷ At the moment, all LDAs, GGAs, and meta-GGAs provided through pyscf can be used in pyflosic; for a thorough list of the functionals implemented in libxc and xcfun, we refer to the libraries’ respective documentations.³⁷ Custom linear combinations of exchange-correlation functionals can also be used within pyscf, further expanding the capabilities of pyflosic.

As an add-on, pyflosic inherits all of the important features of pyscf. Because Pyflosic is implemented as a collection of python modules alike pyscf,¹⁴ those already familiar with pyscf do not have to learn a new package-specific input format to run FLO-SIC calculations with pyflosic. Users of pyflosic can also rely on the full power of python—one of today’s most commonly used and taught programming languages that has a massive and vibrant community—to fulfill their every need. The availability of all python language tools within the input script allows for elaborate work schemes, e.g., the combination of calculation, evaluation, and plotting routines.¹⁴ Calculations can also be done interactively in the python interpreter shell.¹⁴ 

Having introduced and motivated pyflosic, we will next discuss the theory of FLO-SIC in detail in Sec. II. We will start out by motivating the use of self-interaction corrections in general (Subsection II A), continue by presenting the approach of Perdew and Zunger in specific (Subsection II B), introduce Fermi–Löwdin orbitals to undertake the self-interaction correction following the prescription of Perdew and Zunger (Subsection II C), and discuss the mandatory initialization and optimization of the Fermi-orbital descriptors (FOD) that parametrize the Fermi–Löwdin orbitals (Subsection II D). The schemes used in pyflosic for optimizing the FLO-SIC density are discussed in Subsection II E. The code is showcased in Sec. III, starting with example inputs for running a calculation on tetracyanoethylene. The structure of the tetracyanoethylene example follows that of the calculation: first, the FODs are initialized (Subsection III A), and then, the electron density and the FODs are optimized (Subsection III B). A discussion on repeated calculations follows in Subsection III C. As a practical application of the code, in Subsection III D, we perform an in-depth basis set convergence study of the atomization energies of SO₂ and SF₆ that have been found to be challenging cases in the literature, as was already discussed above. The article concludes with a summary and outlook in Sec. IV. Atomic units are used throughout the paper, if not stated otherwise.

We will adopt the notation introduced by Lehtola and Jónsson³⁸ in the following. However, in contrast to Ref. 38, vectors and matrices are distinguished in the presently used notation: all vectors are expressed in bold, non-italic letters (a), and all matrices are given in bold, italic letters (R). All equations are written for real-valued orbitals.

DFT has become one of the standard methods in computational materials science, condensed matter physics, and chemistry, thanks to its combination of reasonable accuracy with computational efficiency.^39,40 However, currently available density functional approximations (DFAs) are well-known to fail in a number of situations.^40–42 Although the description of systems exhibiting significant static correlation remains an open problem, posing limitations on studies of many transition metal complexes and molecules with stretched bonds, challenges exist also in the absence of static correlation. For instance, localized states and negatively charged species such as F⁻ are incorrectly described by pure (i.e., local or semi-local) density functionals.

These shortcomings come from the fact that DFAs include spurious interactions of electrons with themselves, known as the self-interaction error, which causes the electron density to delocalize. Delocalization error has long been identified as a key issue in the application of DFT onto the study of chemical systems; for instance, the barrier height of the simplest hydrogen abstraction reaction H + H₂ ↔ H₂ + H already poses a challenge, as many functionals overestimate the stability of the intermediate H₃ state.⁴³ The error can often be significantly reduced by including a fraction of exact exchange as in (possibly range-separated) hybrid functionals; however, this does not fully remove the problems with self-interaction.

An exchange-correlation functional that fulfills three constraints may be free of self-interaction for many-body systems. It needs to (i) be free of self-interactions for all one-electron densities, (ii) provide the correct piecewise linearity (PWL),⁴⁴ and (iii) recover the correct asymptotic limit of the potential. As no DFAs that fulfill all of these constraints completely are available, self-interaction corrections (SICs) aim to rectify some of the aforementioned issues for the available DFAs. Since the early formulation of SIC,⁴⁵ SICs have been shown to be able to significantly reduce the errors of the underlying exchange-correlation functional for the first two of the aforementioned issues, namely, the delocalization error (i) and the lack of piecewise linearity (ii).^{11,38,46–55} The lack of piecewise linearity is the reason the highest occupied molecular orbital (HOMO) energies from Kohn–Sham (KS) DFAs yield poor estimates for the ionization potential as IP = −ε_HOMO; in contrast, the so-called Koopmans-compliant (KC) functionals^50,56 can deliver especially accurate results for IPs.

Among the multitude of current implementations of SIC, most employ real-valued orbitals (RSIC).^57–64 However, it has been shown that SIC based on complex-valued orbitals (CSIC) has several advantages;^{48,53,65–67} in fact, Lehtola, Head-Gordon, and Jónsson⁶⁸ showed some time ago that the use of complex-valued orbitals is actually mandatory to properly minimize the SIC functional (implementations of RSIC and CSIC that support arbitrary basis sets and exchange-correlation functionals similarly to pyflosic are freely available in the erkale program^69,70).

Another variant of SIC is the FLO-SIC approach mentioned in Sec. I, which has been the subject of various studies during recent years^{9–12,27,28,55,71–86} and is also the topic of the present work; it will be discussed in depth below in Subsections II C–II E. It is first worth mentioning, however, that despite many successes of SIC (some of which were referenced above), there are also several results that show SIC degrading the performance of higher-rung exchange-correlation functionals.^50,62,87 This is the so-called “paradox of self-interaction correction” discussed in the comprehensive summary of Perdew et al.,⁵² which remains a challenge for the future. For instance, it has been pointed out that the Perdew–Zunger approach⁴⁵ does not completely eliminate one-electron self-interaction;⁸⁸ however, the path for rectifying the remaining error is still unclear. Despite these theoretical paradoxes, self-interaction corrections have been shown to be useful in many cases, which is why we will not consider this problem further in this work and will instead carry on with the Perdew–Zunger approach.

The total energy in the KS-DFT formalism is given by

where E_KS[n^α, n^β] is the total KS energy, T_s[n^α, n^β] is the kinetic energy of the non-interacting system, V[n] is the external potential energy, J[n] is the Coulomb energy, K[n^α, n^β] is the exchange-correlation energy, and n^σ is the electron density for spin σ, with α and β representing spin up and spin down, respectively.

The KS functional can be minimized with an SCF procedure⁸⁹ in which the KS-Fock matrix

is iteratively diagonalized using appropriate measures to ensure that the procedure becomes convergent. The first term of Eq. (2), H_core, consists of the matrix elements of the core Hamiltonian that arises from the first two terms of Eq. (1) as

where Z_A is the nuclear charge of atom A at position r_A. The last two terms in Eq. (2), J and K^σ, are the Coulomb and exchange-correlation matrix, respectively.

As only the two first terms in Eq. (1)—or equivalently, the first term of Eq. (2)—are necessary for exactness for one-electron systems,⁴⁵ the exact KS exchange-correlation functional must perfectly cancel out the Coulomb term for any one-electron density $n1σ$ with spin σ,

DFAs violate this condition, leading to a spurious self-interaction (SI) energy

The Perdew–Zunger self-interaction correction⁴⁵ (PZ-SIC) enforces the correct behavior of the DFA by explicitly removing the SI energy orbital by orbital, leading to a corrected total energy functional

where N^σ is the number of electrons with spin σ.

Because of this correction, the PZ-SIC Hamiltonian depends not only on the total spin densities but also on the individual orbital densities, at variance to standard KS-DFT. Due to this, the unitary invariance of KS-DFT⁸⁹ is broken, and the functional is no longer invariant to orbital rotations within the occupied space.

The optimal PZ-SIC orbitals minimize E_PZ. It has been shown that localized orbitals generally deliver lower self-interaction corrected total energies E_PZ than delocalized ones.^45,90,91 Hence, the best way to minimize Eq. (6) is to start out with localized orbitals (see Ref. 92 and references therein for discussion on this topic) produced by the Foster–Boys,⁹³ Edmiston–Ruedenberg,⁹⁴ or Pipek–Mezey⁹⁵ method or generalizations thereof ⁹⁶ and then iteratively rotate the orbitals until the derivative of the energy, the occupied–occupied orbital gradient,³⁸ 

vanishes; κ^σ = 0 is known as the Pederson condition.⁹¹ Here, $fiσ=J(piσ)+K(piσ)$ is the orbital Fock matrix (without the one-electron part) that arises from the density matrix $piσ$ corresponding to the ith occupied orbital of spin σ,

where $ciσ$ are the coefficients for the ith optimal orbital. In addition to the occupied–occupied rotations, the gradient for which was given in Eq. (7), the occupied-virtual rotations also have to be optimized, as discussed in Ref. 68.

Practical calculations based on orbital rotations, such as those with the erkale program, pursue minimization only to a finite numerical threshold; that is, until the orbital gradient is small but non-zero. Still, the minimization of E_PZ by orbital rotations is arduous due to the slow convergence of the orbital optimization. Moreover, as was already mentioned in the Introduction, the correct minimization of Eq. (6) has been shown to require complex-valued orbitals and to exhibit a plethora of local minima.⁶⁸ 

The asymptotic scaling of the orbital rotation approach is determined by the calculation of the orbital gradient in Eq. (7). Assuming K basis functions, c^σ and ${fiσ}i=1Nσ$ are all K × K matrices, with K ≥ N^σ. Disregarding the cost to compute the N^σ orbital Fock matrices (which should overall scale linearly with the system size in an optimal implementation due to the localized nature of the optimal orbitals), the evaluation of κ carries an iterative $O(NσK3)$ cost, while orbital optimization using the approach of Ref. 68 carries an iterative $O(K3)$ cost like conventional Kohn–Sham density functional theory. Still, the bottleneck in practical calculations is typically in the evaluation of the orbital Fock matrices, and alternative approaches that avoid the irrelevant virtual–virtual block can be pursued in the case K ≫ N^σ. However, due to the presence of many local minima and saddle point solutions, the approach of Ref. 68 recommends the use of stability analysis; since there are $O((Nσ)2)$ rotation angles, stability analysis employing iterative diagonalization carries a $O((Nσ)4)$ cost.

At variance to the direct minimization of E_PZ with orbital rotations, Pederson and co-workers proposed using Fermi–Löwdin orbitals (FLOs) together with PZ-SIC, giving rise to the FLO-SIC method.^9–12 FLOs were originally introduced by Luken and co-workers in a series of papers that dealt with the exchange or Fermi hole in the context of many-electron wave functions.^97–100 In the FLO approach, one starts by building Fermi orbitals (FOs) $ϕFO,iσ$ via

where $cFOσ$ and C^σ contain the FO and KS orbital coefficients, respectively. The transformation matrix R^σ is defined as

where $|aiσ$ denotes a position eigenstate localized at the so-called Fermi-orbital descriptor (FOD) $aiσ$⁠, and $ψjσ$ are KS orbitals. The FOs are normalized, but they do not form an orthonormal set in general. As a consequence, Luken proposed using Löwdin’s method of symmetric orthonormalization¹⁰¹ to end up with an orthonormalized set of Fermi–Löwdin orbitals (FLOs) $ϕkσ$ as

where c^σ holds the FLO coefficients, while T^σ and Q^σ contain the eigenvectors and eigenvalues of the FO overlap matrix, respectively.

The usefulness of the FOs and FLOs arises from the key property of the FO. It has the value $⟨aiσ|ϕFO,iσ⟩=nσ(aiσ)$ at its descriptor, meaning that all the electron density (for a given spin channel) at this point in space comes from a single orbital. Thus, the FOs are localized around their corresponding FODs, in contrast to the typically delocalized KS orbitals. As the orthonormalization mixes the FOs, the FLOs end up slightly less localized than the original non-orthogonal FOs.

Another key feature of the FLO approach is that since the FLOs are uniquely determined by Eqs. (10) and (11) for a given set of FODs $aiσ$ and a given set of occupied orbitals C^σ, the FLO-SIC approach turns out to restore the unitary invariance of KS-DFT: rotations of the occupied orbitals in C^σ are countermanded by an inverse rotation occurring in R^σ in Eq. (10). This feature means that the optimization of the FODs and that of the density matrix can be decoupled.

The FODs $aiσ$ formally introduced in Eq. (10) are a key feature of the FLO-SIC approach. In FLO-SIC, the optimization of the N^σ(N^σ − 1)/2 orbital rotation angles within the occupied space (twice that if imaginary rotations are also included as in the CSIC approach) is replaced by the optimization of 3N^σ FOD coordinates.^9–12,27 Although full minimization of the PZ-SIC energy functional used in FLO-SIC requires optimization of the FODs, the interesting part about the model is that since plausible FODs can be generated in an automatic fashion,⁸⁴ qualitative calculations may be performed without having to optimize the FODs, although the predictive power of such calculations is limited in the same way as the use of Foster–Boys orbitals for evaluating the self-interaction correction criticized in Ref. 38, as will be seen in Sec. III D.

Surprisingly, the optimization of this small number of FOD coordinates turns out to be at least as challenging as the optimization of the far more numerous possible orbital rotations. As previously shown in Refs. 10 and 11, the construction of the FOD derivatives from the orbital Fock matrices leads to $O((Nσ)4)$ scaling and $O((Nσ)3)$ data storage cost; this computational scaling and data storage behavior is also realized within pyflosic. The gradient of the energy with respect to the FODs can be expressed as^10,11

where the elements of the Lagrange multiplier matrix $εklk,σ$ are defined by

As always, the FOD optimization starts from some initial values for the FODs. As was discussed in Subsection II B, minimization of the PZ functional by orbital rotation techniques typically starts from localized orbitals. Now, since FLOs are highly localized around their FODs, this procedure for initialization can also be carried over to FLO-SIC by retrieving initial FODs from the centroids of the localized orbitals from the Foster–Boys, Edmiston–Ruedenberg, or (generalized) Pipek–Mezey.⁸⁴ Such a procedure is implemented in the python center of the mass (pycom) module within pyflosic, which employs the second-order orbital localizer¹⁷ in pyscf.

Several other ways to initialize the FODs have also been recently discussed in Ref. 84. They include an electronic force field that yields quasi-classical positions for the electrons to which the FODs can be assigned (pyeff); using Lewis-like bonding information to place FODs accordingly (pylewis) or a Thomson problem-like procedure for obtaining the FODs from Monte–Carlo minimization of a distribution of point charges under the restriction of a certain bond order (fodmc). The initial FODs from any of these alternative generators can be easily read in by pyflosic (see Ref. 84 for further details on automatic FOD initialization).

Regardless of the procedure used to initialize the FODs, the initial FOD geometries should always be visualized; an example will be discussed below in Sec. III. By visualizing the resulting electronic geometry, one should check whether it is in reasonable agreement with Lewis¹⁰² or Linnett double-quartet (LDQ) theory,^103–106 as good-natured FODs should generally correspond to these theories.^84,107

The gradients of these initial FODs are usually non-negligible, which is why FOD optimization is an important part of FLO-SIC calculations. FOD optimization is carried out automatically in pyflosic through an interface to the atomic simulation environment (ASE¹⁰⁸), which has ample algorithms for geometry optimization with forces, including conjugate gradients (CG), the (limited memory) Broyden–Fletcher–Goldfarb–Shanno scheme [(l-)bfgs],^109–116 and the fast inertial relaxation engine (fire).¹¹⁷ In order to use ase, the FOD gradient $∂ESIC/∂aiσ$ is simply expressed in terms of a fictitious force $fiσ=−∂EPZ/∂aiσ$ acting on the FODs, and the optimization is continued until some numerical threshold is reached for these FOD forces.

As was discussed in Subsection II C, the optimization of the energy functional used in FLO-SIC can be split into two subproblems: optimization of the density matrix for a fixed set of FODs and optimization of the FODs for a fixed density matrix, which was already discussed in Subsection II D. Although the PZ-SIC Hamiltonian $FPZσ$ based on Eq. (6) implies that each orbital experiences a different potential, all the orbital-dependent Hamiltonians can be cast in a common form by employing a unified Hamiltonian approach.^38,90,118 Assuming such a unified Hamiltonian, the SCF equations that determine the optimal occupied orbitals for PZ-SIC can be written as

where S is the basis set’s overlap matrix and E^σ is a diagonal matrix that holds the energy eigenvalues.

There are at least two different ways to construct a unified SIC Hamiltonian $FSICσ$⁠,^38,90,118 and both Hamiltonians presented here have been implemented in pyflosic. The first SIC Hamiltonian

only contains operators that project to and from the space of occupied orbitals, as indicated by the identifier OO. This is equivalent to ignoring the frequently small off-diagonal Lagrange multipliers that ensure orbital orthonormality.³⁸ The second SIC Hamiltonian is

where the virtual space projector is defined as

$FSIC,OOOVσ$ also allows for projections between the occupied and virtual orbital spaces, leading to the identifier OOOV; this Hamiltonian does not assume a diagonal Lagrange multiplier matrix. The second SIC Hamiltonian can also be written as a block matrix in the occupied and virtual spaces as

as the operator includes no terms that couple virtual orbitals together.

The features of pyflosic are first showcased using the tetracyanoethylene [C₂(CN)₄] molecule. In addition to specifying the nuclear geometry for the molecule, as is necessary for any electronic structure calculation, the initial FODs have to be specified in order to perform a FLO-SIC calculation, as was discussed in Subsection II D; that is, an electronic geometry needs to be specified as well.

Besides the versatile core functionality inherited from pyscf, pyflosic has unique features specific to FLO-SIC calculations. The first of these is the automatic generation of the electronic geometry, i.e., the FODs. The necessary code to obtain initial FODs with the pycom approach is presented in Fig. 1, and the resulting electronic geometry is visualized in Fig. 2. Note that only two lines of code are needed to generate a reasonable initial guess for the FODs, not counting the import of the necessary python modules.

Alternative to the use of pycom, the FODs can also be initialized in pyflosic by reading in the output of any of the other methods described in Ref. 84. The necessary code for the read-in is shown in Fig. 4 that will be discussed in more detail below.

Optimized FODs have been shown to carry bonding information in the context of Lewis/LDQ theory, which is discussed in detail in Ref. 84. For example, a simple FOD bond order (BO_FOD) can be evaluated via

just by counting the number of FODs m between the bonded atoms and then dividing by the number of atoms n that partake in this bond. The initial set of FODs shown in Fig. 2 appears to be reasonable, as it is in excellent agreement with Lewis/LDQ theory; it can thus be expected that the optimization of these FODs will yield useful results as well. It is, however, important to note that since the exchange-correlation functional affects the electron density, the reasonableness of the pycom guess may depend on the system and the functional, which is why we recommend to always visualize the initial FODs before running calculations. As the FODs change during the optimization, the final geometries should be inspected as well.

Having established the necessary nuclear and electronic geometries, the FLO-SIC calculation can begin. As discussed in Sec. II, FLO-SIC has two classes of degrees of freedom: the occupied orbitals, i.e., the electron density, and the FODs. Both should be optimized in order to minimize the PZ-SIC energy functional; the workflow for such FLO-SIC calculations is illustrated in Fig. 3. Note that the self-interaction correction is re-evaluated at every SCF iteration with the up-to-date density matrix, ensuring the self-consistency of the orbitals and the correction.

However, in some cases, one may want to keep the density or the FODs fixed; this may be useful for exploratory FLO-SIC calculations. The default mode in pyflosic is to optimize both the FODs and the density, but fixed-density and fixed-FOD calculations are also supported.

As was discussed in Sec. II, pyflosic optimizes the FODs via ase, whereas the electron density is optimized through a unified Hamiltonian approach. The example code for the optimization of both the density and the FODs with pyflosic is shown in Fig. 4. For this example, the FODs generated with pycom (see Fig. 2) are used as starting points for the FOD optimization, and the default unified SIC Hamiltonian $FSIC,OOOVσ$ is used for the density optimization within the SCF approach. Note that again, only two function calls are needed to run the calculation. FOD optimization can be analyzed just like any other geometry optimization in ase, e.g., by visualizing the trajectory of the system generated by ase.

Having the optimized FODs for a given level of theory (including the exchange-correlation functional, quadrature grid, and the orbital basis set), calculations for other levels of theory can easily be performed by starting the calculations from the preoptimized FODs. Because the cost of the FLO-SIC calculation is heavily dependent on the size of the orbital basis set, it is recommended to start out with preoptimized FODs from calculations using small basis sets before going to expensive calculations in large but more accurate basis sets.

Note, however, that if the level of theory is changed significantly, e.g., by going from an LDA to a GGA or a meta-GGA, the optimal FOD positions may also experience large changes, decreasing the value of the preoptimization. The changes in the optimal FODs are related not only to the effect of the total density changing with the exchange-correlation functional but also to the FLO single electron densities having much sharper features than the total density,¹²³ which accentuates the sensitivity to the DFA. FLO-SIC may thus be more sensitive to the exchange-correlation functional than KS-DFT. This is also reflected in the functional requirements for PZ-SIC, which differ from KS-DFT. For example, the Perdew–Burke–Ernzerhof functional with adjustments for accuracy for solids¹²⁰ (PBEsol) has been found to afford better accuracy for molecular PZ-SIC calculations than the original functional,¹²⁴ at variance to calculations at the KS-DFT level.^87,125

Having exemplified the use of the novel code, we proceed with a quantitative application. As was already mentioned in the Introduction, reproducing the atomization energies of SO₂ and SF₆ is known to require large and flexible basis sets at the Kohn–Sham level of theory,^30,31,33 making these molecules ideal candidates for a basis set convergence study of FLO-SIC. For simplicity, we chose to use fixed high-level ab initio geometries from the W4-17 database¹²⁷ for the molecules, as fixed nuclear geometries suffice for the present purposes of establishing convergence to the complete basis set limit. We also chose the polarization-consistent pc-n family of basis sets for this study, since these basis sets are designed for achieving optimal convergence to the basis set limit in DFT and Hartree–Fock calculations.^121,122 We furthermore chose to study the PBEsol functional, as it has been found to yield good accuracy in PZ-SIC calculations, as was already mentioned in Sec. III C;^87,125 however, basis set convergence patterns are well-known to be similar for both Hartree–Fock and all density functionals in the lack of post-Hartree–Fock correlation contributions.

However, as it is well-known that SIC methods require much larger quadrature grids than KS-DFT to reach similar levels of convergence,^38,62 a (200,590) quadrature grid was used for all calculations because preliminary tests revealed that such a grid was necessary to converge molecular FLO-SIC total energies of SF₆ and SO₂ roughly to μE_h accuracy. Unpruned grids are used in this work because at variance to KS-DFT, grid pruning leads to significant errors in FLO-SIC calculations: the orbital densities are not spherically symmetric near the nuclei in contrast to the total electron density used in KS-DFT.

For comparison, we also include results for the default basis sets used in nrlmol, i.e., the DFO and DFO+ basis sets for density functional calculations that have been used in a number of FLO-SIC studies in the literature.^{9–12,55,71–83,85,86} The DFO basis set is recommended for general use,²³ while the DFO+ basis set is obtained from the DFO basis by adding further polarization functions aimed to improve the accuracy of polarizability calculations;²³ naturally, the additional functions in DFO+ may also improve the convergence of other properties. For consistency with the previous literature, the DFO and DFO+ basis sets were extracted from NRLMOL for this work; the converted sets are available on github as part of pyflosic.

Four sets of computational models will be considered: the first one being Kohn–Sham DFT and the next three being variants of FLO-SIC. In the guess-pc-0 scheme, the FODs are fixed to the pycom starting guess values computed in the pc-0 basis set. Next, in the opt-pc-0 scheme, the FODs are frozen to the variationally optimized pc-0 values. Finally, in the opt-basis scheme, the FODs are fully optimized. All FOD optimizations are carried out until a force convergence criterion of $maxiσ|fiσ|<10−3Eh/a0$ (fmax = 0.0514 eV/Å) is satisfied.

The atomization energies resulting from the previously described procedures are shown in Table I. The data show remarkable variations of hundreds of kcal/mol when the size of the basis set is increased. Because larger basis sets are successively better at describing polarization effects in the molecules—thus decreasing the total energy of the molecule—the fully variational atomization energies (i.e., KS-DFT and the opt-basis model) increase with increasing basis set size. An acceptable level of convergence (basis set truncation error smaller than 1 kcal/mol) is only reached with the quadruple-ζ pc-3 basis set, highlighting the need to support large basis sets in FLO-SIC calculations that are routinely available in pyflosic. In contrast, the DFO and DFO+ basis sets produce results between the double-ζ pc-1 and the triple-ζ pc-2 basis sets, suggesting that the DFO and DFO+ basis sets are merely of polarized double-zeta quality. The DFO and DFO+ basis sets have a basis set truncation error of roughly 14 kcal/mol and 30 kcal/mol in KS-DFT calculations on SO₂ and SF₆, respectively, while in FLO-SIC calculations, the truncation error for SO₂ increases to 17 kcal/mol, the one for SF₆ remaining 30 kcal/mol. These results indicate that the DFO and DFO+ basis sets are woefully inaccurate for quantitative applications in chemistry, underlining the need to support high-angular momentum basis sets in KS-DFT and FLO-SIC calculations.

An examination into the three flavors of FLO-SIC studied in Table I shows that optimization of the FODs is clearly necessary, the guess-pc-0 atomization energies being very far from the fully optimized FLO-SIC values, guess-pc-0 overestimating the atomization energy of SF₆ by over 300 kcal/mol. In contrast, the atomization energies obtained with FODs frozen to the pc-0 optimal values are surprisingly close to the fully variational results, showing differences below 1 kcal/mol for all values except for the pc-3 value for SO₂ that differs from the fully variational one by 3.2 kcal/mol. These results suggest that using FODs frozen to preoptimized values for a small basis set may be an acceptable alternative to fully variational FLO-SIC calculations if the goal is to estimate the importance of self-interaction corrections to atomization energies. The preoptimized FODs also allow for straightforward basis set convergence studies, as omitting the optimization of the FODs implies a significant speed-up of FLO-SIC calculations in large basis sets.

As a further topic, the contraction error in the pc-n basis sets for DFT calculations and fully optimized FLO-SIC calculations was also studied, since the self-interaction correction may affect the core orbitals significantly. The atomization energies in the uncontracted pc-n (unc-pc-n) basis sets are shown in Table II. As shown by the comparison of the results in Tables I and II, the contraction errors are under 1 kcal/mol for all basis sets larger than pc-1, while the contraction error for the largest pc-3 and pc-4 basis sets is in the order of 0.4 kcal/mol for FLO-SIC calculations, which is (unsurprisingly) several times larger than the contraction error in the KS-DFT calculations. Uncontracting the basis set allows for an improved description of the core orbitals, and the unc-pc-4 results in Table II are our best estimates for the FLO-SIC atomization energies of SO₂ and SF₆ with the PBEsol functional. The W4-17 atomization energies for SO₂ and SF₆ are 260.580 kcal/mol and 485.425 kcal/mol, respectively. As can be seen from Table II, FLO-SIC reduces the error in the PBEsol functional from 42.5 kcal/mol and 83.0 kcal/mol for SO₂ and SF₆ at the KS-DFT level of theory to −36.3 kcal/mol and −33.0 kcal/mol at the FLO-SIC level of theory, respectively, in analogy to the findings for the PZ-SIC level of theory in Ref. 87.

We have presented the implementation of FLO-SIC in pyflosic as an extension to the pyscf code, as well as given instructive examples of code usage. We have also demonstrated the need to be able to use large and flexible basis sets in FLO-SIC calculations with a basis set convergence study of the atomization energies of SO₂ and SF₆, for which we showed that the DFO basis sets that have been commonly used in FLO-SIC calculations in the literature exhibit truncation errors of tens of kcal/mol.

The most important features of pyflosic are summarized in Fig. 5. Our implementation offers the possibility to use FLO-SIC with any GTO basis set, various types of quadrature grids, and hundreds of LDA, GGA, and meta-GGA exchange-correlation functionals available in the libxc and xcfun libraries. FLO-SIC depends on FODs, which need to be initialized and optimized. Importantly, the FODs can be initialized automatically in pyflosic with pycom—the internal routine based on localized orbitals’ centroids—or read in from other FOD generators, such as the pyeff, pylewis, and fodmc routines presented in Ref. 84. FOD optimization is carried out in pyflosic through an interface to ase. The electron density is optimized in pyflosic with an SCF procedure employing a unified Hamiltonian.

The version of pyflosic described in this work has already been used in some applications.^84,123 However, pyscf offers features that have not yet been exploited in FLO-SIC calculations. For instance, pyscf implements solvation models, such as the domain-decomposed conductor-like screening model^128–130 (COSMO) and the domain-decomposed polarizable continuum model^131,132 (PCM), which could be combined with FLO-SIC in pyflosic for studying systems in solution.

For the next iteration of pyflosic, we also aim to extend the program to include conventional RSIC and CSIC approaches based on orbital rotation techniques, which are presently available only in erkale. Since SIC is sensitive to the details of the evaluation of the density functional (e.g., the numerical quadrature), consistent implementations of the various schemes will allow for unbiased comparison of various SIC methods. In addition, FLO-SIC stability analysis along the lines of Ref. 68 will be investigated. Due to the simple interfacing required for SIC calculations and the modularity of pyflosic, additional interfaces to other electronic structure codes such as psi4⁸ could be provided in the future.

An important feature of pyscf is calculations on crystalline systems through the use of periodic boundary conditions.^7,14 Although calculations on periodic solids with exact exchange are tractable within a GTO basis set as in pyscf,¹³³ allowing the use of various hybrid functionals that are less prone to self-interaction error, exact exchange is undesirable in many cases in the study of the solid state due to problems with metallic systems. SIC does not pose such problems, and the application of SIC to solid-state systems has been an active field of research for several decades.^{45,64,67,134–157} pyflosic could be extended to solids in the future as well. The necessity of a localized orbital picture within a periodic system is a complication for pushing FLO-SIC to the solid state, requiring changes to the algorithms and to the initialization of FODs. Although orbital localization methods are less developed for the solid state than for molecular studies, solid-state FODs could be initialized in analogy to the present work either with maximally localized Foster–Boys¹⁵⁸ or generalized Pipek–Mezey Wannier functions.¹⁵⁹ The other methods of Ref. 84 might also be extended for the solid state; for instance, developmental versions of pyeff and fodmc are already able to provide FODs that appear reasonable for simple solids (e.g., solid deuterium¹⁶⁰ or lithium) and metal–organic frameworks (MOFs).^161–165 

The pyflosic code, presented and used within this study, is openly available on github (https://github.com/pyflosic/pyflosic) and can be referenced via https://doi.org/10.5281/zenodo.3948143.

J. Kortus and S. Schwalbe were funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project No. 421663657—KO 1924/9-1. J. Kraus and J. Kortus were funded by the DFG—Project No. 169148856—SFB 920, subproject A04. K. Trepte was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, as part of the Computational Chemical Sciences Program under Award No. DE-SC0018331. S. Lehtola was supported by the Academy of Finland (Suomen Akatemia) through Project No. 311149. We also thank the ZIH in Dresden for computational time and support, Professor Mark R. Pederson for creating the FLO-SIC method that is the core of everything discussed within the manuscript, and Dr. Torsten Hahn for various discussions and his contributions in an earlier stage of the project.