Effect of molecular and electronic geometries on the electronic density in FLO-SIC

Recently, Trepte et al. [J. Chem. Phys., vol. 155, 2021] pointed out the importance of analyzing dipole moments in the Fermi-L\"owdin orbital (FLO) self-interaction correction (SIC) for cyclic, planar molecules. In this manuscript, the effect of the molecular and electronic geometries on dipole moments and polarizabilities is discussed for non-cyclic molecules. Computed values are presented for water, formaldehyde, and nitromethane. Continuing the work of Schwalbe et al. [J. Chem. Phys. vol. 153, (2020)], we reconfirm that systematic numerical parameter studies are essential to obtain consistent results in density functional theory (DFT) and SIC. In agreement with Trepte et al. [J. Chem. Phys., vol. 155, 2021], DFT agrees well with experiment for dipole moments, while SIC slightly overestimates them. A Linnett double-quartet electronic geometry is found to be energetically preferred for nitromethane.

Zunger self-interaction correction (PZ-SIC) [23] approximately removes the oneelectron SIE. It has a long history of successes and failures [24]. In PZ-SIC, the choice of orbitals is important. Lehtola et al. [25] showed that PZ-SIC suffers from the local minima problem. A recent formulation of PZ-SIC utilizes so-called Fermi-Löwdin orbitals (FLO-SIC) [26][27][28][29][30]. FLO-SIC depends on Fermi-orbital descriptors (FODs) [29] to construct the localized orbitals used for PZ-SIC [31]. These FODs can be imagined as semi-classical electron positions which form an electronic geometry. Recently, Trepte et al. [31] showed that one can guide and classify local minima in PZ-SIC with the help of special sets of FODs that reflect chemical bonding theories of Lewis [32] and Linnett [33,34]. The latter is known as Linnett's double-quartet (LDQ) theory. While typically one is interested in the variational total energy of the system, Trepte et al. [31] proposed to additionally monitor the dipole moment to classify PZ-SIC solutions. The dipole moment is one of the most simple descriptor for the electronic density -the key property in any DFA.
In this work we investigate the influence of molecular and electronic geometries as well as a properly chosen parameter space on the quality of densityrelated properties in DFT and FLO-SIC. We show that numerical parameters need to be optimized not only for the total energy but also for, e.g., the electric dipole moments and/or polarizabilities. We discuss the results based on small, illustrative and educationally-valuable molecules.
The manuscript is structured as follows. In the first two sections we outline the theoretical background and the computational details, after which we present the major results. In the last section we summarize and conclude our findings.
where 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 functional, K XC [n α , n β ] is the exchange-correlation (XC) functional, n is the electron density, and α and β indicate spin channels.
To compute the XC functional one needs to evaluate an explicit density integral using a numerical quadrature, see Section 3. Here, ε hom XC [n α , n β ] is the XC energy-density of the homogeneous electron gas and F XC is an XC enhancement factor.
Several approximations exist for the XC enhancement factor, many of which are available in Libxc [35]. These approximations lead to artificial interactions of electrons with themselves; this is called self-interaction (SI). The corresponding SI energy comes from an incomplete cancellation of the exchange-correlation energy and the Coloumb energy for one-electron densities n σ In PZ-SIC, the total E KS is corrected orbital-by-orbital as A novel flavor of PZ-SIC is FLO-SIC. This formulation utilizes FODs to construct Fermi orbitals (FO). These FOs are then orthogonalized to become FLOs. The FODs need to be optimized in the employed numerical parameter space using the respective analytical gradients [36]. With energy and gradient expressions at hand one can study, e.g., ionization potentials, atomization energies or barrier heights. However, guided by Trepte et al. [31] our focus is not only on energies, i.e., Eq. (1) and Eq. (4), but on properties characterizing the density. Thus, having introduced energy expressions for DFT and PZ-SIC, we continue to discuss dipole moments and polarizabilities as fingerprints of the electron density.
Density-related properties can be analyzed using small applied electric fields. The total energy of a system under an external electrical field ε can be written as Here, E 0 refers to a ground state energy, e.g, KS-DFT E KS (see Eq. (1)) or PZ-SIC E PZ (see Eq. (4)). From this energy expression we can derive the electric dipole moment as Commonly, the dipole moment is directly calculated from the electronic density where Z A , R A , and n(r) are nuclear charges and positions and the total electron density, respectively. Note, we only discuss the electric dipole moment in this work and refer to it simply as dipole moment. The dipole moment is a measure for the polarity of a system and tells us about the charge separation in this system. As an example, the molecular electrostatic potential [39] (MEP) and Fig. 2. Illustration of the dipole moment µ in the H2O molecule using the geometry [37] from CCCBDB [38]. The molecular electrostatic potential (MEP) [39] was calculated using DFT with LDA-PW, aug-pc-3, and a grid= (200,1454) in PySCF. The visualization was done using Jmol [40].
the dipole moment of the H 2 O molecule is visualized in Fig. 2. The dipole points from the more electronegative O atom to the less electronegative H atoms. This can be directly seen in the coloring scheme of the MEP.
Having the possibility to calculate the total energy of the system under the influence of an external electric field allows to study dipole moments using, e.g., a 2-point finite difference (FD) stencil One can also derive the electric polarizability α ij from Eq. (5) as The electric polarizability can be calculated using analytical approaches, e.g., solving the coupled perturbed Hartree-Fock (CPHF) equation [41][42][43]. It describes the tendency of a system to acquire an induced dipole moment in the presence of an external electric field. Similar to the dipole moment we will refer to the electric polarizability as polarizability.
With the dipole moment at hand, i.e., Eq. (7) or Eq. (8), one can calculate the directional components of the polarizability tensor as vector components which is a row in Eq. (11).
To further simplify the characterization of the density, we introduce scalar values. The vectorial dipole moment µ = (µ x , µ y , µ z ) will be represented as µ = |µ|, while the tensorial polarizability will be represented as α = Tr(α ij )/3. As shown by Trepte et al. [31], dipole moments are sensitive to the chemical bonding described by the electronic geometry. Moving FODs changes the dipole moment in FLO-SIC. Dipole moments in FLO-SIC provide insights into whether the electronic density respects the molecular symmetry or not [31].
In practical FLO-SIC calculations, the quality of the density is determined by the used basis set, numerical quadrature, molecular and electronic geometry. Computational details and the used electronic structure codes are discussed next.

