Edoardo Di Napoli

Abhiram Badrinarayanan, Davor Davidovic, Edoardo Di Napoli et al.

Simulating large molecular systems comprising thousands of atoms requires highly scalable methodologies. While modern Density Functional Theory (DFT) codes exhibit linear scaling, solving the associated large, sparse generalized eigenproblems remains a critical computational bottleneck on exascale architectures. In the context of the LimitX project, we propose a data-driven framework to accelerate these calculations. By shifting the machine learning target from discrete eigenvalues to the coefficients of an interpolating Chebyshev polynomial, and by comparing both all-atom and fragment-based structural representations, we successfully overcome the dimensionality constraints of large-scale spectral prediction. We investigate three machine learning models (Kernel Ridge Regression, Graph Neural Networks, and Random Forests) trained on a novel 2 TB dataset of protein dimers. The predicted spectra provide initial guesses that effectively bypass early Self-Consistent Field (SCF) iterations in BigDFT. Ultimately, these spectral predictors will be deployed to dynamically optimize upcoming rational filter-based eigensolvers, such as FrASE, which is currently in initial development.

Edoardo Di Napoli, Eric Polizzi, Yousef Saad

Estimating the number of eigenvalues located in a given interval of a large sparse Hermitian matrix is an important problem in certain applications and it is a prerequisite of eigensolvers based on a divide-and-conquer paradigm. Often an exact count is not necessary and methods based on stochastic estimates can be utilized to yield rough approximations. This paper examines a number of techniques tailored to this specific task. It reviews standard approaches and explores new ones based on polynomial and rational approximation filtering combined with a stochastic procedure.

Lukas Krämer, Edoardo Di Napoli, Martin Galgon et al.

We analyze the FEAST method for computing selected eigenvalues and eigenvectors of large sparse matrix pencils. After establishing the close connection between FEAST and the well-known Rayleigh-Ritz method, we identify several critical issues that influence convergence and accuracy of the solver: the choice of the starting vector space, the stopping criterion, how the inner linear systems impact the quality of the solution, and the use of FEAST for computing eigenpairs from multiple intervals. We complement the study with numerical examples, and hint at possible improvements to overcome the existing problems.

Edoardo Di Napoli, Xinzhe Wu, Thomas Bornhake et al.

In the last decade, the use of Machine and Deep Learning (MDL) methods in Condensed Matter physics has seen a steep increase in the number of problems tackled and methods employed. A number of distinct MDL approaches have been employed in many different topics; from prediction of materials properties to computation of Density Functional Theory potentials and inter-atomic force fields. In many cases the result is a surrogate model which returns promising predictions but is opaque on the inner mechanisms of its success. On the other hand, the typical practitioner looks for answers that are explainable and provide a clear insight on the mechanisms governing a physical phenomena. In this work, we describe a proposal to use a sophisticated combination of traditional Machine Learning methods to obtain an explainable model that outputs an explicit functional formulation for the material property of interest. We demonstrate the effectiveness of our methodology in deriving a new highly accurate expression for the enthalpy of formation of solid solutions of lanthanides orthophosphates.

Jan Winkelmann, Edoardo Di Napoli

Rational filter functions can be used to improve convergence of contour-based eigensolvers, a popular family of algorithms for the solution of the interior eigenvalue problem. We present a framework for the optimization of rational filters based on a non-convex weighted Least-Squares scheme. When used in combination with the FEAST library, our filters out-perform existing ones on a large and representative set of benchmark problems. This work provides a detailed description of: (1) a set up of the optimization process that exploits symmetries of the filter function for Hermitian eigenproblems, (2) a formulation of the gradient descent and Levenberg-Marquardt algorithms that exploits the symmetries, (3) a method to select the starting position for the optimization algorithms that reliably produces effective filters, (4) a constrained optimization scheme that produces filter functions with specific properties that may be beneficial to the performance of the eigensolver that employs them.

Edoardo Di Napoli, Clément Richefort, Xinzhe Wu

Studying the optoelectronic structure of materials can require the computation of several thousands of the smallest positive eigenpairs of a pseudo-hermitian Hamiltonian. Iterative eigensolvers may be preferred over direct methods for this task since their complexity is a function of the desired fraction of the spectrum. In addition, they generally rely on highly optimized and scalable kernels such as matrix-vector multiplications that leverage the massive parallelism and the computational power of modern exascale systems. The Chebyshev Accelerated Subspace iteration Eigensolver (ChASE) is able to compute several thousands of the most extreme eigenpairs of dense hermitian matrices with proven scalability over massive parallel accelerated clusters. This work presents an extension of ChASE to solve for a portion of the smallest positive eigenpairs of pseudo-hermitian Hamiltonians as they appear in the treatment of excitonic materials. By exploiting the numerical structure and spectral properties of the Hamiltonian matrix, we preserve the characteristic positive-negative symmetry in the treatment of the eigenvectors and propose an oblique variant of Rayleigh-Ritz projection that features quadratic convergence of the Ritz values with no explicit construction of the dual basis. Additionally, we introduce a parallel implementation of the recursive matrix-product operation appearing in the Chebyshev filter with limited amount of global communications. Our development is supported by a full numerical analysis and experimental tests.

