Full self-consistency in Fermi-orbital self-interaction correction

TL;DR

This paper introduces a new algorithm for FLOSIC that achieves full self-consistency, improving accuracy and lowering total energies in density functional calculations.

Contribution

A novel algorithm for FLOSIC that ensures full self-consistency with minimal additional computational cost.

Findings

Total energies are lower with the new self-consistent method.

The new algorithm maintains computational efficiency.

Improved accuracy in systems with strongly changed occupied manifolds.

Abstract

The Perdew-Zunger self-interaction correction cures many common problems associated with semilocal density functionals, but suffers from a size-extensivity problem when Kohn-Sham orbitals are used in the correction. Fermi-L\"{o}wdin-orbital self-interaction correction (FLOSIC) solves the size-extensivity problem, allowing its use in periodic systems and resulting in better accuracy in finite systems. Although the previously published FLOSIC algorithm [J. Chem. Phys. 140, 121103 (2014)] appears to work well in many cases, it is not fully self-consistent. This would be particularly problematic for systems where the occupied manifold is strongly changed by the correction. In this paper we demonstrate a new algorithm for FLOSIC to achieve full self-consistency with only marginal increase of computational cost. The resulting total energies are found to be lower than previously reported non-self-consistent results.

Tables

Table 1. Table 1: Total energies and atomization energies of 12 molecules calculated with SC-FLOSIC and non-SC-FLOSIC. We use the experimental geometries in all calculations. The experimental atomization energies are obtained from Ref. CCC . Zero-point vibrational energies have been removed from the experimental values. The mean absolute error (MAE) and the mean absolute relative error (MARE) for the atomization energies are also listed.

Mol.	Total (Hartree)			Atomization (eV)
	LSDA	non-SC	SC	LSDA	non-SC	SC	Exp.
N2	-108.692	-109.842	-109.856	11.53	10.22	10.19	9.76
O2	-149.332	-150.736	-150.757	7.48	5.10	5.14	5.12
CO	-112.471	-113.635	-113.648	12.91	11.18	11.10	11.11
CO2	-187.273	-189.096	-189.121	20.37	16.27	16.24	16.56
C2H2	-76.625	-77.594	-77.607	19.91	18.97	18.92	16.86
LiF	-106.702	-107.716	-107.734	6.75	5.64	5.72	5.97
H2	-1.1251	-1.1745	-1.1758	4.57	4.81	4.79	4.48
Li2	-14.724	-15.050	-15.053	1.02	1.02	0.96	1.04
CH4	-40.109	-40.667	-40.676	19.77	19.64	19.62	17.02
NH3	-56.107	-56.761	-56.772	14.60	14.45	14.45	12.00
H2O	-75.909	-76.665	-76.677	11.53	10.69	10.69	9.51
CH3OH	-114.84	-116.13	-116.15	25.30	24.85	24.83	20.85
MAE				2.13	1.16	1.14
MARE				18%	8.5%	8.7%

Equations

$$E_{\mathrm{XC}}^{\text{SIC}}[n_{\uparrow },n_{\downarrow }]=E_{\mathrm{XC}}[n_{\uparrow },n_{\downarrow }]-\sum _{\sigma }\sum _i^{N_{\sigma }}\left\{U[n_{i\sigma }]+E_{\mathrm{XC}}[n_{i\sigma },0]\right\},$$

$$⟨{\varphi }_{j\sigma }\left|V_{j\sigma }^{\text{SIC}}-V_{i\sigma }^{\text{SIC}}\right|{\varphi }_{i\sigma }⟩=0$$

$$\left({\hat{H}}_{\sigma }^{\text{KS}}+V_{i\sigma }^{\text{SIC}}\right){\varphi }_{i\sigma }=\sum _j^{N_{\sigma }}{\lambda }_{ji\sigma }{\varphi }_{j\sigma },$$

$${\varphi }_{i\sigma }^{\text{FO}}(\mathbf{r})=\frac{n_{\sigma }({\mathbf{a}}_{i\sigma },\mathbf{r})}{\sqrt{n_{\sigma }({\mathbf{a}}_{i\sigma })}},$$

$${\hat{H}}_{\sigma }={\hat{H}}_{\sigma }^{\text{KS}}+\sum _{ij}^{N_{\sigma }}\frac{V_{ij}^{(i\sigma )}+V_{ij}^{(j\sigma )}}{2}|{\varphi }_{i\sigma }⟩⟨{\varphi }_{j\sigma }|,$$

$$H_{mn}^{(i\sigma )}=⟨{\varphi }_{m\sigma }\left|H_{\sigma }^{\text{KS}}+V_{i\sigma }^{\text{SIC}}\right|{\varphi }_{n\sigma }⟩.$$

