Black-box parameter selection and cost–accuracy trade-offs for contour-integral eigenstate extraction in O(N) density functional theory

Linear-scaling ($O(N)$) density functional theory, implemented in the CONQUEST code, can treat systems of $10^{4}$–$10^{5}$ atoms, and its 2017 combination with the contour-integral Sakurai–Sugiura eigensolver (SSM) made Kohn–Sham eigenstates of such systems computable — but only with eigensolver parameters (quadrature-point count $N_{\mathrm{q}}$, moment count $M$, source-block size $L_{\mathrm{s}}$, energy-window partitioning, quadrature rule) chosen manually per system. The open methodological problem left by Nakata et al. (J. Chem. Theory Comput. 13, 4146, 2017) is to establish a generally applicable, tuning-free selection of these parameters for arbitrary large systems with certified accuracy and predictable cost, together with a systematic quantification of the reconstruct-then-SSM trade-offs (primitive-basis reconstruction restores unoccupied-state accuracy at the price of a much larger Hamiltonian and costlier inner linear solves) and of hybrid MPI/OpenMP parallel best practice. No generally applicable scheme for this problem is known.

Background

Kohn–Sham density functional theory (DFT) computes the ground-state electronic structure of a material by solving a matrix eigenvalue problem for the single-particle (Kohn–Sham) Hamiltonian; in conventional implementations the dense diagonalization cost grows as the cube of system size, which restricts routine first-principles calculations to roughly a thousand atoms. Linear-scaling, or $O(N)$, methods avoid this barrier by exploiting the spatial decay of the single-particle density matrix (the matrix of orbital occupations) in systems with an energy gap, so that cost grows only linearly with atom count $N$; for metallic (gapless) systems the same decay of the density matrix is recovered at a finite electronic temperature, so that the linear-scaling premise extends to the material classes targeted here. The code CONQUEST implements $O(N)$ DFT with sparse matrices expressed in localized support functions, which are linear combinations of pseudo-atomic orbitals (PAOs, the atom-centered primitive basis); multi-site support functions (MSSFs) combine the PAOs of an atom and its neighbours into a compact basis that is accurate for occupied states while greatly reducing the matrix dimension. Two structural limitations remain: $O(N)$ methods deliver the total energy and density but none of the individual Kohn–Sham eigenstates, and the compact MSSF basis loses accuracy for unoccupied states — while eigenstates and unoccupied levels are exactly what is needed for density-of-states analysis or simulated scanning-tunnelling-microscopy (STM) spectra of large, complex systems. Eigenstates in a chosen finite energy window can be obtained without full diagonalization from interior-eigenvalue solvers based on contour integration, of which the Sakurai–Sugiura method (SSM) and FEAST are the main examples: the resolvent $(z\mathbf{S}-\mathbf{H})^{-1}$ (for the Hamiltonian $\mathbf{H}$ and overlap matrix $\mathbf{S}$ of a non-orthogonal basis) is integrated along a closed contour enclosing the window, the resulting complex moments span a small subspace containing the wanted eigenvectors, and the small projected problem is solved exactly. In practice the contour integral is replaced by numerical quadrature, and accuracy and cost are governed by user-set parameters: the number of quadrature points $N_{\mathrm{q}}$, the number of complex moments $M$, the number of columns $L_{\mathrm{s}}$ of the random source matrix, the partitioning of the target energy range into sub-windows, and the quadrature rule; heuristic guidelines exist but the choice is made case by case. Nakata et al. (2017) combined these pieces into one workflow: a CONQUEST calculation (self-consistent field minimization or $O(N)$ geometry optimization) in the compact MSSF basis; reconstruction of the Hamiltonian in the larger primitive-PAO basis using the converged MSSF charge density, which restores the accuracy of unoccupied states without a full primitive-basis self-consistent calculation; and a one-shot SSM solve for eigenstates in the window of interest. They demonstrated accuracy against dense ScaLAPACK diagonalization (eigenvalues of a 3,439-atom hydrated DNA system agreeing to below $10^{-9}$ eV), simulated STM images for P2 molecules adsorbed on a 3,074-atom Si(001) surface region with a 5.5–7.4-fold wall-clock advantage over the ScaLAPACK dense solver, and — beyond the reach of any dense solver — eigenstates near the Fermi level for Ge hut clusters (strained germanium islands with (105)-oriented facets that self-assemble on Si(001) during heteroepitaxy) containing 23,737 and 194,573 atoms, with Hamiltonian dimensions 213,633 and 1,751,157, computed in 146 s on 64 nodes and 2,399 s on 6,400 nodes of the K computer. In every case the SSM parameters were fixed empirically ($N_{\mathrm{q}}=16$, $M=8$, $L_{\mathrm{s}}=64$ or 100, equal-width windows, Chebyshev quadrature), with implementation details and parameter-setting guidance delegated to the prior methodology literature (Sakurai, Tadano and Futamura, 2013) rather than derived from the system under study. The hybrid parallel balance between MPI (the Message Passing Interface, which distributes work across compute nodes) and OpenMP (shared-memory threading within a node) differed between the compared solvers and was reported only as a benchmark description (4 MPI processes with 5 OpenMP threads per node for SSM against 1 process with 20 threads for the dense solver), and the reconstruction step enlarges the Hamiltonian (from 4,774 to 27,883 dimensions for the DNA system) and hence its memory and inner-linear-system cost; the detailed physical analysis of the huge hut-cluster systems was deferred to a separate future publication.

