From Hartree Product to Kohn-Sham and Beyond: Exploring Self-Interaction in Self-Consistent Field Methods

Abstract

Self-interaction error (SIE) is a commonly known problem that most Kohn-Sham density functional theory (KS-DFT) approximate functionals display to varying extents. It originates from the incomplete cancellation of the Coulomb self-repulsion by the approximate exchange functionals. This is one of the major challenges for DFT, and therefore increasing our understanding of it could have great benefits for future use of DFT. Herein we advance techniques to dissect, understand, and textit{avoid} SIE in new ways. Considering that KS-DFT requires solving the self-consistent field (SCF) equations, we first present a robust and economical SCF solver - the "Quasi-Newton Unitary Optimization with Trust region"(QUOTR) solver.[Slattery, et. al. textit{Phys. Chem. Chem. Phys.}, textbf{2024}, 26, 6557-6573] To be robust, the solver is a direct-minimization solver equipped with a trust region (TR); to be economical, the solver uses an L-BFGS approximate Hessian and a physically-relevant preconditioner. Coupling these two aspects together is a solver for the TR subproblem that exploits the low-rank structure of the L-BFGS Hessian. We demonstrate that QUOTR is useful, not only for obtaining KS-DFT wave functions in difficult cases, but also for solving for Hartree-Fock (HF) orbitals in challenging chemical systems containing Cr or Fm. Although not able to beat the low cost of traditional Roothaan-Hall (RH) solvers with acceleration, QUOTR is robust in its convergence at only a modest increase in computational cost. The many examples of SCF convergence problems when using semi-local KS-DFT functionals are known to be the result of a vanishing HOMO-LUMO gap, which is further the result of SIE. A major motivation for developing QUOTR came from our desire to understand the "true" (albeit unphysical) ground state solutions in cases where KS-DFT could not be converged by a traditional diagonalization-based SCF solver. We reinvestigate the relationship between the vanishing HOMO-LUMO gap and SCF non-convergence using our QUOTR solver. A set of difficult biological systems that had previously been shown to display convergence problems [Rudberg, et. al. textit{J. Phys.: Condens. Matter}, textbf{2012}, 24, 072202] was selected for deeper analysis. In addition to being able to obtain converged solutions, we analyze the resulting densities matrices in comparison to HF. The source of the vanishing HOMO-LUMO gaps is demonstrated to be incompatible eigenspectrums of spatially distant fragments in the peptides. We show that by using a local solver (QUOTR) with an appropriate initial guess, that a non-Aufbau filled stationary point can be found for vacuum-separated charged fragments. A systematic scan of all 20 natural amino acids for some common DFAs is used to examine the prevalence of predicted non-Aufbau filling. We find that hybrid functionals improve upon GGAs more than meta-generalized gradient approximations (GGAs) do, and range-separated functionals are much better - though not completely solving the problem. Having addressed SIE in biomolecules in terms of where charges comes from and where it goes, we finally analyze SIE on an orbital-by-orbital basis. We define the textit{genuine} exchange energy as the difference between the HF energy and the (self-interaction free) orthogonal Hartree product wave function energy. We propose that the Edmiston-Ruedenberg [Edmiston, et. al. textit{Rev. Mod. Phys.}, textbf{1963}, 35, 457-464] localized HF orbitals are the most appropriate HF frame for this analysis, due to their connection with the Hartree product wave function. Although the use of HF orbitals to quantify the genuine exchange energy is an approximation, we demonstrate that the error of total exchange energy is approximately 10-15% for a set of small molecules. The good performance of two popular GGAs is shown to arise with considerable error cancellation between orbitals (particularly core and valence). We also examine two orbital-dependent DFAs: the Perdew-Zunger self-interaction correction (PZ-SIC), and a generalization of the Hartree-Fock-Gopinathan (HFG) method.