$${\tilde{H}}_{mn\sigma }=⟨{\varphi }_{m\sigma }\left|H_{\sigma }^{\text{KS}}+V_{m\sigma }^{\text{SIC}}\right|{\varphi }_{n\sigma }⟩,$$

$${\tilde{H}}_{i\le N_{\sigma },j>N_{\sigma },\sigma }$$

$${\tilde{H}}_{i>N_{\sigma },j>N_{\sigma },\sigma }$$

$${\hat{\varphi }}_{i\sigma }=\sum _j^{N_{\sigma }}{C}_{ij}{\varphi }_{j\sigma },i\le N_{\sigma }.$$

$$\begin{array}{rl}{\varphi }_{i\sigma }^{\prime }& =\cos t{\varphi }_{i\sigma }+\sin t{\varphi }_{j\sigma },\\ {\varphi }_{j\sigma }^{\prime }& =-\sin t{\varphi }_{i\sigma }+\cos t{\varphi }_{j\sigma }.\end{array}$$

$$t=\frac{1}{2}\arctan \left(\frac{2{\tilde{H}}_{ij\sigma }}{{\tilde{H}}_{ii\sigma }-⟨{\varphi }_{j\sigma }\left|H_{\sigma }^{\text{KS}}+V_{i\sigma }^{\text{SIC}}\right|{\varphi }_{j\sigma }⟩}\right).$$

