Trace resetting density matrix purification in O(N) self-consistent-field theory

Calculation of Vibrational Circular Dichroism Spectra Employing Nuclear Velocity Perturbation or Magnetic Field Perturbation Theory Using an Atomic-Orbital-Based Linear Response Approach

We present the implementation of an AO-based solver for vibrational circular dichroism (VCD) spectra calculations, employing nuclear velocity perturbation theory using velocity-gauge included in atomic orbitals or magnetic field perturbation theory using gauge-including atomic orbitals. The implementations are done within the Gaussian and plane waves framework in the CP2K package. The previously implemented approaches in CP2K are based on MO-solvers, which solve the Sternheimer equation for linear response. Our AO-solver implementations were validated against the MO-solver by performing VCD calculations for the R-enantiomer of mirtazapine. Additionally, we extended the AO-based solver implementation to nuclear displacement perturbation theory in order to calculate infrared absorption spectra. The AO-based solver produced spectra that matched exactly with the MO-based results, confirming the accuracy of the implementation. This allows for efficient calculations of vibrational properties, further extending the capabilities of CP2K for large molecular systems.

Reducing Numerical Precision Requirements in Quantum Chemistry Calculations

The abundant demand for deep learning compute resources has created a renaissance in low-precision hardware. Going forward, it will be essential for simulation software to run on this new generation of machines without sacrificing scientific fidelity. In this paper, we examine the precision requirements of a representative kernel from quantum chemistry calculations: the calculation of the single-particle density matrix from a given mean-field Hamiltonian (i.e., Hartree-Fock or density functional theory) represented in an LCAO basis. We find that double precision affords an unnecessarily high level of precision, leading to optimization opportunities. We show how an approximation built from an error-free matrix multiplication transformation can be used to potentially accelerate this kernel on future hardware. Our results provide a roadmap for adapting quantum chemistry software for the next generation of high-performance computing platforms.

Susceptibility formulation of density matrix perturbation theory

Density matrix perturbation theory based on recursive Fermi-operator expansions provides a computationally efficient framework for time-independent response calculations in quantum chemistry and materials science. From a perturbation in the Hamiltonian, we can calculate the first-order perturbation in the density matrix, which then gives us the linear response in the expectation values for some chosen set of observables. We present an alternative, dual formulation, where we instead calculate the static susceptibility of an observable, which then gives us the linear response in the expectation values for any number of different Hamiltonian perturbations. We show how the calculation of the susceptibility can be performed with the same expansion schemes used in recursive density matrix perturbation theory, including generalizations to fractional occupation numbers and self-consistent linear response calculations, i.e., similar to density functional perturbation theory. As with recursive density matrix perturbation theory, the dual susceptibility formulation is well suited for numerically thresholded sparse matrix algebra, which has linear scaling complexity for sufficiently large sparse systems. Similarly, the recursive computation of the susceptibility also seamlessly integrates with the computational framework of deep neural networks used in artificial intelligence (AI) applications. This integration enables the calculation of quantum response properties that can leverage cutting-edge AI-hardware, such as NVIDIA Tensor Cores or Google Tensor Processing Units. We demonstrate performance for recursive susceptibility calculations using NVIDIA Graphics Processing Units and Tensor Cores.

Graph-based quantum response theory and shadow Born-Oppenheimer molecular dynamics

Graph-based linear scaling electronic structure theory for quantum-mechanical molecular dynamics simulations [A. M. N. Niklasson et al., J. Chem. Phys. 144, 234101 (2016)] is adapted to the most recent shadow potential formulations of extended Lagrangian Born-Oppenheimer molecular dynamics, including fractional molecular-orbital occupation numbers [A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020) and A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021)], which enables stable simulations of sensitive complex chemical systems with unsteady charge solutions. The proposed formulation includes a preconditioned Krylov subspace approximation for the integration of the extended electronic degrees of freedom, which requires quantum response calculations for electronic states with fractional occupation numbers. For the response calculations, we introduce a graph-based canonical quantum perturbation theory that can be performed with the same natural parallelism and linear scaling complexity as the graph-based electronic structure calculations for the unperturbed ground state. The proposed techniques are particularly well-suited for semi-empirical electronic structure theory, and the methods are demonstrated using self-consistent charge density-functional tight-binding theory both for the acceleration of self-consistent field calculations and for quantum-mechanical molecular dynamics simulations. Graph-based techniques combined with the semi-empirical theory enable stable simulations of large, complex chemical systems, including tens-of-thousands of atoms.

