Diagonalization without Diagonalization: A Direct Optimization Approach for Solid-State Density Functional Theory

Abstract

We present a novel approach to address the challenges of variable occupation numbers in direct optimization of density functional theory (DFT). By parameterizing both the eigenfunctions and the occupation matrix, our method minimizes the free energy with respect to these parameters. As the stationary conditions require the occupation matrix and the Kohn-Sham Hamiltonian to be simultaneously diagonalizable, this leads to the concept of ``self-diagonalization,'' where, by assuming a diagonal occupation matrix without loss of generality, the Hamiltonian matrix naturally becomes diagonal at stationary points. Our method incorporates physical constraints on both the eigenfunctions and the occupations into the parameterization, transforming the constrained optimization into an fully differentiable unconstrained problem, which is solvable via gradient descent. Implemented in JAX, our method was tested on aluminum and silicon, confirming that it achieves efficient self-diagonalization, produces the correct Fermi-Dirac distribution of the occupation numbers and yields band structures consistent with those obtained with SCF methods in Quantum Espresso.

Figures (4)

• Figure 1: The self-diagonalization process of the Kohn-Sham Hamiltonian matrix $\hat{h}$ for an aluminum crystal at the $\Gamma$ point, at a temperature of $T=0.01$ Ha. Left: A visualization of the Hamiltonian matrix at steps 100, 300, 500, 1000, 2000 and 5000. The color represents the matrix elements of $\hat{h}$: blue indicates values below zero down to $-5$ Ha, red indicates values above zero up to $5$ Ha, and white indicates close to zero, as represented in the colorbar. The plot on the right shows how the variance of the off-diagonal elements of $\hat{h}$ decreases as the optimization progresses. It is evident that as the free energy (Eq. \ref{['eq:free_energy']}) is minimized, $\hat{h}$ gradually self-diagonalizes.
• Figure 2: A visualization of the change of occupation numbers during the optimization for a FCC aluminum crystal. Top: Occupation numbers as a function of Hamiltonian diagonal matrix elements relative to the Fermi level ($\varepsilon_i({\bf k}) - \varepsilon_{\text{fermi}}$). Each point represents a diagonal element-occupation number pair for a potentially occupied orbital. The color indicates the value of the occupation number, with red representing 1 and blue representing 0. A theoretical Fermi-Dirac distribution is shown as a dotted line. The occupation number distributions are displayed at steps 0, 1000, 2000, and 3000. The rightmost figure focuses on the eigenvalues near the Fermi level within a narrow energy range (x-axis) at step 3000. Bottom: An illustration of the band structure of a metal (aluminum based on an FCC conventional unit-cell) and its relation to the occupation numbers.
• Figure 3: Comparison of the electronic band structures of Aluminum and Silicon calculated using the proposed method and Quantum Espresso. We use $3\times 3 \times 3$ k-point mesh, cutoff energy of $100$ Ha, smearing/temperature values of $T=0.01Ha$ for both materials. As our method is an all-electron method, we have tuned the pseudo-potentials in our Quantum Espresso calculations to allow for an all-electron calculation. Quantum Espresso employs the conventional self-consistent field (SCF) Kohn-Sham DFT with LDA_X and a plane-wave basis in the pseudopotential projector-augmented wave formalism. The pseudopotentials used in Quantum Espresso are modified such that all-electron calculations are conducted to match our method for evaluation purposes.
• Figure 4: Scaling comparison between the proposed method and full Hamiltonian matrix method. The x-axis represents the number of orbitals (dimension of the Hamiltonian matrix), while the y-axis shows the total training time of 1000 iterations.