Problem Statement

The problem is to make the demonstrated CONQUEST + Sakurai–Sugiura workflow a generally applicable, black-box procedure for eigenstate extraction from $O(N)$ density functional theory, by resolving three linked aspects of one objective — that a qualified user can obtain all eigenstates in a stated energy interval of a system at or beyond the $10^{4}$-atom scale, to a stated accuracy, without expert per-system tuning. The three aspects are all required, and each is separately resolvable — completing one leaves the other two open:

(1) Generally applicable parameter selection. Given the converged sparse Hamiltonian and overlap matrices $(\mathbf{H},\mathbf{S})$ from an $O(N)$ DFT calculation, a target energy interval, and an accuracy tolerance, the number of quadrature points $N_{\mathrm{q}}$, the number of complex moments $M$, the source-matrix block size $L_{\mathrm{s}}$, the partitioning of the interval into sub-windows, and the quadrature rule must be chosen automatically, with a priori or a posteriori justification (error and cost bounds or validated models), such that all eigenpairs in the interval are certified to the requested accuracy at predictable cost and memory — for arbitrary large systems across material classes (insulating, semiconducting, metallic), not by the case-by-case manual tuning used in the 2017 demonstrations. This aspect is resolved when a single automatic rule fixes all five parameter choices from $(\mathbf{H},\mathbf{S})$, the target interval and the accuracy tolerance alone, justified by a rigorous error-and-cost bound or by a validated error-and-cost model, with the certified accuracy, cost and memory holding for insulating, semiconducting and metallic systems at and beyond the $10^{4}$-atom scale.

(2) Reconstruct-then-SSM trade-offs. The interplay among reconstructing the primitive-PAO-basis Hamiltonian (which restores unoccupied-state accuracy but enlarges the matrix dimension and the memory and inner-linear-system cost), the SSM parameter choices, and the achievable spectral resolution must be systematically quantified, so that it is determined when reconstruction is preferable to alternatives such as enlarged support-function bases. This aspect is resolved when that comparison is quantified — in accuracy, matrix dimension, memory, inner-linear-system cost and achievable spectral resolution, over representative systems — and yields an explicit criterion for when reconstruction is preferable.

(3) Parallel best practice. The hybrid MPI/OpenMP configuration and the distribution of quadrature points and energy windows over processes must be established as documented, machine-aware best practice, so that wall-clock and resource-efficiency comparisons are not artifacts of a single benchmark's parallel balance. This aspect is resolved when the configuration and distribution recommendations, and the measurements that support them, are set out so that those comparisons are reproducible and not tied to the benchmark hardware of the 2017 demonstrations.

The detailed physical analysis of the specific Ge hut-cluster systems (surface-versus-bulk state mapping, strain-dependent electronic features, spectroscopic signatures) that the source deferred to a separate publication is an application study outside this methodological problem.

References