A variational formulation of the Harris functional as a correction to approximate Kohn-Sham density functional theory

Accurate descriptions of intermolecular interactions are of great importance in simulations of molecular liquids. We present an electronic structure method that combines the accuracy of the Harris functional approach with the computational efficiency of approximately linear-scaling density functional theory (DFT). The non-variational nature of the Harris functional has been addressed by constructing a Lagrangian energy functional, which restores the variational condition by imposing stationarity with respect to the reference density. The associated linear response equations may be treated with linear-scaling efficiency in an atomic orbital based scheme. Key ingredients to describe the structural and dynamical properties of molecular systems are the forces acting on the atoms and the stress tensor. These first-order derivatives of the Harris Lagrangian have been derived and implemented in consistence with the energy correction. The proposed method allows for simulations with accuracies close to the Kohn-Sham DFT reference. Embedded in the CP2K program package, the method is designed to enable ab initio molecular dynamics simulations of molecular solutions for system sizes of several thousand atoms. Available subsystem DFT methods may be used to provide the reference density required for the energy correction at near linear-scaling efficiency. As an example of production applications, we applied the method to molecular dynamics simulations of the binary mixtures cyclohexane-methanol and toluene-methanol, performed within the isobaric-isothermal ensemble, to investigate the hydrogen bonding network in these non-ideal mixtures.

Q-Next: A Fast, Parallel, and Diagonalization-Free Alternative to Direct Inversion of the Iterative Subspace

As computer systems dedicated to scientific calculations become massively parallel, the poor parallel performance of the Fock matrix diagonalization becomes a major impediment to achieving larger molecular sizes in self-consistent field (SCF) calculations. In this Article, a novel, highly parallel, and diagonalization-free algorithm for the accelerated convergence of the SCF procedure is presented. The algorithm, called Q-Next, draws on the second-order SCF, quadratically convergent SCF, and direct inversion of the iterative subspace (DIIS) approaches to enable fast convergence while replacing the Fock matrix diagonalization SCF bottleneck with higher parallel efficiency matrix multiplications. Performance results on both parallel multicore CPU and GPU hardware for a variety of test molecules and basis sets are presented, showing that Q-Next achieves a convergence rate comparable to the DIIS method while being, on average, one order of magnitude faster.

Subsystem Density Functional Theory Augmented by a Delta Learning Approach to Achieve Kohn-Sham Accuracy

Simulations based on electronic structure theory naturally include polarization and have no transferability problems. In particular, Kohn-Sham density functional theory (KS-DFT) has become the method of reference for ab initio molecular dynamics simulations of condensed matter systems. However, the high computational cost often poses strict limits on the affordable system size as well as on the extension of sampling (number of configurations). In this work, we propose an improvement to the subsystem density functional theory approach, known as the Kim-Gordon (KG) scheme, thus enabling the sampling of configurations for condensed molecular systems keeping the KS-DFT level accuracy at a fraction of computer time. Our scheme compensates the known KG shortcomings of the electronic kinetic energy term by adding a simple correction and can match KS-DFT accuracy in energies and forces. The computationally cheap correction is determined by means of a machine learning procedure. The proposed KG scheme is applied within a linear scaling self-consistent field formalism and is assessed by a series of molecular dynamics simulations of liquid water under different conditions. Although system-dependent, the correction is transferable between system sizes and temperatures.

Parallel Implementation of Large-Scale Linear Scaling Density Functional Theory Calculations With Numerical Atomic Orbitals in HONPAS