Christoph Conrads, Edoardo Di Napoli

This paper presents our research software engineering (RSE) experiences over two years with libNEGF, a quantum transport code. We describe practical approaches to code quality assurance--including continuous integration, automated testing, and compiler warning correction--and performance engineering through continuous benchmarking. Our systematic application of these practices revealed critical defects: uninitialized memory reads, out-of-bounds writes, and notably, a misunderstood mathematical model in our boundary condition handling. We also document how continuous benchmarking exposed performance regressions caused by HPC system configuration changes. Our findings provide data points suggesting that a dangerous class of defects--equivalent to undefined behavior in C/C++ and processor-dependent behavior in Fortran--is as prevalent in Fortran scientific codes as elsewhere. While libNEGF is implemented in Fortran, most recommendations are applicable to scientific software regardless of implementation language, and they can be implemented selectively or in their entirety for both new and existing projects.

Edoardo Di Napoli, Alessandro Pecchia, Gustavo Ramirez-Hidalgo

The simulation of quantum transport in nanodevices requires the solution of the Dyson and Keldysh equations, a task dominated by the inversion of massive, block-tridiagonal matrices. While the Recursive Green's Function (RGF) method has long been the standard $O(N)$ solver for quasi-1D systems, its formulation has typically been restricted to sequential execution and nearest-neighbor interactions. In this work, we carefully reformulate RGF through the lens of Domain Decomposition and Schur Complement theory. This allows us to extend the recursive formalism to block $n$-diagonal systems (handling higher-order stencils) and to derive a parallel algorithm, Domain-Decomposition based RGF (DDRGF), which stitches macroscopic domains via reduced interface systems. We explore data dependencies in DDRGF in detail, by means of block-sparse structures and tracing back to the desired output as a block tridiagonal approximation, giving a clear, reproducible and extensible formulation. We validate these algorithms using \texttt{LibNEGF.jl}, a Julia-based implementation, demonstrating that the structural insights of domain decomposition provide a robust pathway for high-performance quantum transport simulations on modern multi-core clusters. The theory presented here lays down the base for tackling the Keldysh problem, to be similarly handled in future stages of our work. Although the target here is the acceleration of kernels in the non-equilibrium Green's function method, the algorithms and the implementations presented can be immediately used in any application involving block $n$-diagonal systems.

Edoardo Di Napoli, Xinzhe Wu

Chebyshev filtered subspace iteration is a well-known algorithm for the solution of (symmetric/Hermitian) algebraic eigenproblems which has been implemented in several application codes~\cite{Kronik:2006ff, abinit:2020} or in stand alone libraries~\cite{ChASE}. An essential part of the algorithm is the QR-factorization of the array of vectors spanning the active subspace that have been filtered by the Chebyshev filter. Typically such an array has an a-priori unknown high condition number that directly influences the choice of QR-factorization algorithm. In this work we show how such condition number can be bound from above with precise and inexpensive estimates. We then proceed to use these estimates to implement a mechanism for the choice of QR-factorization in the ChASE library. We show how such mechanism enhance the performance of the library without compromising on its accuracy.

Matthias Petschow, Edoardo Di Napoli, Paolo Bientinesi

One of the most used approaches in simulating materials is the tight-binding approximation. When using this method in a material simulation, it is necessary to compute the eigenvalues and eigenvectors of the Hamiltonian describing the system. In general, the system possesses few explicit symmetries. Due to them, the problem has many degenerate eigenvalues. The ambiguity in choosing a orthonormal basis of the invariant subspaces, associated with degenerate eigenvalues, will result in eigenvectors which are not invariant under the action of the symmetry operators in matrix form. A meaningful computation of the eigenvectors needs to take those symmetries into account. A natural choice is a set of eigenvectors, which simultaneously diagonalizes the Hamiltonian and the symmetry matrices. This is possible because all the matrices commute with each other. The simultaneous eigenvectors and the corresponding eigenvalues will be in a parametrized form in terms of the lattice momentum components. This functional dependence of the eigenvalues is the dispersion relation and describes the band structure of a material. Therefore it is important to find this functional dependence in any numerical computation related to material properties.