Computational details
All scripts to produce the data presented in the manuscript are available at https: //gitlab.com/opensic/dippo [44]. The calculations were performed with the allelectron Gaussian-type orbital (GTO) codes PySCF [42] and PyFLOSIC2.
PyFLOSIC2, see https://gitlab.com/opensic/pyflosic2, is the successor of PyFLOSIC [30]. It offers a cleaner and more modular code structure and can now easily be installed via the Python package manager pip. For calculations of real-and complex-valued SIC, i.e., RSIC and CSIC [25,45,46], we used the ERKALE code [47]. In previous studies [30,31] we observed that the pc-n basis sets [48,49] perform well for DFT as well as SIC calculations. Therefore, for all calculations in this work we use pc-n basis set variants. All codes use the Libxc [35] library, offering access to a vast variety of exchange-correlation functionals. From this library we access LDA-PW [50], PBEsol [51], and r 2 SCAN [22]. The used codes are Open-Source codes, meaning they are freely available to anyone [52,53]. Open-Source codes enable faster code development, re-usable concepts, and versatile tool-boxes.
PySCF and PyFLOSIC2 are written in Python, where only numerically demanding parts in PySCF are written in C. Python is simple and elegant, has a friendly and helpful community, and provides various well-maintained libraries. These are only some reasons why it is easy for students or non-programmers to start coding with Python. This allows to solve even non-trivial tasks, like writing a DFT code from scratch [54] in the limited time of a master thesis when guided and educated with novel strategies [55].
A numerical quadrature [56] is needed to evaluate XC properties in DFT and SIC, see Eq. (2) in Section 2. We will refer to it simply as grid. A typical grid consists of a radial and an angular part. Its size is given as a pair of numbers, i.e., the number of radial shells and the number of angular points. SIC requires significantly finer grids than DFT [45,46,57]. In analogy to [30,31], we prune the used grids neither for DFT nor for SIC. This is done as the orbital densities evaluated in SIC are not as smooth as the total density used in DFT [58], thus requiring a finer resolution [46].
FLO-SIC has two major variational degrees of freedom, the density matrix (DM) and the FODs. All FLO-SIC calculations in this work are realized with a two-step FLO-SIC SCF cycle, which follows the idea proposed by [59]. In FLO-SIC, the initial DM and initial molecular coefficients are typically the ones from a DFT calculation. The initial FODs can be generated with various procedures, e.g., Python-based center of mass PyCOM, Fermi-orbital descriptor Monte-Carlo fodMC, or other so-called FOD generators [29]. All initial FOD configurations used in this work were generated using the fodMC.
Within a full FLO-SIC calculation, the FODs are fully optimized in an inner FOD loop for a given DM. All SIC properties are calculated for this DM and the respective optimized FODs. Based on this, the unified Hamiltonian [30,[60][61][62] is updated which then provides the next DM in the outer DM loop. The initial FODs in the inner FOD loop are optimized using the SciPy L-BFGS-B [63][64][65][66][67][68][69][70] optimization algorithm with a maximum force component threshold of f max,tol = 2 · 10 −4 E h /a 0 . This two-step procedure is repeated until the FOD forces reach f max,tol and the DM is not changing anymore.
The computational methods introduced in this section have been applied to calculate dipole moments and polarizabilities. The results for water, formaldehyde, and nitromethane are discussed in the next section.

Sisyphus rock: The importance of grid and basis set size
For the correct description of density-related properties, it is important to converge the used numerical parameters space consisting of grid and basis set. Systematic parameters studies have been performed utilizing CCCBDB molecular geometries [38], and FODs generated with the fodMC in the case of FLO-SIC to exemplify this. The determined trends and optimal values should be transferable to other molecular geometries and other FOD arrangements.
We carried out systematic grid convergence tests varying the number of radial shells with a fixed number of angular points and vice versa. Detailed information can be found at https://gitlab.com/opensic/dippo [44]. The DFT values converge at smaller grids than the respective SIC values, see Table 1. For LDA-PW, a value of N rad = 200 for the radial shells gives converged results for both DFT and SIC. The angular dependency for DFT as well as for FLO-SIC is converged at N ang = 590. However, we use a value of N ang = 1454 to resolve all one-electron and total densities accurately. Table 1. Optimal numerical parameters for DFT and FLO-SIC regarding the molecular geometries reported in CCCBDB [38] for H2O [37], CH2O [71] and CH3NO2 [72].

System
DFT FLO-SIC n rad nang basis set n rad nang basis set H2O 100 302 aug-pc-3 150 590 aug-pc-3 CH2O 100 302 aug-pc-3 200 590 aug-pc-3 CH3NO2 150 590 aug-pc-3 200 590 aug-pc-3 Note that we investigated the convergence of the density in DFT using the dipole moment as well as the polarizability. Both density fingerprints deliver the same optimal parameters for the tested molecules. Thus, for FLO-SIC we only used the dipole moment to determine optimal parameters.
Having an optimal grid of (200,1454) enabled us to determine a suitable basis set. The convergence of the basis set for water is shown in Fig. 3 using DFT and in Fig. 4 using FLO-SIC. The aug-pc-3 basis set shows convergence w.r.t. the basis set size. Those parameters are optimal for water, formaldehyde, and nitromethane. However, such convergence checks need to be done for any molecule -a Sisyphus work. Otherwise, the meaning of absolute values for dipole moments or polarizabilities are questionable. Convergence of the total energy Etot, the dipole moment µ and the polarizability α w.r.t. increasing basis set size for LDA-PW DFT using PySCF. We used pc-n, aug-pc-n, and unc-aug-pc-n with n=0-4 [48,49]. Each plot shows the difference to the largest used basis set. We used pc-n, aug-pc-n, and unc-aug-pc-n with n=0-4 [48,49]. Each plot shows the difference to the largest used basis set. Only the density matrix was optimized, while the FODs were fixed.
Having established a suitable numerical parameter space, we continue to discuss the influence of molecular geometries on dipole moments and polarizabilities.

Pandora's box: The quality of molecular geometries matters
A molecular geometry is needed to perform electronic structure theory calculations; in case of DFT see Fig. 1. While such geometries can be optimized within most theories, it is not uncommon to use a fixed molecular geometry to be comparable to other approximations or to simply save computational time.
For quantum chemical calculations there exist a vast a variety of seemingly promising databases, such as CCCBDB [38], ChemSpider [73,74], PubChem [75], and many more. However, the quality of the geometries in these databases can vary [76]. CCCBDB offers access to a vast variety of molecular geometries. In this work, we used the CCCBDB geometries which can be found in the experimental section. Note that those geometries are not necessarily experimental ones. For example, for water [37] the geometry is derived semi-empirically utilizing experimental reference values. PubChem provides molecular geometries calculated using the MMFF94s [77] force field. For other databases such as ChemSpider it is even not that trivial to find the quality of the geometries.
The question of the quality of the molecular geometry might be important for other fields as well. For example, machine learning [78] models might be trained on low quality geometries, which could affect the predictability of the resulting models. For SIC calculations, the quality of the molecular geometry is of great importance, as the orbital densities are sensitive to the underlying molecular geometry. Molecular geometry optimizations are a standard task for commonly used approaches like DFT. However, for more computational demanding methods like FLO-SIC, full geometry optimization require high computational effort. This is caused by the coupled degrees of freedoms of nuclei and FODs [31].
The following results are based on the small, educational systems H 2 O, CH 2 O, and CH 3 NO 2 . We summarized some essential molecular information for those molecules in Table 2, including experimental reference values for dipole moments and polarizabilities.   Table 2.
The effect of molecular geometries on density fingerprints, i.e., dipole moments and polarizabilities, is significant even at the DFT level (see Fig. 5). For the used test systems the molecular geometries from CCCBDB provide the best molecular geometries; the dipole moments and polarizabilities are close to the values obtain from a DFT optimized geometry. The basis set size affects the density-related properties significantly, and only the aug-pc-3 basis set provides converged results. The molecular geometries from PubChem and ChemSpider should be used with care, as the dipole moments deviate significantly from values obtained from a DFT optimized geometry.
For all calculations, the finite difference approximations, see Eq. (8) and Eq. (10), agree well with the analytical results. The mean error for the dipole moments is 0.00 D, while the respective mean errors for the polarizabilities are given in Table 3. Accordingly, for the investigated systems and the employed method the chosen value of 10 −7 a.u. for the magnitude of ε regarding the 2-point finite difference approximation delivers reliable numerical results. Note that this finding might not be reproducible for other systems using the same value. Having examined the dependence on the molecular geometry, in Section 4.3 we proceed to investigate electronic degrees of freedom in FLO-SIC, i.e., DM and FODs.

The sword of Damocles: Curse and blessing of approximations
Approximations are often needed and can be useful to enable the treatment or computation of a specific property at a certain level of theory. However, each approximation needs to be carefully investigated regarding the limits of its predictive power.
Effect of initial FODs Continuing the work of Trepte et al. [31], we show for CH 3 NO 2 that it is possible to find several FOD configurations following chemical bonding theories. One possible FOD configuration can be based on Lewis theory of bonding. There, the FODs form one double N=O and one single N-O bond in the -NO 2 group, see  In electronic structure theories there exist several possibilities to treat the spin of the system. In FLO-SIC, we can do restricted calculations where all electrons are paired, N α = N β , and unrestricted calculations where N α and N β can vary. In restricted FLO-SIC only one set of FODs is needed, while in unrestricted FLO-SIC two sets of FODs are required.
Given that the calculations for CH 3 NO 2 are more computational demanding, we use a grid=(200,1454) and the aug-pc-2 basis set. As seen in Section 4.1, using this basis set comes with errors in the order of mE h and mD w.r.t. to the basis set limit. However, the energy difference of the considered FOD configurations, see Fig. 6, are in the order of 10 −2 E h . The differences in the dipole are in the order of 10 −1 D. Thus, using the aug-pc-2 basis set should deliver reliable trends.
We performed restricted FLO-SIC calculations for the Lewis configuration of CH 3 NO 2 . The FODs converge to an electronic geometry which does not follow any bonding theory. The double bond FODs are not lying on the N-O bond axis. Instead, they are distorted towards the respective O atoms. Given this non-symmetric arrangement of FODs, the density of CH 3 NO 2 becomes nonsymmetric. This leads to an energy of E final = −245.80939 E h and a nonsymmetric dipole moment with an absolute value of µ = 4.06 D.
The effect of various unrestricted FOD configurations for CH 3 NO 2 is shown in Fig. 6. The differences in these configurations can already be seen for the initial FODs. Only optimizing the density shows significantly different dipole moments, and only the LDQ value is close to the experimental value. The energy for the initial LDQ arrangement is also the lowest. Upon full optimization of the density matrix as well as the FODs, the configuration based on LDQ deliver the lowest energy. In the unrestricted calculations, both the initial LDQ and the initial Lewis configuration converge to a final FOD arrangement which can be characterized via LDQ. Interestingly, the other FODs stay in the their configuration, but the final energy is higher than for the LDQ FODs.
LDQ not only gives the lowest energy but also the best SIC dipole moment for CH 3 NO 2 . As a highlight, the center of mass (COM) for the optimal RSIC localized orbitals using ERKALE also reflect the LDQ chemical bonding motif, see the RSIC box in Fig. 6. In this section, we reconfirm that the choice of the initial FODs significantly influences dipole moments in FLO-SIC. Next, we describe how the FOD optimization itself can influence dipole moments.

Effect of FOD optimization
We continue the discussion for H 2 O and CH 2 O. These two molecules only have one meaningful FOD configuration, see Fig. 7. Having one FOD configuration may allow finding the same local minima in a reproducible fashion. Thus, those molecules are promising candidates for systematic FOD convergence studies.
In FLO-SIC it is rather easy to make various approximations. A common approximation is to use fixed FODs which are not optimized for the numerical parameter space of the respective calculation. Optimizing the FODs is, however, critical to obtain reasonable FLO-SIC solutions. Such optimizations are carried out until a specific threshold for the maximum FOD force, f max,tol , is reached. The influence of the FOD optimization on the total energy and the dipole moment for water and formaldehyde is shown in Fig. 8. FLO-SIC values for the total energy and the dipole moment change drastically when going from an initial set of FODs to optimized FODs. Here, optimized FODs are characterized by f max,tol = 2 · 10 −4 E h /a 0 . Stopping the optimization too early, i.e., at f max,tol = 5·10 −3 E h /a 0 can lead to insufficiently converged densities and energies. This furthermore leads to an incorrect prediction of trends; increasing the basis set size should smoothly converge the dipole moments, see Fig. 5. However, with insufficiently optimized FODs these trends can be predicted incorrectly, as can be seen in Fig. 8.
Accordingly, global statements and generalization of trends are only valid for optimized FODs in combination with a sufficient basis set and grid. In the tested cases, the aug-pc-3 basis set with a grid=(200,1454) deliver converged results, in analogy to Sections 4.1 and 4.2.
As proposed by the authors [31], monitoring the dipole moment is important to classify and analyze SIC solutions. We showed here that converging the total energy is necessary, but might not be sufficient when one aims to study densityrelated properties, i.e., dipole moments or polarizabilities. All calculations in the previous sections utilized the LDA-PW functional. Next, we discuss how changing the exchange-correlation functional influences dipole moments.
Effect of exchange-correlation functional For H 2 O and CH 2 O, we find that already LDA-PW describes dipole moments qualitatively correctly for DFT and SIC. However, comparing three exchange-correlation functionals the absolute dipole moments differ, see Table 4. Note that for CH 3 NO 2 this comparison was not carried, see Effect of initial FODs for more information. All three pure DFAs agree quite well with the experimental dipole reference values, see Table 2. FLO-SIC tends to overshot the dipole moment in order of 10 −1 D w.r.t. the DFT values. This trend has also been observed in the literature, see [31]. For DFT, LDA-PW performs best for H 2 O while r 2 SCAN agrees the most for CH 2 O. In case of FLO-SIC, the best dipoles are given by PBEsol for H 2 O and by r 2 SCAN for CH 2 O.
In the previous sections we have discussed results for FLO-SIC. To investigate the influence of a specific flavor of SIC on dipole moments, we compare FLO-SIC, RSIC, and CSIC in the next section.
Effect of SIC methods To verify the used numerical parameter space, i.e., the aug-pc-3 basis set and a grid=(200,1454), we carried out RSIC and CSIC calculations using ERKALE. There is no significant difference between FLO-SIC and RSIC values, see Table 5. This is noteworthy as the results are calculated with two independent electronic structure codes. Thus, the used numerical parameter space is sufficient to deliver reproducible results. SIC suffers from the multiple local-minima problem [25], see Section 1. Accordingly, we recommend to verify FLO-SIC results with independent SIC methods such as RSIC. For example for nitromethane, see Effect of initial FODs, only the LDQ FLO-SIC dipole agrees with the RSIC solution. FLO-SIC, RSIC and CSIC together deliver a consistent SIC description. The advantage of FLO-SIC is the access to bonding information. This information allows to easily classify and further analyze PZ-SIC solutions [31]. In the next section we summarize and conclude our findings.

Summary and Conclusion
As shown in Section 4.1, and already stated in earlier works [30,31,45], SIC needs finer numerical quadrature meshes in comparison to DFT calculations. Global statements about the predictive power of SIC are only meaningful using very accurate numerical parameter spaces. In this work we show how density-related properties, i.e., dipole moments and polarizabilities, can help to determine an appropriate numerical parameter space for DFT and SIC. For our investigated molecules, the aug-pc-3 basis set and a numerical grid of (200,1454) is such an appropriate parameter space. Note, even trends can clearly change using other numerical parameters. While it is mandatory to converge the energy it is not necessarily sufficient for the study of density-related properties.
Furthermore, using water, formaldehyde, and nitromethane we show that these density-related properties are not only sensitive to the used numerical parameter space. They are also significantly influenced by the used molecular geometry, see Section 4.2. Molecular geometries from common chemical database, i.e., CCCBDB, PubChem, or ChemSpider, can deliver very different dipole moments and polarizabilities. Only molecular geometries from the CCCBDB database deliver reasonable trends in comparison to optimized geometries.
Continuing the work of [31], for FLO-SIC we showed that density-fingerprints are also sensitive to the chemical bonding situation introduced by FODs. We demonstrate that the numerical quality of the FODs, represented by their gradients, as well as the choice of initial FODs clearly influences the dipole moments. For molecules with non-trivial bonding situations, e.g., nitromethane, it is highly recommended to use various initial FOD configurations. This allows to reasonably sample the FOD configuration space to determine the most reasonable FOD configuration. Such FOD configurations typically follow chemical bonding theories [31], i.e., Lewis and LDQ. In case where several FOD configurations are possible, FODs based on LDQ are often superior in FLO-SIC. This has been shown in [31] and verified here again for nitromethane in Section 4.3.
When computational time matters, the simplest density functional approximation, i.e., a local density approximation like LDA-PW, provides reasonable trends for dipole moments in FLO-SIC, as seen in Section 4.3. If computational time is not limited we recommend to use higher rung functional like PBEsol or r 2 SCAN to verify and further analyze determined trends.
We want to emphasis the importance of Open-Science and Open-Source developments, as this work would not be possible without them.

Acknowledgement
S. Schwalbe has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Project ID 421663657 -KO 1924/9-2. S. Liebing wants to express his gratitude for Prof. Jens Kortus enabling the authors to work on this topic. All authors thank M. Sc. Wanja T. Schulze for his contributions to the PyFLOSIC2 code. We thank both referees, as their comments significantly improved the content and shape of our manuscript. This work is part of our OpenSIC project, and we thank all members for fruitful discussions. We thank the ZIH Dresden for computational time and support.