Linear-scaling density functional theory (DFT) is an efficient method to describe the electronic structures of molecules, semiconductors, and insulators to avoid the high cubic-scaling cost in conventional DFT calculations. Here, we present a parallel implementation of linear-scaling density matrix trace correcting (TC) purification algorithm to solve the Kohn-Sham (KS) equations with the numerical atomic orbitals in the HONPAS package. Such a linear-scaling density matrix purification algorithm is based on the Kohn's nearsightedness principle, resulting in a sparse Hamiltonian matrix with localized basis sets in the DFT calculations. Therefore, sparse matrix multiplication is the most time-consuming step in the density matrix purification algorithm for linear-scaling DFT calculations. We propose to use the MPI_Allgather function for parallel programming to deal with the sparse matrix multiplication within the compressed sparse row (CSR) format, which can scale up to hundreds of processing cores on modern heterogeneous supercomputers. We demonstrate the computational accuracy and efficiency of this parallel density matrix purification algorithm by performing large-scale DFT calculations on boron nitrogen nanotubes containing tens of thousands of atoms.

CP2K: An electronic structure and molecular dynamics software package - Quickstep: Efficient and accurate electronic structure calculations

CP2K is an open source electronic structure and molecular dynamics software package to perform atomistic simulations of solid-state, liquid, molecular, and biological systems. It is especially aimed at massively parallel and linear-scaling electronic structure methods and state-of-the-art ab initio molecular dynamics simulations. Excellent performance for electronic structure calculations is achieved using novel algorithms implemented for modern high-performance computing systems. This review revisits the main capabilities of CP2K to perform efficient and accurate electronic structure simulations. The emphasis is put on density functional theory and multiple post-Hartree-Fock methods using the Gaussian and plane wave approach and its augmented all-electron extension.

Machine Learning Adaptive Basis Sets for Efficient Large Scale Density Functional Theory Simulation

It is chemically intuitive that an optimal atom centered basis set must adapt to its atomic environment, for example by polarizing toward nearby atoms. Adaptive basis sets of small size can be significantly more accurate than traditional atom centered basis sets of the same size. The small size and well conditioned nature of these basis sets leads to large saving in computational cost, in particular in a linear scaling framework. Here, it is shown that machine learning can be used to predict such adaptive basis sets using local geometrical information only. As a result, various properties of standard DFT calculations can be easily obtained at much lower costs, including nuclear gradients. In our approach, a rotationally invariant parametrization of the basis is obtained by employing a potential anchored on neighboring atoms to ultimately construct a rotation matrix that turns a traditional atom centered basis set into a suitable adaptive basis set. The method is demonstrated using MD simulations of liquid water, where it is shown that minimal basis sets yield structural properties in fair agreement with basis set converged results, while reducing the computational cost in the best case by a factor of 200 and the required flops by 4 orders of magnitude. Already a very small training set yields satisfactory results as the variational nature of the method provides robustness.

Efficient Computation of Sparse Matrix Functions for Large-Scale Electronic Structure Calculations: The CheSS Library

We present CheSS, the "Chebyshev Sparse Solvers" library, which has been designed to solve typical problems arising in large-scale electronic structure calculations using localized basis sets. The library is based on a flexible and efficient expansion in terms of Chebyshev polynomials and presently features the calculation of the density matrix, the calculation of matrix powers for arbitrary powers, and the extraction of eigenvalues in a selected interval. CheSS is able to exploit the sparsity of the matrices and scales linearly with respect to the number of nonzero entries, making it well-suited for large-scale calculations. The approach is particularly adapted for setups leading to small spectral widths of the involved matrices and outperforms alternative methods in this regime. By coupling CheSS to the DFT code BigDFT, we show that such a favorable setup is indeed possible in practice. In addition, the approach based on Chebyshev polynomials can be massively parallelized, and CheSS exhibits excellent scaling up to thousands of cores even for relatively small matrix sizes.

Communication: Generalized canonical purification for density matrix minimization

A Lagrangian formulation for the constrained search for the N-representable one-particle density matrix based on the McWeeny idempotency error minimization is proposed, which converges systematically to the ground state. A closed form of the canonical purification is derived for which no a posteriori adjustment on the trace of the density matrix is needed. The relationship with comparable methods is discussed, showing their possible generalization through the hole-particle duality. The appealing simplicity of this self-consistent recursion relation along with its low computational complexity could prove useful as an alternative to diagonalization in solving dense and sparse matrix eigenvalue problems.

O(N) methods in electronic structure calculations

