PyFLOSIC: Python-based Fermi-L\"owdin orbital self-interaction correction

TL;DR

PyFLOSIC is an open-source Python implementation of the Fermi-L"owdin orbital self-interaction correction, enabling flexible and efficient quantum chemistry calculations with various basis sets and functionals.

Contribution

It introduces a modular, Python-based implementation of FLO-SIC integrated with PySCF, supporting automatic initialization and optimization of Fermi-orbital descriptors.

Findings

Supports various basis sets and functionals within PySCF

Enables automatic initialization and optimization of descriptors

Facilitates applications and further development of FLO-SIC

Abstract

We present PyFLOSIC, an open-source, general-purpose Python implementation of the Fermi-L\"owdin orbital self-interaction correction (FLO-SIC), which is based on the Python simulation of chemistry frame-work (PySCF) electronic structure and quantum chemistry code. Thanks to PySCF, PyFLOSIC can be used with any kind of Gaussian-type basis set, various kinds of radial and angular quadrature grids, and all exchange-correlation functionals within the local density approximation (LDA), generalized-gradient approximation (GGA), and meta-GGA provided in the Libxc and XCFun libraries. A central aspect of FLO-SIC are Fermi-orbital descriptors, which are used to estimate the self-interaction correction. Importantly, they can be initialized automatically within PyFLOSIC and optimized with an interface to the atomic simulation environment, a Python library which provides a variety of powerful gradient-based algorithms for geometry optimization. Although PyFLOSIC has already facilitated applications of FLO-SIC to chemical studies, it offers an excellent starting point for further developments in FLO-SIC approaches, thanks to its use of a high-level programming language and pronounced modularity.

Tables

Table 1. Table 1: Basis set convergence for the atomization energy of SO 2 and SF 6 in kcal/mol , calculated with an unpruned (200,590) quadrature grid and the PBEsol exchange-correlation functional. W4-17 Karton, Sylvetsky, and Martin ( 2017 ) nuclear geometries were used.

	KS-DFT		FLO-SIC (guess-pc-0)		FLO-SIC (opt-pc-0)		FLO-SIC (opt-basis)
basis	SO2	SF6	SO2	SF6	SO2	SF6	SO2	SF6
pc-0	155.876	396.551	298.089	574.227	56.087	241.725	56.087	241.725
pc-1	268.312	512.775	412.687	729.099	184.844	388.415	185.488	388.740
DFO	289.147	538.089	434.999	763.010	207.368	422.575	207.406	422.623
DFO+	289.261	538.638	435.792	763.500	207.205	423.141	207.626	423.183
pc-2	294.159	556.651	440.800	786.006	215.682	442.813	214.234	443.019
pc-3	302.640	568.079	448.222	792.663	220.803	452.028	224.012	452.425
pc-4	303.179	568.534	448.727	792.518	224.647	452.562	224.629	452.898

Table 2. Table 2: Atomization energies of SO 2 and SF 6 in kcal/mol with uncontracted basis sets, calculated with an unpruned (200,590) quadrature grid and the PBEsol exchange-correlation functional. W4-17 nuclear geometries were used.

	KS-DFT		FLO-SIC (opt-basis)
basis	SO2	SF6	SO2	SF6
unc-pc-0	158.414	400.937	59.557	245.947
unc-pc-1	269.065	515.849	186.428	395.124
unc-pc-2	294.361	557.464	215.115	443.889
unc-pc-3	302.614	568.003	224.295	452.023
unc-pc-4	303.126	568.411	224.235	452.456

Equations

$$E_{\text{KS}}[n^{\alpha },n^{\beta }]=T_{\text{s}}[n^{\alpha },n^{\beta }]+V[n]+J[n]+K[n^{\alpha },n^{\beta }],$$

$${\mathbf{F}}_{\text{KS}}^{\sigma }={\mathbf{H}}_{\text{core}}+\mathbf{J}+{\mathbf{K}}^{\sigma },$$

$${\mathbf{H}}_{\text{core}}=-\frac{1}{2}{\nabla }^2-\sum _A^{\text{nuclei}}\frac{Z_A}{|\text{r}-{\text{r}}_A|},$$

$$K[n_1^{\sigma },0]=-J[n_1^{\sigma }].$$

$$E_{\text{SI}}[n_1^{\sigma }]=K[n_1^{\sigma },0]+J[n_1^{\sigma }].$$

$$E_{\text{PZ}}$$

$$=E_{\text{KS}}[n^{\alpha },n^{\beta }]-\sum _{\sigma }\sum _i^{N^{\sigma }}{E}_{\text{SI}}[n_i^{\sigma }],$$