$$\tilde{\phi}_{i\sigma}=\left\{\begin{array}[]{ll}\sum_{j=1}^{N_{\sigma}}\left<\phi_{j\sigma}\right.\left|\tilde{\psi}_{i\sigma}\right>\left|\phi_{j\sigma}\right>&i\leq N_{\sigma},\\ \sum_{j=N_{\sigma}+1}^{M}\left<\phi_{j\sigma}\right.\left|\tilde{\psi}_{i\sigma}\right>\left|\phi_{j\sigma}\right>&i>N_{\sigma}.\end{array}\right.$$

Full text

I Introduction

The Kohn-ShamKohn and Sham (1965) (KS) density-functional theoryKohn and Sham (1965); Hohenberg and Kohn (1964); Fiolhais et al. (2003) (DFT) is a widely used computational method in physics and chemistry. To perform a DFT calculation, one needs to approximate the so-called exchange-correlation (xc) energy, a functional of the electronic density. Most of the available approximations are semilocal approximations, where the xc energy density at a certain point only depends on the density and derivatives of the density at that point. These approximations have been extensively used thanks to their computational efficiency and they generally work well, but they all contain the self-interaction error (SIE), which means that the sum of the Hartree interaction energy and the approximated xc energy does not properly vanish for all one-electron systems. SIE causes a wide range of problems, such as the xc potential decaying too fast asymptotically, the orbital energies of occupied orbitals lying too high in a non-systematic way, and the total energy varying in a strongly nonlinear way between adjacent integer electron numbers. These problems show up in DFT calculations as missing Rydberg states, mis-ordering of states for some systems, incorrect description of stretched bonds, problems with band gaps, unstable anions, and similar uncertainties in modeling processes that depend upon such phenomena.

The self-interaction correctionPerdew and Zunger (1981) (SIC), proposed by Perdew and Zunger (PZ), cures these problems by introducing orbital-dependent corrections to the xc energy functional:

[TABLE]

where $\sigma$ is the spin index, $N_{\sigma}$ is the number of occupied orbitals of spin $\sigma$, and $n_{i\sigma}=|\psi_{i\sigma}|^{2}$ is the orbital density of the KS orbital $\psi_{i\sigma}$. Eq. (1) can be directly applied to any xc energy functional, and in principle an approximated density functional could be developed within a systematic self-interaction corrected formalism. The correction vanishes when $E_{\scriptscriptstyle\rm XC}$ is the exact xc functional. Eq. (1) corrects the SIE, but it can also deviate from the KS scheme of DFT, which only allows a local multiplicative potential for the system. While through the optimized effective potential (OEP)Sharp and Horton (1953); Talman and Shadwick (1976) technique the SIC calculation can be done within the KS schemeKörzdörfer et al. (2008), the cost of doing so can outweigh the benefit. Also, a recent studyYang et al. (2016) supports the conclusion that properties such as the band gap are inherently inaccurate in the KS scheme, so there is less reason to use the OEP procedure in practical calculations. The commonly used SIC functional is thus defined in the generalized KS schemeSeidl et al. (1996). The computational cost of PZSIC is higher than that of regular DFT, since the problem has yet to be reduced to a simple eigenvalue problem.

PZSIC works well for small molecules, but there is a flaw that the correction for extended systems can vanish in the bulk limitPerdew and Zunger (1981) if the SIC energy is evaluated naively with KS orbitals. The reason is easy to see. For periodic systems, the KS orbitals all contain a normalization constant of $1/\sqrt{V_{\text{crystal}}}$. When using KS orbitals in Eq. (1), the magnitude of the correction per electron decreases as the size of the crystal increases, and eventually vanishes in the limit of an infinitely large crystal. The problem also shows up in stretched bonds and other situationsKlüpfel et al. (2011). This is the so-called size-extensivity problem, that the energy of separated fragments does not equal the sum of energies of the fragments as it should. For this reason, Perdew and Zunger suggested that the SIC functional could be constructed in terms of a set of possibly-localized orbitals that minimize the SIC energy. However, this only solves the problem when each SIC term in the energy is negativePerdew and Zunger (1981) after localization. This is a problem with the original PZSIC formulation.

The total energy of a regular DFT calculation is invariant with respect to unitary rotations of KS orbitals, but this is no longer the case for SIC calculations, since the SIC term depends explicitly on orbitals. By applying a unitary rotation that transforms occupied KS orbitals to localized orbitals, one can often lower the energy and achieve a finite correction in periodic systems. Localization schemes have been explored thoroughly by the Lin groupHeaton et al. (1983); Pederson et al. (1984, 1985); Pederson and Lin (1988), where the localized orbitals that minimize the energy are obtained by solving the following localization equationPederson et al. (1984, 1985)

[TABLE]

together with the KS-SIC equation of localized orbitals:

[TABLE]

where $\phi$ denotes a localized orbital, $\sigma$ is the spin index, $V^{\text{SIC}}_{i\sigma}({\bf r})=-\delta\left\{U[n]+E_{\scriptscriptstyle\rm XC}[n,0]\right\}/\delta n({\bf r})|_{n({\bf r})=\left|\phi_{i\sigma}({\bf r})\right|^{2}}$ for occupied orbitals, and $V^{\text{SIC}}_{i\sigma}({\bf r})=0$ for virtual orbitals. We denote the canonical orbitals as $\psi$, which diagonalize the Lagrange multiplier matrix $\lambda_{\sigma}$ (the diagonal elements are orbital energies in SIC) and are the analog of KS orbitals in SIC. The sets of $\psi$ and $\phi$ are related to each other by a unitary transformation. Eq. (3) can also be formulated in terms of canonical orbitals, but it is computationally preferable to solve for the localized orbitals directly instead of calculating the canonical orbitals first. The set of localized orbitals that satisfy Eq. (2) and (3) is not unique, and in practice one has to be careful in choosing the initial guess when solving these equations to avoid converging to a local minimum. The localization equation method works generally well, but the computation cost is prohibitively high, as it scales as $O(N^{6})$in the small N limitPederson et al. (2014). As a comparison, the coupled cluster singles doubles (CCSD)Purvis and Bartlett (1982) method has a similar scaling, but it is typically more accurate than DFT methods. Therefore the localization equation method has only been applied to small systems, and it has little advantage compared to methods of similar scaling.

The Fermi-Löwdin-orbital (FLO) SIC Pederson et al. (2014) is a recently proposed localization scheme, which circumvents the need for solving Eq. (2). Instead of searching over all possible unitary transformations to minimize the energy, FLOSIC restricts the transformation to a very specific form, where the localized orbitals are derived directly from the so-called Fermi orbitals:

[TABLE]

where $n_{\sigma}(\mathbf{a}_{i\sigma},{\bf r})=\sum_{j}^{N_{\sigma}}\psi_{j\sigma}^{*}(\mathbf{a}_{i\sigma})\psi_{j\sigma}({\bf r})$ is the single-particle density matrix of the KS system, and $\mathbf{a}_{i\sigma}$ is the Fermi orbital descriptor (FOD). The localized orbitals $\{\phi_{i\sigma}\}$ are the Löwdin-orthogonalizedPederson et al. (1985); Löwdin (1962) Fermi orbitals $\{\phi_{i\sigma}^{\text{FO}}\}$. The unitary transformation between canonical orbitals and localized orbitals is determined with only $N_{\sigma}$ parameters ($\mathbf{a}_{i\sigma}$) instead of $N_{\sigma}^{2}$ parameters (the unitary transformation matrix). The FODs can be seen as quasi-classical electron positions, and they usually can be reliably guessed based on physical and symmetry arguments. Therefore, even though the set of FODs that minimizes the energy can also be non-unique, the situation is much better than in the localization equation method.

Motivated by earlier work Heaton et al. (1983), the FLOSIC calculation in its original paper Ref. Pederson et al. (2014) is realized by adding an extra matrix to the Kohn-Sham Hamiltonian $\hat{H}_{\sigma}^{\text{KS}}$ and diagonalizing the resulting matrix $\hat{H}_{\sigma}$:

[TABLE]

where $V_{ij}^{(i\sigma)}=\left\langle\phi_{i\sigma}\left|V_{i\sigma}^{\text{SIC}}\right|\phi_{j\sigma}\right\rangle$. This physically motivated Hermitian hamiltonian contains most of the effects due to the self-interaction-corrected functional and has been shown to yield canonical-orbital eigenvalues that are not too different from the exact canonical-orbital eigenvalues used in earlier versions of SIC, but it is an approximation Pederson et al. (1985); Pederson and Lin (1988). This algorithm has the advantage of its simplicity for implementation, since the problem remains an eigenvalue problem. However, transforming the orbital-dependent SIC to a regular eigenvalue problem introduces non-self-consistency, which is evident from Eq. (5): the canonical orbitals obtained from solving Eq. (5) are restricted inside the occupied Hilbert space of the previous self-consistent-field (SCF) iteration. Even though the orbitals can still drift away from the initially occupied space since $\hat{H}_{\sigma}^{\text{KS}}$ changes between SCF iterations due to density change, the algorithm is not fully self-consistent (SC).

While non-SC FLOSIC energies, including cohesive energies, are often quite accurate if wave functions from a nearly identical functional are used, SC FLOSIC is needed to make sure that the results are reliable in difficult cases such as molecular magnetsPostnikov et al. (2006), and to have a Hellmann-Feynman theorem which allows for the efficient optimization of molecular and unit-cell geometries.

The non-interacting spin-density matrix $n_{\sigma}({\bf r}^{\prime},{\bf r})$ is the matrix in the position representation of the projection operator onto the space or manifold of the occupied orbitals of spin $\sigma$. Thus it is invariant under unitary transformation of those orbitals. We can calculate the total energy from $n_{\sigma}({\bf r}^{\prime},{\bf r})$, and we can also find the electron spin-density $n_{\sigma}({\bf r})=n_{\sigma}({\bf r},{\bf r})$ from its diagonal. We will exploit these facts here to find the optimized FLOs, for a given set of FODs directly. Before the first iteration, we will have approximate canonical occupied orbitals, which generate a set of initial FLOs, and canonical unoccupied orbitals. Then we will mix the occupied and unoccupied orbitals to lower the SIC energy, with a reconstruction of localized FLOs after each iteration. We will not need to generate canonical occupied orbitals after the iteration begins, but they can be obtained easily if needed.

II Algorithm

In this work we set out to find the SC solution of Eq. (3), where the localized orbitals are constrained to be symmetrically orthogonalized Fermi orbitals which are referred to as Fermi-Löwdin Orbitals (FLO). Both the localized orbitals and the Lagrange multiplier matrices need to be determined. Since the Lagrange multiplier matrices are in principle different for hamiltonians defined by different orbitals, solving Eq. (3) directly would mean determining $N_{\sigma}^{3}$ variables for a given set of FODs. Even though we only need to determine $N_{\sigma}$ parameters in the FLOSIC method instead of $N_{\sigma}^{2}$ parameters in the localization equation method, the computational cost is still very high. Instead of trying to solve Eq. (3) directly, we develop a Jacobi-like method for solving Eq. (3) iteratively. Even though we employ several approximations to reduce the computational cost, the final result obtained from the iterative procedure remains an exact solution of Eq. (3) within the chosen basis set.

In the following, we represent all matrices in the basis set formed by the FLOs (indexed $1$ to $N_{\sigma}$) and the virtual KS orbitals (indexed $N_{\sigma}+1$ to $M$, $M$ being the size of the basis set) of the previous SCF iteration, all of which are denoted as $\phi_{i\sigma}$. Note that all $\phi_{i\sigma}$’s are orthogonal to each other. The FLOs are calculated with pre-determined FODs, which are fixed during the SCF loop (refer to Fig. 1). The Hamiltonian matrix of occupied orbital $i\sigma$ in this representation would be

[TABLE]

To reduce the amount of calculations needed, instead of calculating Eq. (6) exactly, we construct a matrix for the operator $\tilde{H}_{\sigma}$

[TABLE]

so that the $i$-th row of $\tilde{H}_{\sigma}$ is the same as the $i$-th row of $H^{(i\sigma)}$. According to Eq. (3), the upper-left $N_{\sigma}\times N_{\sigma}$ corner of $\tilde{H}_{\sigma}$ is the transpose of the Lagrange multiplier matrix $\lambda_{\sigma}$.

We have the following relations if the group of $\phi_{i\sigma}$ solves Eq. (3):

[TABLE]

Eq. (8) holds since the Hilbert spaces spanned by the occupied and virtual orbitals are orthogonal, and Eq. (9) holds since $V^{\text{SIC}}_{i\sigma}=0$ for virtual orbitals.

Eq. (8) only holds when the group of $\phi_{i\sigma}$ solves Eq. (3). Consider the case where trial functions $\hat{\phi}_{i\leq N_{\sigma},\sigma}$ are unitarily-transformed $\phi_{i\sigma}$:

[TABLE]

Unless $C_{ij}=\delta_{ij}$, applying $\hat{H}_{\sigma}^{\text{KS}}+V_{i\sigma}^{\text{SIC}}$ on $\hat{\phi}_{i\sigma}$ will generally yield $\sum_{j}^{N_{\sigma}}\hat{\lambda}_{ji\sigma}\phi_{j\sigma}+\Delta_{i}$, where the residue $\Delta_{i}$ belongs to the virtual Hilbert space. Therefore instead of trying to solve Eq. (3) directly, we can solve Eq. (3) by requiring that Eqs. (8) and (9) be satisfied. This is done by finding a set of $\phi^{\prime}_{i\sigma}$ that makes the upper-right $(M-N_{\sigma})\times N_{\sigma}$ corner of $\tilde{H}$ vanish in the matrix representation where the group of $\phi^{\prime}_{i\sigma}$ forms the basis set. We achieve this by pairwise mixing the FLOs with the virtual KS orbitals. For every pair of FLO $\phi_{i\leq N_{\sigma},\sigma}$ and virtual orbital $\phi_{j>N_{\sigma},\sigma}$, we define $\phi^{\prime}_{i\sigma}$ and $\phi^{\prime}_{j\sigma}$ as a rotation of the original pair:

[TABLE]

We then require that $\left\langle\phi^{\prime}_{i\sigma}\left|H_{\sigma}^{\text{KS}}+V^{\text{SIC}}_{i\sigma}\right|\phi^{\prime}_{j\sigma}\right\rangle=0$ and obtain the equation for the rotation angle $t$:

[TABLE]

The angle $t$ being in the correct direction is more important than its actual value. Not using its exact value only impacts the speed of convergence. We also found that at times convergence can be hastened by stepping a shorter distance than determined from Eq. (12). Such rotations are done for all FO-virtual orbital pairs. After each rotation, we replace $\phi_{i\sigma}$ and $\phi_{j\sigma}$ with $\phi^{\prime}_{i\sigma}$ and $\phi^{\prime}_{j\sigma}$, so that the later rotations will take the earlier rotations into account. Convergence is checked after all pairs of orbitals are rotated, and it is reached when the biggest matrix element in the upper-right corner of $\tilde{H}$ is smaller than a pre-determined criterion. The rotation procedure is repeated until convergence is reached. We find that the updated occupied orbitals are still primarily localized to their original regions and the updated virtual orbitals are still primarily delocalized.

This Jacobi-like method for zeroing-out the upper-right corner of $\tilde{H}$ shows strong dependence on the initial set of $\phi_{i\sigma}$. For some starting points more than 1000 iterations are needed for convergence. To accelerate the convergence, we precondition the starting orbitals by diagonalizing the averaged matrix $(\tilde{H}+\tilde{H}^{T})/2$, which yields a unitary transformation between approximate canonical orbitals of the current iteration and the occupied FLO and unoccupied delocalized orbitals (UDO) of the previous iteration. We then transform back to an FLO/UDO representation before performing the Jacobi-like iterative update. This method for preconditioning always converges within 20 iterations for all systems checked. This can be seen as an equivalent to applying a pre-conditioner to $\tilde{H}$. The upper-left corner of the resulting $\tilde{H}$ matrix(the new Lagrange multiplier matrix $\lambda_{\sigma}$) is also found to be almost diagonal, meaning that the group of $\phi$ after rotation is very close to canonical orbitals. As discussed in greater detail below, this preconditioning approach can be, and for the results here has been, used to converge to approximate canonical orbitals or exact FLO depending on whether a back transformation scheme is used after the exact diagonalization. Working within the FLO representation will most likely be more efficient with respect to the use of order-N methods and has definite advantages for applications of SIC to meta-stable donor-acceptor complexes as demonstrated for LiF in the applications section.

Though the SIC equation only applies to occupied orbitals, it is possible to introduce orbital-dependent potentials for virtual orbitalsHarrison et al. (1983); Hunt and Goddard (1969); Schneider and Collins (1984); Baruah and Pederson (2009); Jackson and Lin (1988); Pederson and Klein (1988). This would not change the complexity of the method, but it can affect the convergence speed. Since the relative mixing of the virtual orbitals with the occupied orbitals changes as the process iterates to self-consistency, the energy gap between the occupied and unoccupied orbitals affects the rate of self-consistency. In an early paper, Harrison, Heaton and Lin Harrison et al. (1983) suggested that the unoccupied orbitals should move in a potential that resembles the SIC potential of the lowest unoccupied localized orbital since that potential is similar to the potential that would arise if the particle-hole interaction is included. This method is similar to the improved-virtual-orbital (IVO) method of Schneider and Collins Schneider and Collins (1984) and a recent rendition by Pederson and Baruah Baruah and Pederson (2009). Approaches such as these lead to more meaningful estimates of localized excitation energies (such as excitions) and allow for a pre-screening of electronically interesting geometries. They will likely decrease the number of self-consistent iterations when the gap is small but may increase the number of self-consistent iterations when the gap is big. Early efforts to use IVO-like approaches to simultaneously determine the ground-state and excited state energies in defects include work by Jackson and LinJackson and Lin (1988) and Pederson and Klein Pederson and Klein (1988).

The method can be modified to produce a group of $\phi_{i\sigma}$ that is close to FLOs instead of canonical orbitals, which can be useful for ensuring consistent convergence to the same ground state in stretched bonds (Fig. 5 is an example) and other difficult cases. In the previously mentioned trick, the starting orbitals are the eigenvectors of $(\tilde{H}+\tilde{H}^{T})/2$, and we denote them as $\tilde{\psi}_{i\sigma}$. We construct auxiliary orbitals $\tilde{\phi}_{i\sigma}$ by projection into the occupied and virtual spaces of the previous SCF iteration:

[TABLE]

The new starting orbitals for the Jacobi-like method are obtained by applying the Löwdin symmetric orthogonalizationPederson et al. (1985); Löwdin (1962) separately to $\{\tilde{\phi}_{i\leq N_{\sigma},\sigma}\}$ and $\{\tilde{\phi}_{i>N_{\sigma},\sigma}\}$. By truncating the part of $\tilde{\psi}_{i\sigma}$ that does not lie within the occupied space of the previous iteration when $i\leq N_{\sigma}$, we make the starting orbitals as close to the FLOs of the previous iteration as possible, and we find the results after the rotations are close to the exact FLOs. With this choice of starting orbitals, the diagonal elements of the Lagrange multiplier matrix are not approximations to orbital energies anymore, as the Lagrange multiplier matrix is no longer near-diagonal. However, the total energy does not depend on the choice of the starting orbitals of the Jacobi-like method.

The algorithm described above only comprises a complete calculation when optimal FODs are known. The optimal FODs can be guessed in very simple systems based on symmetry arguments (such as in $\mathrm{H}_{2}$), but for most of the systems they have to be optimized to ensure that the total energy is at the lowest point. The initial FODs need to be guessed before the first SCF step. Since it is found that FODs often represent chemical bonds and lone pairsPederson et al. (2014); Hahn et al. (2015), making the initial guess is usually simple in first- and second-row elements. We calculate the forces on the FODs after the SCF loop for electronic structures as described in Ref. Pederson (2015). These forces are then used as input of the conjugate-gradient method to generate new FODs, which are used in the next SCF step. This optimization of the FODs is repeated until the forces converge.

Symmetry can be used in the optimization of the FODs to simplify the process and to improve the stability of the optimization. Although this is not a general fact, the optimal FODs for all the systems we tested transform according to a sub-group of the molecular point group. One can easily guess the point group of the FODs for simple systems such as the ones in Table 1, but currently there does not exist a completely systematic way of determining the symmetry of FODs for more complex systems before the calculation.

In regular KS calculations, the KS potential is usually mixed with that of the previous iteration for numerical stability. The non-SC FLOSIC implementationPederson et al. (2014) employs the potential mixing as well, but it was not applied to the SIC potentials. We extend the mixing to the SIC part in this work in order to improve the robustness of the algorithm. Instead of potential mixing, we choose to apply the Broyden mixingJohnson (1988) to the $\tilde{H}$ matrix directly for its simplicity. We obtain similar results and number of SCF iterations for both the potential mixing and the Hamiltonian mixing methods.

III Results

All the calculations present in this paper are done with a modified version of NRLMOLPederson and Jackson (1990); Jackson and Pederson (1990); Pederson et al. (2000), which is an electronic structure code based on Gaussian-type basis sets. We use the local spin density approximation (LSDA)Kohn and Sham (1965); Perdew and Wang (1992) for exchange and correlation because of its simplicity, but there is no difficulty in using other energy functionals together with SC-FLOSIC. The forces on the FODs are converged to 0.001 a.u.

To demonstrate the effect of self-consistency, we apply SC-FLOSIC to the 11 small molecules that were calculated using non-SC-FLOSIC in Ref. Pederson et al. (2014) plus CH3OH. The results are listed in Table 1. FODs in Ref. Pederson et al. (2014) are determined manually and not fully optimized, so some of the non-SC-FLOSIC results listed in Table 1 are different from those in Ref. Pederson et al. (2014).

The self-consistency lowers the total energies systematically as expected. The effect of the self-consistency is small for the small molecules in Table 1. SC-FLOSIC does not further improve the atomization energies over those of non-SC-FLOSIC. To obtain more accurate atomization energies, one has to use more sophisticated energy functionals than the LSDA.

Spin-polarized calculations are needed for O2 and all the atoms in the molecules of Table 1. The FODs for spin-polarized calculations are more difficult to guess, and a poor initial choice usually leads to numerical instabilities. Fig. 2(a)(b) shows the the FODs for N2 and O2 molecules. For the closed-shell N2, the FODs and the shapes of the corresponding FLOs in Fig. 3 have clear chemical importances, as they correspond to the 1s electrons, the lone pairs and the triple N-N bond. As pointed out in Ref. Pederson et al. (2014), the FLOs can be seen as linear combinations of molecular orbitals. For example, the banana orbitals corresponding to the triple bond are linear combinations of the $\sigma_{g}$, $\pi_{x}$ and $\pi_{y}$ molecular orbitals. For open-shell systems such as O2, the link between FODs and the chemical bonds is more elusive. The FODs shown in Fig. 2(b) do not appear to have clear chemical importance. Most of the spin-up and spin-down FODs are not paired as they are in the spin-unpolarized cases, and the FLOs of different spin in Fig. 4 have very different shapes from each other.

For atoms of the first three rows of the periodic table, the FODs are connected to the hybridization of the $s$-$p^{n}$ orbitals with n= 1,2, or 3 depending upon $p$ filling.Pederson and Lin (1988) As for the case of most generalized-gradient approximations, FLOSIC as well as the versions of the PZ functional which account for unitary transformations Heaton et al. (1983); Pederson et al. (1984, 1985); Pederson and Lin (1988); Fois et al. (1993); Szotek et al. (1993); Arai and Fujiwara (1995); Garza et al. (2000); Patchkovskii et al. (2001); Vydrov and Scuseria (2004, 2005); Pemmaraju et al. (2007); Körzdörfer et al. (2008); Pemmaraju et al. (2008); Klüpfel et al. (2011); Hofmann et al. (2012a, b) favors integer occupied $p$ valences over perfectly spherical atoms with fractionally occupied $p$ valences. However, especially when starting new calculations, it is necessary to account for the fact that starting points derived from spherical potentials will always provide fractionally occupied starting orbitals in open-shell systems. To address this issue, we use starting configurations with non-spherical symmetry in calculations. As in the case of antiferromagnetic starting potentials, which are commonly used in molecular magnet calculations and molecular bond-breaking calculations, symmetric solutions would still emerge if they are indeed the lowest-energy solution. Fig. 2(c) shows the FODs of the F atom. The three spin-up $p$ orbitals are fully occupied, and the spin-up FODs form a tetrahedron according to the $sp^{3}$ hybridization. Two spin-down $p$ orbitals are occupied, and the spin-down FODs form a triangle according to the $sp^{2}$ hybridization. The shapes of the corresponding localized orbitals also resemble $sp^{3}$ and $sp^{2}$ hybridized orbitals respectively.

For the systems in Table 1, the optimal FODs of the SC calculation and the non-SC calculation have the same general spatial distribution. The FODs corresponding to core electrons and chemical bonds do not change much from the SC optimal positions in the non-SC calculation, but the FODs corresponding to lone pairs for N2 and CO are strongly affected by self-consistency. For the CO molecule, the distance between the non-SC FOD and the SC FOD corresponding to the lone pair on the O atom can be as long as 0.68Å, and the corresponding non-SC and SC FLOs have visible differences.

The geometry optimization with SC FLOSIC is possible, but deriving and implementing the correct Hellmann-Feynman forces on atomic positions has not yet been completed. Nevertheless, we tested geometry optimization using LSDA forces while using the SC FLOSIC total energy to check convergence. We have limited success for very simple molecules such as H2 and Li2, since the LSDA forces coincides with FLOSIC forces in these simple systems. The correct FLOSIC forces are needed for more general cases. We also find that the optimization of FODs can be merged with the geometry optimization to be more efficient. The initial FODs have to be close to the correct positions in such a combined optimization, or it may not converge.

The SC FLOSIC dissociation curve of the LiF molecule is plotted in Fig. 5. Within most approximations of density functional theory, it is generally a challenge to calculate the charge transfer energy between two fragments without introducing a constraint. For example, constrained orthogonality calculations have been introduced by Baruah and PedersonBaruah and Pederson (2009) and by Ziegler and collaboratorsZiegler et al. (2009) which allow one to calculate an excited charge-transfer state, provided that the ground-state of the system within the approximation to DFT is qualitatively correct. However, in many cases, systems composed of two fragments are found to be electronically unstable for either two neutral separated fragments or two oppositely charged cation/anion fragment pairsSlater (1974); Bally and Sastry (1997); Chermette et al. (2001); Gräfenstein et al. (2004); Ruzsinszky et al. (2006); Constantin et al. (2010), and rare cases where this is not a qualitative problem have been demonstrated in Ref. Ossowski et al. (2003). The LiF molecular potential energy curve provides a good example of this problemOssowski et al. (2003); Vydrov et al. (2007); Makmal et al. (2011). In this case, LSDA total energies as functions of charge on the separated systems predict that the ground state of a separate Li and separated F has approximately 2.3 units of charge on the Li and 9.7 units of charge on the F (this result depends slightly on choice of functional). Here we address the calculation of the full charge-transfer state of LiF (cation/anion pair) which can be automatically determined within FLOSIC. In conjunction with this method and the constrained charge-transfer method of Baruah and PedersonBaruah and Pederson (2009), it then becomes possible to calculate the lower-energy neutral state as well.

It was difficult to generate the curve in Fig. 5 in non-SC FLOSIC, since the calculation can converge to either the neutral Li and F atoms, or to Li+ and F- ions, and making sure the ground states at all separations are consistent is hard. This problem is alleviated in SC FLOSIC since we are able to solve for localized orbitals directly, so that it is easy to check the charged state of the atoms in each iteration.

IV Conclusion

SIC is a straightforward way to fix many of the flaws of semilocal xc energy functionals. Localized orbitals are a must to properly carry out SIC calculations. The localization equation Eq. (2) is the most general way of determining the localized orbitals, but it requires an overwhelming amount of calculation. FLOSIC simplifies the problem by reducing the number of parameters by an order of magnitude. Despite this simplification, the computational cost of FLOSIC is still high if one needs to solve Eq. (3) directly, since the Lagrange multiplier matrix $\lambda_{\sigma}$ is no longer diagonal as in KS equations, and the Hamiltonian for each orbital is different. Ref. Pederson et al. (2014) circumvents this problem by effectively converting the problem to a regular eigenvalue problem, but the results obtained this way are not fully self-consistent.

In this work we solve (3) with a novel Jacobi-like numerically exact iterative algorithm, and the results don’t suffer from the non-self-consistent problem associated with the method of Ref. Pederson et al. (2014). The newly developed algorithm only requires $N_{\sigma}\times M$ extra calculations of matrix elements compared with the KS equation, instead of $N_{\sigma}\times M^{2}$ as suggested in Eq. (3). Compared with the $N_{\sigma}\times N_{\sigma}$ extra matrix elements used in non-SC FLOSIC, we managed to include self-consistency with marginal increase in the computation cost.

The effect of the self-consistency is small for the total energies of the simple molecules tested in this work, and the accuracy of the FLOSIC method is confirmed. The changes in the SIC total energy from the non-SC values is about $1\%\sim 2\%$, and the changes in the atomization energies are smaller. The FODs mostly retain their positions as in the non-SC FLOSIC, but in some cases they show significant changes, indicating that properties depending on the density matrix may also have changed significantly. The atomization energies for both the non-SC and SC FLOSIC show vast improvement due to the removal of self-interaction, and higher accuracy can be expected from semilocal energy functionals specifically designed for SIC, or perhaps from strongly-constrained semilocal functionals like SCANSun et al. (2015).

This paper has developed a self-consistent approach for solving the FLOSIC equations in the representation of the FLO. We assumed the existence of canonical orbitals throughout the paper, but we must point out that, as discussed in Ref. Pederson et al. (2016), a rigorous derivation of the Kohn-Sham-like Hamiltonian for the canonical orbitals of the FLOSIC formulation is still needed. By observation, we find that Eq. (7) appears to be a good approximation to this Hamiltonian, but we also find that the magnitude of its eigenvalue for the highest occupied orbital overestimates the experimental ionization potential. This might be due to the approximation to the Hamiltonian, but it can also be due to the approximated energy functional. Further study is needed to solve this problem.

Acknowledgements

ZY is currently supported by Science Challenge Project No. TZ2016003 (China). ZY and JP were supported in part by the National Science Foundation(Grant Nos. DMR-1305135 and DMR-1607868, with contributions from Chemical Theory, Modeling and Computation). Support from the Office of Naval Research grant N00014-16-1-2464 graciously acknowledged.