Linear-scaling methods, or O(N) methods, have computational and memory requirements which scale linearly with the number of atoms in the system, N, in contrast to standard approaches which scale with the cube of the number of atoms. These methods, which rely on the short-ranged nature of electronic structure, will allow accurate, ab initio simulations of systems of unprecedented size. The theory behind the locality of electronic structure is described and related to physical properties of systems to be modelled, along with a survey of recent developments in real-space methods which are important for efficient use of high-performance computers. The linear-scaling methods proposed to date can be divided into seven different areas, and the applicability, efficiency and advantages of the methods proposed in these areas are then discussed. The applications of linear-scaling methods, as well as the implementations available as computer programs, are considered. Finally, the prospects for and the challenges facing linear-scaling methods are discussed.

Assessment of density matrix methods for linear scaling electronic structure calculations

Purification and minimization methods for linear scaling computation of the one-particle density matrix for a fixed Hamiltonian matrix are compared. This is done by considering the work needed by each method to achieve a given accuracy in terms of the difference from the exact solution. Numerical tests employing orthogonal as well as non-orthogonal versions of the methods are performed using both element magnitude and cutoff radius based truncation approaches. It is investigated how the convergence speed for the different methods depends on the eigenvalue distribution in the Hamiltonian matrix. An expression for the number of iterations required for the minimization methods studied is derived, taking into account the dependence on both the band gap and the chemical potential. This expression is confirmed by numerical tests. The minimization methods are found to perform at their best when the chemical potential is located near the center of the eigenspectrum. The results indicate that purification is considerably more efficient than the minimization methods studied even when a good starting guess for the minimization is available. In test calculations without a starting guess, purification is more than an order of magnitude more efficient than minimization.

Linear scaling density matrix perturbation theory for basis-set-dependent quantum response calculations: An orthogonal formulation

Linear scaling density matrix perturbation theory [A. M. N. Niklasson and M. Challacombe, Phys. Rev. Lett. 92, 193001 (2004)] is extended to basis-set-dependent quantum response calculations for a nonorthogonal basis set representation. The generalization is achieved by a perturbation-dependent congruence transform, derived from the factorization of the inverse overlap matrix, which transforms the generalized eigenvalue problem to an orthogonal, standard form. With this orthogonalization transform the basis-set-dependent perturbation in the overlap matrix is included in the orthogonalized Hamiltonian, which is expanded in orders of the perturbation. In this way density matrix perturbation theory developed for an orthogonal representation can be applied also to basis-set-dependent response calculations. The method offers an alternative to the previous solution of the basis-set-dependent response problem, based on a nonorthogonal generalization of the density matrix perturbation theory, where the calculations are performed within a purely nonorthogonal setting [A. M. N. Niklasson et al., J. Chem. Phys. 123, 44107 (2005)].

Linear scaling calculation of band edge states and doped semiconductors

Linear scaling methods provide total energy, but no energy levels and canonical wave functions. From the density matrix computed through the density matrix purification methods, we propose an order-N [O(N)] method for calculating both the energies and wave functions of band edge states, which are important for optical properties and chemical reactions. In addition, we also develop an O(N) algorithm to deal with doped semiconductors based on the O(N) method for band edge states calculation. We illustrate the O(N) behavior of the new method by applying it to boron nitride (BN) nanotubes and BN nanotubes with an adsorbed hydrogen atom. The band gap of various BN nanotubes are investigated systematically and the acceptor levels of BN nanotubes with an isolated adsorbed H atom are computed. Our methods are simple, robust, and especially suited for the application in self-consistent field electronic structure theory.

A hierarchic sparse matrix data structure for large-scale Hartree-Fock/Kohn-Sham calculations

A hierarchic sparse matrix data structure for Hartree-Fock/Kohn-Sham calculations is presented. The data structure makes the implementation of matrix manipulations needed for large systems faster, easier, and more maintainable without loss of performance. Algorithms for symmetric matrix square and inverse Cholesky decomposition within the hierarchic framework are also described. The presented data structure is general; in addition to its use in Hartree-Fock/Kohn-Sham calculations, it may also be used in other research areas where matrices with similar properties are encountered. The applicability of the data structure to ab initio calculations is shown with help of benchmarks on water droplets and graphene nanoribbons.