$${\kappa }_{ij}^{\sigma }=({\text{c}}_j^{\sigma }{)}^{\text{T}}({\mathbf{f}}_i^{\sigma }-{\mathbf{f}}_j^{\sigma }){\text{c}}_i^{\sigma }$$

$${\mathbf{p}}_i^{\sigma }={\text{c}}_i^{\sigma }({\text{c}}_i^{\sigma }{)}^{\text{T}},$$

$${\mathbf{c}}_{\text{FO}}^{\sigma }={\mathbf{C}}^{\sigma }{\mathbf{R}}^{\sigma },$$

$$R_{ji}^{\sigma }=\frac{⟨{\psi }_j^{\sigma }⟩{\text{a}}_i^{\sigma }}{\sqrt{n^{\sigma }({\mathbf{a}}_i^{\sigma })}}.$$

$${\mathbf{c}}^{\sigma }={\mathbf{c}}_{\text{FO}}^{\sigma }[{\mathbf{T}}^{\sigma }({\mathbf{Q}}^{\sigma }{)}^{-1/2}({\mathbf{T}}^{\sigma }{)}^{\text{T}}],$$

$$\frac{\partial E_{\text{PZ}}}{\partial {\text{a}}_i^{\sigma }}=\frac{\partial E_{\text{SIC}}}{\partial {\text{a}}_i^{\sigma }}=\sum _{k,l}^{N^{\sigma }}{\epsilon }_{kl}^{k,\sigma }\left(⟨\frac{\partial {\varphi }_k^{\sigma }}{\partial {\text{a}}_i^{\sigma }}⟩{\varphi }_l^{\sigma }+⟨{\varphi }_l^{\sigma }⟩\frac{\partial {\varphi }_k^{\sigma }}{\partial {\text{a}}_i^{\sigma }}\right),$$

$${\epsilon }_{ij}^{k,\sigma }=-({\text{c}}_i^{\sigma }{)}^{\text{T}}{\mathbf{f}}_k^{\sigma }{\text{c}}_j^{\sigma }.$$

$${\mathbf{F}}_{\text{PZ}}^{\sigma }{\mathbf{C}}^{\sigma }=({\mathbf{F}}_{\text{KS}}^{\sigma }+{\mathbf{F}}_{\text{SIC}}^{\sigma }){\mathbf{C}}^{\sigma }=\mathbf{S}{\mathbf{C}}^{\sigma }{\mathbf{E}}^{\sigma },$$

$${\mathbf{F}}_{\text{SIC,OO}}^{\sigma }=-\frac{1}{2}\sum _i^{N^{\sigma }}({\mathbf{f}}_i^{\sigma }{\mathbf{p}}_i^{\sigma }\mathbf{S}+\mathbf{S}{\mathbf{p}}_i^{\sigma }{\mathbf{f}}_i^{\sigma })$$

$${\mathbf{F}}_{\text{SIC,OOOV}}^{\sigma }=-\mathbf{S}\sum _i^{N^{\sigma }}({\mathbf{p}}_i^{\sigma }{\mathbf{f}}_i^{\sigma }{\mathbf{p}}_i^{\sigma }+{\mathbf{v}}^{\sigma }{\mathbf{f}}_i^{\sigma }{\mathbf{p}}_i^{\sigma }+{\mathbf{p}}_i^{\sigma }{\mathbf{f}}_i^{\sigma }{\mathbf{v}}^{\sigma })\mathbf{S},$$

$${\mathbf{v}}^{\sigma }=\sum _i^{\text{virtual}}{\text{C}}_i^{\sigma }({\text{C}}_i^{\sigma }{)}^{\text{T}}.$$

$${\mathbf{F}}_{\text{SIC,OOOV}}^{\sigma }$$

$$=\left(\begin{array}{cc}{\mathbf{F}}_{\text{SIC,OO}}^{\sigma }& {\mathbf{F}}_{\text{SIC,OV}}^{\sigma }\\ {\mathbf{F}}_{\text{SIC,VO}}^{\sigma }& 0\end{array}\right),$$

$${\text{BO}}_{\text{FOD}}=m_{\text{electrons}}/n_{\text{centers}}$$

Full text