1. Ayako Nakata, Yasunori Futamura, Tetsuya Sakurai, David R. Bowler, Tsuyoshi Miyazaki, "Efficient Calculation of Electronic Structure Using O(N) Density Functional Theory", Journal of Chemical Theory and Computation 13(9), 4146–4153 (2017). DOI: 10.1021/acs.jctc.7b00385, https://doi.org/10.1021/acs.jctc.7b00385 (arXiv:1704.03203, https://arxiv.org/abs/1704.03203)
2. Tetsuya Sakurai, Hiroshi Sugiura, "A projection method for generalized eigenvalue problems using numerical integration", Journal of Computational and Applied Mathematics 159(1), 119–128 (2003). DOI: 10.1016/S0377-0427(03)00565-X, https://doi.org/10.1016/S0377-0427(03)00565-X
3. Tetsuya Sakurai, Yasunori Futamura, Hiroto Tadano, "Efficient Parameter Estimation and Implementation of a Contour Integral-Based Eigensolver", Journal of Algorithms & Computational Technology 7(3), 249–269 (2013). DOI: 10.1260/1748-3018.7.3.249, https://doi.org/10.1260/1748-3018.7.3.249
4. Eric Polizzi, "Density-matrix-based algorithm for solving eigenvalue problems", Physical Review B 79, 115112 (2009). DOI: 10.1103/PhysRevB.79.115112, https://doi.org/10.1103/PhysRevB.79.115112
5. Sergiu Arapan, David R. Bowler, Tsuyoshi Miyazaki, "A linear scaling DFT study of the growth of a new {105} facet layer on a Ge hut cluster", arXiv:1510.00526 (2015). https://arxiv.org/abs/1510.00526
6. Ayako Nakata, Jack S. Baker, Shereif Y. Mujahed, Jack T. L. Poulton, Sergiu Arapan, Jianbo Lin, Zamaan Raza, Sushma Yadav, Lionel Truflandier, Tsuyoshi Miyazaki, David R. Bowler, "Large scale and linear scaling DFT with the CONQUEST code", The Journal of Chemical Physics 152(16), 164112 (2020). DOI: 10.1063/5.0005074, https://doi.org/10.1063/5.0005074 (arXiv:2002.07704, https://arxiv.org/abs/2002.07704)
7. Ayako Nakata, David R. Bowler, Tsuyoshi Miyazaki, "Large-Scale DFT Methods for Calculations of Materials with Complex Structures", Journal of the Physical Society of Japan 91(9), 091011 (2022). DOI: 10.7566/JPSJ.91.091011, https://doi.org/10.7566/JPSJ.91.091011 (arXiv:2204.12894, https://arxiv.org/abs/2204.12894)
8. Tsuyoshi Miyazaki, Ayako Nakata, David R. Bowler, "Large-Scale First-Principles Calculation Technique for Nanoarchitectonics: Local Orbital and Linear-Scaling DFT Methods with the CONQUEST Code", NIMS Monographs (Springer), 303–317 (2022). DOI: 10.1007/978-4-431-56912-1_18, https://doi.org/10.1007/978-4-431-56912-1_18
9. Takahiro Yano, Yasunori Futamura, Akira Imakura, Tetsuya Sakurai, "Performance evaluation of the Sakurai-Sugiura method with a block Krylov subspace linear solver for large dense Hermitian-definite generalized eigenvalue problems", JSIAM Letters 10, 77–80 (2018). DOI: 10.14495/jsiaml.10.77, https://doi.org/10.14495/jsiaml.10.77
10. Yasunori Futamura, Tetsuya Sakurai, "Efficient Contour Integral-based Eigenvalue Computation Using an Iterative Linear Solver with Shift-Invert Preconditioning", International Conference on High Performance Computing in Asia-Pacific Region (HPC Asia 2021), 90–99 (2021). DOI: 10.1145/3432261.3432269, https://doi.org/10.1145/3432261.3432269
11. Akira Imakura, Keiichi Morikuni, Akitoshi Takayasu, "Verified eigenvalue and eigenvector computations using complex moments and the Rayleigh–Ritz procedure for generalized Hermitian eigenvalue problems", Journal of Computational and Applied Mathematics 424, 114994 (2023). DOI: 10.1016/j.cam.2022.114994, https://doi.org/10.1016/j.cam.2022.114994 (arXiv:2110.01822, https://arxiv.org/abs/2110.01822)
12. Edoardo Di Napoli, Eric Polizzi, Yousef Saad, "Efficient estimation of eigenvalue counts in an interval", Numerical Linear Algebra with Applications 23(4), 674–692 (2016). DOI: 10.1002/nla.2048, https://doi.org/10.1002/nla.2048
13. Stefan Güttel, Eric Polizzi, Ping Tak Peter Tang, Gautier Viaud, "Zolotarev Quadrature Rules and Load Balancing for the FEAST Eigensolver", SIAM Journal on Scientific Computing 37(4), A2100–A2122 (2015). DOI: 10.1137/140980090, https://doi.org/10.1137/140980090
14. Vassilis Kalantzis, James Kestyn, Eric Polizzi, Yousef Saad, "Domain decomposition approaches for accelerating contour integration eigenvalue solvers for symmetric eigenvalue problems", Numerical Linear Algebra with Applications 25(5), e2154 (2018). DOI: 10.1002/nla.2154, https://doi.org/10.1002/nla.2154
15. William Dawson, Eisuke Kawashima, Laura E. Ratcliff, Muneaki Kamiya, Luigi Genovese, Takahito Nakajima, "Complexity reduction in density functional theory: Locality in space and energy", The Journal of Chemical Physics 158(16), 164114 (2023). DOI: 10.1063/5.0142652, https://doi.org/10.1063/5.0142652 (arXiv:2301.06313, https://arxiv.org/abs/2301.06313)