Fast Generation of Pipek-Mezey Wannier Functions via the Co-Iterative Augmented Hessian Method

TL;DR

This work extends the co-iterative augmented Hessian (CIAH) approach to k-point sampling for Pipek--Mezey Wannier function localization, introducing k-CIAH. By leveraging an efficient Hessian--vector product, the method attains $O(N_k^2 n^3)$ scaling and demonstrates fast, robust convergence across diverse solids, offering 2–3× speedups over first-order k-space methods and orders of magnitude improvements over Γ-point CIAH when localizing thousands of orbitals. The PMWFs produced enable accurate Wannier interpolation of electronic bands, with interpolation errors well below 0.1 eV for representative systems and faster real-space Fock-matrix decay than non-PMWF bases. Overall, k-CIAH provides a scalable, reliable workflow for PMWF localization in periodic systems, with potential extensions to real-space formulations and other second-order periodic orbital-optimization problems.

Abstract

We report a $k$-point extension of the second-order co-iterative augmented Hessian (CIAH) algorithm, termed $k$-CIAH, for Pipek-Mezey (PM) localization of Wannier functions (WFs). By exploiting an efficient evaluation of the Hessian-vector product, $k$-CIAH achieves $O(N_k^2 n^3)$ scaling in both CPU time and memory, matching that of previously reported first-order $k$-space approaches while improving upon the $O(N_k^3 n^3)$ scaling of $Γ$-point CIAH, where $N_k$ denotes the number of $k$-points sampling the first Brillouin zone and $n$ characterizes the unit-cell size. Benchmark calculations on a diverse set of solids -- including insulators, semiconductors, metals, and surfaces -- demonstrate the fast and robust convergence of $k$-CIAH-based PMWF optimization, which yields an overall computational efficiency approximately 2-3 -fold higher than first-order $k$-space methods and orders of magnitude higher than $Γ$-point CIAH for localizing 1000-5000 orbitals. The quality of the resulting PMWFs is further validated by accurate electronic band structures obtained via PMWF-based Wannier interpolation.

Figures (3)

• Figure 1: Convergence of $k$-CIAH and $k$-BFGS measured by the decay of the norm of gradient for (A) SiO2, (B) CO/MgO(001), and (C) aluminum. The convergence threshold ($10^{-5}$) is denoted by the black horizontal line.
• Figure 2: CPU time (seconds) for PMWF optimization using $k$-CIAH, $k$-BFGS, and $\Gamma$-CIAH for (A) h-BN and (B) CO/MgO(001). The Brillouin zone is sampled with uniform $n\times n\times 1$$k$-meshes, with $n=9$--$35$ for h-BN and $n=2$--$8$ for CO/MgO(001). Fewer data points are shown for $\Gamma$-CIAH at large $n$ due to memory limitations. The $k$-CIAH and $k$-BFGS timings are fitted to an $O(N_{\textrm{orb}}^{2})$ scaling, whereas the $\Gamma$-CIAH data are fitted to an $O(N_{\textrm{orb}}^{3})$ scaling.
• Figure 3: (A--D) Electronic band structure of h-BN obtained from reference non-SCF calculations (gray) and from Wannier interpolation using PMWFs generated by $k$-CIAH (blue) for increasing SCF $k$-mesh sizes: (A) $3\times3$, (B) $5\times5$, (C) $7\times7$, and (D) $9\times9$. For comparison, bands interpolated using Kohn--Sham-based WFs without PM localization (KSWFs) are shown in red. In all cases, four occupied and two virtual bands are included in the interpolation. The high-symmetry points in the first Brillouin zone are indicated on the right. (E) Mean absolute error of the interpolated highest-occupied band (HOB; filled markers) and lowest-unoccupied band (LUB; hollow markers) as a function of SCF $k$-mesh size for h-BN. (F--G) Frobenius norm of the real-space Fock matrix $\|F(\bm{R})\|_{F}$ in the basis of (F) PMWFs and (G) KSWFs, visualized as two-dimensional heat maps for h-BN using a $21\times21$ SCF $k$-mesh. Black lines delineate the boundaries of the Wigner--Seitz cells.