I Introduction

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 sub-problemsSun (2015); Kim et al. (2018); Lehtola et al. (2018); Herbst et al. (2019); Koval, Barbry, and Sánchez-Portal (2019); Iskakov et al. (2019) that can be easily combined to form a code that is greater than the sum of its parts; see e.g. 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).Pederson, Ruzsinszky, and Perdew (2014); Pederson (2015); Pederson and Baruah (2015); Yang, Pederson, and Perdew (2017) The core routines of PyFLOSIC were developed during Lenz Fiedler’s master’s thesis,Fiedler (2018) and the code is freely available on GitHub (https://github.com/pyflosic/pyflosic). PyFLOSIC builds on the Python simulation of chemistry framework (PySCF),Sun et al. (2018, 2020) 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.Sun et al. (2018, 2020) For example, self-consistent field (SCF) approaches like Hartree–Fock and density functional theoryHohenberg and Kohn (1964); Kohn and Sham (1965) (DFT) are available with various choices for spin treatment, density fitting techniques, as well as numerical quadrature methods. Second-order orbital optimization is also available in PySCF.Sun (2016)

Since PySCF is highly modular, and its parts can be imported like any other Python module, it is rather easy to add new functionality to it or to combine it with existing software,Sun et al. (2018) as has already been demonstrated by various authors.Kim et al. (2018); Koval, Barbry, and Sánchez-Portal (2019); Iskakov et al. (2019) 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).Pederson and Jackson (1990); Pederson, Jackson, and Pickett (1991); Pederson and Jackson (1991); Perdew et al. (1992); Porezag and Pederson (1996); Porezag (1997); Porezag and Pederson (1999); Kortus and Pederson (2000); Pederson et al. (2000) There is also a FLO-SIC implementation in a developer branch of the NWChem program,Aquino and Wong (2018); Aquino, Shinde, and Wong (2020) but to the best of our knowledge this implementation is not currently publicly available.

Similarly to NRLMOL and NWChem, PySCF—and therefore PyFLOSIC—uses a Gaussian-type orbital (GTO) basis set.Sun et al. (2018, 2020) 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 as well as effective core potential GTO basis sets are already available within PySCF, and additional ones can be downloaded from, e.g., the Basis Set ExchangePritchard et al. (2019) and parsed in from several formats.

The importance of the tractability of calculations within large basis sets is underlined by recent fully numerical benchmark studies that have shown that the reproduction of atomization energies even within DFT may require several shells of polarization functions.Jensen et al. (2017); Jensen (2017); Feller and Dixon (2018); Lehtola (2020a) PySCF contains fast routines for the evaluation of molecular integrals via the Libcint library,Sun (2015) which are also used in PyFLOSIC. Moreover, PyFLOSIC inherits OpenMP parallelization from PySCF, enabling the efficient use of multi-core computation architectures.Sun et al. (2020) 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 of even molecules for which very large and flexible basis sets are required, such as \ceSO2 and \ceSF6,Jensen et al. (2017); Jensen (2017); Lehtola (2020a) as will also be demonstrated later in this work. Recently developed powerful approaches to handle significant linear dependencies in the underlying molecular basis setLehtola (2019, 2020b) are also available in PyFLOSIC through PySCF, enabling accurate FLO-SIC calculations even in pathologically over-complete basis sets.Lehtola (2019, 2020b, 2020a)

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 LibxcLehtola et al. (2018) and XCFunEkström et al. (2010) 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, non-local correlation functionals, and range-separated hybrids.Sun et al. (2020) 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.fun 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 like PySCF,Sun et al. (2018) 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.Sun et al. (2018) Calculations can also be done interactively in the Python interpreter shell.Sun et al. (2018)

Having introduced and motivated PyFLOSIC, we will next discuss the theory of FLO-SIC in detail in Section II. We will start out by motivating the use of self-interaction corrections in general (Subsection II.1), continue by presenting Perdew and Zunger’s approach in specific (Subsection II.2), introduce Fermi–Löwdin orbitals to undertake the self-interaction correction following Perdew and Zunger’s prescription (Subsection II.3), and discuss the mandatory initialization and optimization of the Fermi-orbital descriptors that parametrize the Fermi–Löwdin orbitals (Subsection II.4). The schemes used in PyFLOSIC for optimizing the FLO-SIC density are discussed in Subsection II.5. The code is showcased in Section 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.1), and then the electron density and the FODs are optimized (Subsection III.2). A discussion on repeated calculations follows in Subsection III.3. As a practical application of the code, in Subsection III.4 we perform an in-depth basis set convergence study of the atomization energies of \ceSO2 and \ceSF6 that have been found to be challenging cases in the literature, as was already discussed above. The article concludes in a summary and outlook in Section IV. Atomic units are used throughout the paper, if not stated otherwise.

