Stable, Fast, and Accurate Kohn-Sham Inversion in Gaussian Basis for Open Shell Molecular and Condensed Phase Systems via Density Matrix Penalization

Abstract

Here we present a density matrix based KS inversion method formulated entirely within a Gaussian basis representation to optimize a KS potential matrix that reproduces a target electron density. Inverse Kohn-Sham (KS) density functional theory (DFT) aims to determine the effective local KS potential that reproduces a target electron density, and is important both for electronic structure analysis and for the development of orbital based correction methods. In finite Gaussian basis implementations, however, conventional inverse KS-DFT approaches such as the Zhao-Morrison-Parr (ZMP) method often become poorly constrained and inefficient, because the real space penalty potential is projected onto a limited number of Gaussian basis matrix elements, which can strongly coarse-grain its spatial variation. In the present method, the density matrix mismatch is defined in a Lowdin orthogonalized basis, which yields a penalty energy invariant under unitary rotations in that basis. The corresponding penalty potential contribution to the KS Hamiltonian is derived analytically in the original nonorthogonal Gaussian basis. Across a wide range of penalty strengths, the self consistent field (SCF) optimization remains robust and efficient for various open shell systems, while progressively tightening the penalty drives the electron density into accurate agreement with the target. Benchmarks on molecules and condensed phase systems show that the method achieves substantially smaller attainable density deviations than the conventional ZMP method. The method provides a fast and accurate route to KS inversion in finite Gaussian basis sets and may also be useful for future orbital based correction schemes.

Figures (4)

• Figure 1: Schematic illustration of the penalty strength iteration and KS-DFT SCF optimization procedure. The inner loop performs SCF iterations at fixed $\epsilon$ values, using an initial guess from either atomic calculations or the previously converged electron density and KS orbitals. The outer loop gradually decreases $\epsilon$ until satisfactory density matrix agreement is achieved.
• Figure 2: (a) The maximum absolute deviations among all matrix elements of $\Delta \mathbf{P}'^{\uparrow}$ and $\Delta \mathbf{P}'^{\downarrow}$, and (b) the maximum absolute deviations over all real space grid points of $\rho^{\uparrow}(\mathbf{r})-\rho_{\mathrm{target}}^{\uparrow}(\mathbf{r})$ and $\rho^{\downarrow}(\mathbf{r})-\rho_{\mathrm{target}}^{\downarrow}(\mathbf{r})$, for each test system at a given value of $\epsilon$. When $\epsilon$ is smaller than $1\times10^{-7}$, omitted bars indicate calculations that did not finish within the wall-time limit. We note that, for each test system, once $\epsilon$ becomes smaller than the value listed below, the SCF convergence criterion can no longer be satisfied within 6000 SCF iterations. The values shown in the figure therefore correspond to those obtained at the 6000th SCF iteration: $\mathrm{O}_2$ ($\epsilon=1\times10^{-10}$), $\mathrm{NO}$ ($\epsilon=1\times10^{-10}$), $\mathrm{CuCl}_2$ ($\epsilon=1\times10^{-11}$), $(\mathrm{C}_6\mathrm{H}_6)_2^{+}$ ($\epsilon=1\times10^{-11}$), rutile $\mathrm{TiO}_2(110)$ polaron ($\epsilon=1\times10^{-4}$), rutile $\mathrm{TiO}_2(110)+\mathrm{OH}$ ($\epsilon=1\times10^{-3}$), $\mathrm{NiO}$ (AFM) ($\epsilon=1\times10^{-11}$), $\mathrm{CoO}$ (AFM) ($\epsilon=1\times10^{-11}$), bulk $\mathrm{CeO}_2$ polaron ($\epsilon=1\times10^{-11}$), and $32\,\mathrm{H_2O}+\mathrm{Ag}^{2+}$ ($\epsilon=1\times10^{-11}$).
• Figure 3: (a) Smallest attainable maximum absolute deviation between the real space electron density and the target density on the grid for all test systems as $\epsilon$ decreases. (b) Corresponding $\epsilon$ values. Results obtained using the density matrix penalization method, the Coulomb-based ZMP method, and the Yukawa-based ZMP method are shown in green, blue, and yellow, respectively. Solid and dashed lines denote calculations that converged within 6000 SCF iterations and calculations that reached the maximum limit of 6000 iterations without satisfying the convergence criterion, respectively. Detailed numerical results for the Coulomb-based and Yukawa-based ZMP methods are listed in Tables S3 and S4.
• Figure 4: Number of SCF iterations required to achieve SCF convergence for all test systems in inverse KS-DFT calculations along the decreasing sequence of $\epsilon$. Results obtained using the density matrix penalization method, the Coulomb-based ZMP method, and the Yukawa-based ZMP method are shown in the top, middle, and bottom panels, respectively.