II Fermi–Löwdin orbital self-interaction correction

We will adopt the notation introduced by Lehtola and Jónsson (2014) 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 ($\bm{R}$).

II.1 Self-interaction correction in density functional theory

DFT has become one of the standard methods in computational materials science, condensed matter physics, as well as chemistry thanks to its combination of reasonable accuracy with computational efficiency.Sherrill (2010); Burke (2012) However, currently available density functional approximations (DFAs) are well-known to fail in a number of situations.Cohen, Mori-Sánchez, and Yang (2008); Perdew et al. (2009); Burke (2012) 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 \ceF- 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 been long 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 \ceH + H2 <-> H2 + H already poses a challenge, as many functionals overestimate the stability of the intermediate \ceH3 state.Johnson et al. (1994) 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),Kraisler and Kronik (2013) 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,Perdew and Zunger (1981) 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).Pederson and Lin (1988); Vydrov and Scuseria (2005); Klüpfel, Klüpfel, and Jónsson (2011); Gudmundsdóttir et al. (2013); Borghi et al. (2014); Gudmundsdóttir et al. (2014); Lehtola and Jónsson (2014); Pederson and Baruah (2015); Perdew et al. (2015); Cheng et al. (2016); Zhang, Weber, and Jónsson (2016); Schwalbe et al. (2018) 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 $\text{IP}=-\varepsilon_{\text{HOMO}}$; in contrast, the so-called Koopmans-compliant (KC) functionalsBorghi et al. (2014, 2015) can deliver especially accurate results for IPs.

Among the multitude of current implementations of SIC, most employ real-valued orbitals (RSIC).Garza, Nichols, and Dixon (2000); Garza et al. (2001); Patchkovskii, Autschbach, and Ziegler (2001); Patchkovskii and Ziegler (2002a, b); Vydrov and Scuseria (2004); Vydrov et al. (2006); Pemmaraju et al. (2007) However, it has been shown that SIC based on complex-valued orbitals (CSIC) has several advantages;Klüpfel, Klüpfel, and Jónsson (2011, 2012); Valdes et al. (2012); Gudmundsdóttir, Jónsson, and Jónsson (2015); Cheng et al. (2016) in fact, Lehtola, Head-Gordon, and Jónsson (2016) 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.Lehtola et al. (2012); S. Lehtola. ERKALE – HF/DFT from Hel (2016))

Another variant of SIC is the FLO-SIC approach mentioned in Section I, which has been the subject of various studies during recent years,Pederson, Ruzsinszky, and Perdew (2014); Hahn et al. (2015); Pederson and Baruah (2015); Pederson (2015); Pederson et al. (2016); Hahn et al. (2017); Yang, Pederson, and Perdew (2017); Kao et al. (2017); Sharkas et al. (2018); Schwalbe et al. (2018); Joshi et al. (2018); Withanage et al. (2018); Aquino and Wong (2018); Johnson et al. (2019); Jackson et al. (2019); Zope et al. (2019); Withanage et al. (2019); Santra and Perdew (2019); Trepte et al. (2019); Schwalbe et al. (2019); Yamamoto et al. (2019); Aquino, Shinde, and Wong (2020); Vargas et al. (2020) and is also the topic of the present work; it will be discussed in depth below in Subsections II.3, II.4, and II.5. It is first worth mentioning, however, that in spite of the 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.Vydrov and Scuseria (2004); Borghi et al. (2014); Lehtola, Jónsson, and Jónsson (2016) This is the so-called "paradox of self-interaction correction" discussed in the comprehensive summary of Perdew et al. (2015), which remains a challenge for the future. For instance, it has been pointed out that the Perdew–Zunger approachPerdew and Zunger (1981) does not completely eliminate one-electron self-interaction;Lundin and Eriksson (2001) 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.

II.2 Perdew–Zunger self-interaction correction

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

[TABLE]

where $E_{\text{KS}}[n^{\alpha},n^{\beta}]$ is the total KS energy, $T_{\text{s}}[n^{\alpha},n^{\beta}]$ is the kinetic energy of the non-interacting system, $V[n]$ is the external potential energy, $J[n]$ is the Coulomb energy, $K[n^{\alpha},n^{\beta}]$ is the exchange-correlation energy, and $n^{\sigma}$ is the electron density for spin $\sigma$, with $\alpha$ and $\beta$ representing spin up and spin down, respectively.

The KS functional can be minimized with a SCF procedure,Lehtola, Blockhuys, and Van Alsenoy (2020) in which the KS-Fock matrix

[TABLE]

is iteratively diagonalized, using appropriate measures to ensure the procedure becomes convergent. The first term of equation (2), $\bm{H}_{\text{core}}$, is the core Hamiltonian that arises from the first two terms of equation (1) as

[TABLE]

where $Z_{A}$ is the nuclear charge of atom $A$ at position $\textbf{r}_{A}$. The last two terms in equation (2), $\bm{J}$ and $\bm{K}^{\sigma}$, are the Coulomb and exchange-correlation matrix, respectively.

As only the two first terms in equation (1)—or equivalently, the first term of equation (2)—are necessary for exactness for one-electron systems,Perdew and Zunger (1981) the exact KS exchange-correlation functional must perfectly cancel out the Coulomb term for any one-electron density $n_{1}^{\sigma}$ with spin $\sigma$

[TABLE]

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

[TABLE]

The Perdew–Zunger self-interaction correctionPerdew and Zunger (1981) (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

[TABLE]

where $N^{\sigma}$ is the number of electrons with spin $\sigma$.

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-DFTLehtola, Blockhuys, and Van Alsenoy (2020) is broken, and the functional is no longer invariant to orbital rotations within the occupied space.

The optimal PZ-SIC orbitals minimize $E_{\text{PZ}}$. It has been shown that localized orbitals generally deliver lower self-interaction corrected total energies $E_{\text{PZ}}$ than delocalized ones.Perdew and Zunger (1981); Harrison, Heaton, and Lin (1983); Pederson, Heaton, and Lin (1984) Hence, the best way to minimize equation (6) is to start out with localized orbitals (see ref. 92 and references therein for discussion on this topic) produced, e.g., by the Foster–Boys,Boys (1960) Edmiston–Ruedenberg, Edmiston and Ruedenberg (1963) or Pipek–Mezey Pipek and Mezey (1989) method or generalizations thereof,Lehtola and Jónsson (2014) and then iteratively rotate the orbitals until the derivative of the energy, the occupied-occupied orbital gradient Lehtola and Jónsson (2014)

[TABLE]

vanishes; $\bm{\kappa}^{\sigma}=\bm{0}$ is known as the Pederson condition. Here, $\bm{f}_{i}^{\sigma}=\bm{J}(\bm{p}_{i}^{\sigma})+\bm{K}(\bm{p}_{i}^{\sigma})$ is the orbital Fock matrix (without the one-electron part) that arises from the density matrix $\bm{p}_{i}^{\sigma}$ corresponding to the $i$:th occupied orbital of spin $\sigma$

[TABLE]

where $\textbf{c}_{i}^{\sigma}$ are the coefficients for the $i$:th optimal orbital. In addition to the occupied-occupied rotations, the gradient for which was given in equation (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 still non-zero. Still, the minimization of $E_{\text{PZ}}$ by orbital rotations is arduous due to slow convergence of the orbital optimization. Moreover, as was already mentioned in the Introduction, the correct minimization of equation (6) has been shown to require complex-valued orbitals and to exhibit a plethora of local minima.Lehtola, Head-Gordon, and Jónsson (2016)

The asymptotic scaling of the orbital rotation approach is determined by the calculation of the orbital gradient in equation (7). Assuming $K$ basis functions, $\mathbf{c}^{\sigma}$ and $\{\mathbf{f}_{i}^{\sigma}\}_{i=1}^{N^{\sigma}}$ are all $K\times K$ matrices, with $K\geq N^{\sigma}$. Disregarding the cost to compute the $N^{\sigma}$ orbital Fock matrices (which should overall scale linearly with system size in an optimal implementation due to the localized nature of the optimal orbitals), the evaluation of $\bm{\kappa}$ carries an iterative $\mathcal{O}(N^{\sigma}K^{3})$ cost, while orbital optimization using the approach of ref. 68 carries an iterative $\mathcal{O}(K^{3})$ 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\gg N^{\sigma}$. 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 $\mathcal{O}((N^{\sigma})^{2})$ rotation angles, stability analysis employing iterative diagonalization carries a $\mathcal{O}((N^{\sigma})^{4})$ cost.

II.3 Fermi–Löwdin orbitals

At variance to the direct minimization of $E_{\text{PZ}}$ with orbital rotations, Pederson and coworkers proposed using Fermi–Löwdin orbitals (FLOs) together with PZ-SIC, giving rise to the FLO-SIC method.Pederson, Ruzsinszky, and Perdew (2014); Pederson (2015); Pederson and Baruah (2015); Yang, Pederson, and Perdew (2017) FLOs were originally introduced by Luken and coworkers in a series of papers that dealt with the exchange or Fermi hole in the context of many-electron wave functions.Luken and Culberson (1982); Luken and Beratan (1982); Luken (1984); Luken and Culberson (1984) In the FLO approach, one starts by building Fermi orbitals (FOs) $\phi_{\text{FO},i}^{\sigma}$ via

[TABLE]

where $\bm{c}_{\text{FO}}^{\sigma}$ and $\bm{C}^{\sigma}$ contain the FO and KS orbital coefficients, respectively. The transformation matrix $\bm{R}^{\,\sigma}$ is defined as

[TABLE]

where $|\textbf{a}_{i}^{\sigma}\rangle$ denotes a position eigenstate localized at the so-called Fermi-orbital descriptor (FOD) $\textbf{a}_{i}^{\sigma}$, and $\psi_{j}^{\sigma}$ 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 Löwdin (1950) to end up with an orthonormalized set of Fermi–Löwdin orbitals (FLOs) $\phi_{k}^{\sigma}$ as

[TABLE]

where $\bm{c}^{\sigma}$ holds the FLO coefficients, while $\bm{T}^{\sigma}$ and $\bm{Q}^{\sigma}$ contain the eigenvectors and eigenvalues of the FO overlap matrix, respectively.

The usefulness of the FOs and FLOs arises from the property of the FO. It has the value $\braket{\textbf{a}^{\sigma}_{i}}{\phi_{\text{FO},i}^{\sigma}}=\sqrt{n^{\sigma}({\bf a}_{i}^{\sigma})}$ 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 are 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 equations (10) and (11) for a given set of FODs $\bm{a}_{i}^{\sigma}$ and a set of occupied orbitals $\bm{C}^{\sigma}$, the FLO-SIC approach turns out to restore the unitary invariance of KS-DFT: rotations of the occupied orbitals in $\bm{C}^{\sigma}$ are countermanded by an inverse rotation occurring in $\bm{R}^{\sigma}$ in equation (10). This feature means that the optimization of the FODs and that of the density matrix can be decoupled.

II.4 Initialization and optimization of Fermi-orbital descriptors

The FODs $\textbf{a}^{\sigma}_{i}$ formally introduced in equation (10) are a key feature of the FLO-SIC approach. In FLO-SIC, the optimization of the $N^{\sigma}(N^{\sigma}-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^{\sigma}$ FOD coordinates.Pederson, Ruzsinszky, and Perdew (2014); Pederson (2015); Pederson and Baruah (2015); Yang, Pederson, and Perdew (2017); Aquino and Wong (2018) 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,Schwalbe et al. (2019) 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, e.g., Foster–Boys orbitals for evaluating the self-interaction correction criticized in ref. 38.

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 $\mathcal{O}((N^{\sigma})^{4})$ scaling and $\mathcal{O}((N^{\sigma})^{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 Pederson (2015); Pederson and Baruah (2015)

[TABLE]

where the elements of the Lagrange multiplier matrix $\varepsilon^{k,\sigma}_{kl}$ are defined by

[TABLE]

As always, the FOD optimization starts from some initial values for the FODs. As was discussed in Subsection II.2, minimization of the PZ functional by orbital rotation techniques is typically started from localized orbitals. Now, since FLOs are highly localized around their FODs, this 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.Schwalbe et al. (2019) Such a procedure is implemented in the Python center of mass (PyCOM) module within PyFLOSIC, which employs the second-order orbital localizerSun (2016) in PySCF.

Several other ways to initialize the FODs have also been recently discussed in ref. 84. They include an electronic force field which 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 of 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. Please refer to 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 Section III. By visualizing the resulting electronic geometry, one should check whether it is in reasonable agreement with Lewis Lewis (1916) or Linnett double-quartet (LDQ) theory,Linnett (1960, 1961, 1964); Luder (1964) as good-natured FODs should generally correspond to these theories.Kraus (2017); Schwalbe et al. (2019)

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 (ASELarsen et al. (2017)), which has ample algorithms for geometry optimization with forces, including conjugate gradients (CG), the (limited memory) Broyden–Fletcher–Goldfarb–Shanno scheme ((L-)BFGS),Broyden (1970); Fletcher (1970); Goldfarb (1970); Shanno (1970); Nocedal (1980); Liu and Nocedal (1989); Byrd et al. (1995); Zhu et al. (1997) and the fast inertial relaxation engine (FIRE).Bitzek et al. (2006) In order to use ASE, the FOD gradient $\partial E_{\text{SIC}}/\partial\,\textbf{a}^{\sigma}_{i}$ is simply expressed in terms of a fictitious force $\textbf{f}^{\,\sigma}_{i}=-\partial E_{\text{PZ}}/\partial\,\textbf{a}^{\sigma}_{i}$ acting on the FODs, and the optimization is continued until some numerical threshold is reached for these FOD forces.

II.5 Optimization of the density with FLO-SIC

As was discussed in Subsection II.3, the optimization of the energy functional used in FLO-SIC can be split into two sub-problems: 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.4. Although the PZ-SIC Hamiltonian $\bm{F}^{\sigma}_{\text{PZ}}$ based on equation (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.Heaton, Harrison, and Lin (1982); Harrison, Heaton, and Lin (1983); Lehtola and Jónsson (2014) Assuming such a unified Hamiltonian, the SCF equations that determine the optimal occupied orbitals for PZ-SIC can be written as

[TABLE]

where $\bm{S}$ is the basis set’s overlap matrix, and $\bm{E}^{\sigma}$ is a diagonal matrix that holds the energy eigenvalues.

There are at least two different ways to construct a unified SIC Hamiltonian $\bm{F}^{\sigma}_{\text{SIC}}$,Heaton, Harrison, and Lin (1982); Harrison, Heaton, and Lin (1983); Lehtola and Jónsson (2014) and both Hamiltonians presented here have been implemented in PyFLOSIC. The first SIC Hamiltonian

[TABLE]

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.Lehtola and Jónsson (2014) The second SIC Hamiltonian is

[TABLE]

where the virtual space projector is defined as

[TABLE]

Therefore, $\bm{F}_{\text{SIC,OOOV}}^{\sigma}$ 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

[TABLE]

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

III Example calculations with PyFLOSIC

The features of PyFLOSIC are first showcased using the tetracyanoethylene (\ceC2(CN)4) 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.4; that is, an electronic geometry needs to be specified as well.

III.1 Automatic generation of electronic geometry

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 figure 1, and the resulting electronic geometry is visualized in figure 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.

Alternatively 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 figure 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${}_{\text{FOD}}$) can be evaluated via

[TABLE]

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 figure 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.

III.2 FOD and density optimization

Having established the necessary nuclear and electronic geometries, the FLO-SIC calculation can begin. As was discussed in Section 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 figure 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, e.g., 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 Section II, PyFLOSIC optimizes the FODs via ASE, whereas the electron density is optimized through a unified Hamiltonian approach. Example code for the optimization of both the density and the FODs with PyFLOSIC is shown in figure 4. For this example, the FODs generated with PyCOM (see figure 2) are used as starting points for the FOD optimization, and the default unified SIC Hamiltonian $\bm{F}^{\sigma}_{\text{SIC,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.

III.3 Repeated calculations

Having the optimized FODs for a given level of theory (including, e.g., 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,Shahi et al. (2019) 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 solidsPerdew et al. (2008) (PBEsol) has been found to afford better accuracy for molecular PZ-SIC calculations than the original functional,Perdew, Burke, and Ernzerhof (1996) at variance to calculations at the KS-DFT level.Jónsson, Lehtola, and Jónsson (2015); Lehtola, Jónsson, and Jónsson (2016)

III.4 Basis set convergence of the atomization energy of \ceSO2 and \ceSF6

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 \ceSO2 and \ceSF6 is known to require large and flexible basis sets at the Kohn–Sham level of theory,Jensen et al. (2017); Jensen (2017); Lehtola (2020a) 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 databaseKarton, Sylvetsky, and Martin (2017) 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.Jensen (2001, 2002) 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 the Introduction;Jónsson, Lehtola, and Jónsson (2015); Lehtola, Jónsson, and Jónsson (2016) 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,Vydrov and Scuseria (2004); Lehtola and Jónsson (2014) 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 \ceSF6 and \ceSO2 to roughly $\mu E_{\text{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 which have been used in a number of FLO-SIC studies in the literature.Pederson, Ruzsinszky, and Perdew (2014); Hahn et al. (2015); Pederson and Baruah (2015); Pederson (2015); Pederson et al. (2016); Hahn et al. (2017); Yang, Pederson, and Perdew (2017); Kao et al. (2017); Sharkas et al. (2018); Schwalbe et al. (2018); Joshi et al. (2018); Withanage et al. (2018); Johnson et al. (2019); Jackson et al. (2019); Zope et al. (2019); Withanage et al. (2019); Santra and Perdew (2019); Trepte et al. (2019); Yamamoto et al. (2019); Vargas et al. (2020) The DFO basis set is recommended for general use,Porezag (1997) while the DFO+ basis set is obtained from the DFO basis by adding further polarization functions aimed to improve the accuracy of polarizability calculations;Porezag (1997) naturally, the additional functions in DFO+ may also improve the convergence of other properties. For consistency with 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. Last, in the opt-basis scheme, the FODs are fully optimized. All FOD optimizations are carried out until a force convergence criterion of $\max_{i\sigma}|\textbf{f}^{\,\sigma}_{i}|<10^{-3}~{}E_{\text{h}}/a_{0}$ is satisfied.

The atomization energies resulting from the previously described procedures are shown in table 1. 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-$\zeta$ pc-3 basis set, highlighting the need to support large basis sets in FLO-SIC calculations which is routinely available in PyFLOSIC. In contrast, the DFO and DFO+ basis sets produce results between the double-$\zeta$ pc-1 and the triple-$\zeta$ pc-2 basis sets, suggesting 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 \ceSO2 and \ceSF6, respectively, while in FLO-SIC calculations the truncation error for \ceSO2 increases to 17 kcal/mol, the one for \ceSF6 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 as well as FLO-SIC calculations.

An examination into the three flavors of FLO-SIC studied in table 1 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 \ceSF6 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 \ceSO2 that differs from the fully variational one by 3.2 kcal/mol. These results suggest 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 e.g. 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 were 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 2. As shown by the comparison of the results in tables 1 and 2, 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 2 are our best estimates for the FLO-SIC atomization energies of \ceSO2 and \ceSF6 with the PBEsol functional. The W4-17 atomization energies for \ceSO2 and \ceSF6 are 260.580 kcal/mol and 485.425 kcal/mol. As can be seen from table 2, FLO-SIC reduces the error in the PBEsol functional from 42.5 kcal/mol and 83.0 kcal/mol for \ceSO2 and \ceSF6 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, analogously to the findings for the PZ-SIC level of theory in ref. 87.

IV Summary and Outlook

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 \ceSO2 and \ceSF6, 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 figure 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.Shahi et al. (2019); Schwalbe et al. (2019) 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 modelCancès, Maday, and Stamm (2013); Lipparini et al. (2013, 2014) (COSMO) as well as the domain-decomposed polarizable continuum modelStamm et al. (2016); Lipparini and Mennucci (2016) (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 fully variational FLO-SIC, as well as 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 PSI4Smith et al. (2020) could be provided in the future.

An important feature of PySCF are calculations on crystalline systems through the use of periodic boundary conditions.Sun et al. (2018, 2020) Although calculations on periodic solids with exact exchange are tractable within a GTO basis set as in PySCF,Pisani, Dovesi, and Roetti (1988) 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, e.g., 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.Perdew and Zunger (1981); Heaton, Harrison, and Lin (1983); Heaton and Lin (1984); Svane and Gunnarsson (1988a, b); Erwin and Lin (1988); Szotek, Temmerman, and Winter (1990a, b); Svane and Gunnarsson (1990a, b); Svane (1992); Rieger and Vogl (1995); Vogel, Krüger, and Pollmann (1995); Arai and Fujiwara (1995); Vogel, Krüger, and Pollmann (1996); Svane (1996); Svane et al. (2000); Filippetti and Spaldin (2003); Bylaska, Tsemekhman, and Jónsson (2004); Lüders et al. (2005); Bylaska, Tsemekhman, and Gao (2006); Pemmaraju et al. (2007); Hourahine et al. (2007); Stengel and Spaldin (2008); Däne et al. (2009); Gudmundsdóttir, Jónsson, and Jónsson (2015); Nguyen et al. (2018) 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 as well as 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–BoysMarzari and Vanderbilt (1997) or generalized Pipek–Mezey Wannier functions.Jónsson et al. (2017) 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 Su and Goddard (2007) or lithium) and metal-organic frameworks (MOFs).Trepte, Schwalbe, and Seifert (2015); Schwalbe et al. (2016); Trepte et al. (2017, 2018); Trepte and Schwalbe (2019)

Data availability

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.

Acknowledgements.

J. Kortus and S. Schwalbe have been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)

• Project-ID 421663657 - KO 1924/9-1. J. Kortus and J. Kraus have been funded by the DFG - Project-ID 169148856 - SFB 920, subproject A04. K. Trepte has been 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 Number #DE-SC0018331. S. Lehtola has been supported by the Academy of Finland (Suomen Akatemia) through project number 311149. We also thank the ZIH in Dresden for computational time and support, Prof. 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.