Chebyshev expansion methods for electronic structure calculations on large molecular systems

Robustness of Local Predictions in Atomistic Machine Learning Models

Machine learning (ML) models for molecules and materials commonly rely on a decomposition of the global target quantity into local, atom-centered contributions. This approach is convenient from a computational perspective, enabling large-scale ML-driven simulations with a linear-scaling cost and also allows for the identification and posthoc interpretation of contributions from individual chemical environments and motifs to complicated macroscopic properties. However, even though practical justifications exist for the local decomposition, only the global quantity is rigorously defined. Thus, when the atom-centered contributions are used, their sensitivity to the training strategy or the model architecture should be carefully considered. To this end, we introduce a quantitative metric, which we call the local prediction rigidity (LPR), that allows one to assess how robust the locally decomposed predictions of ML models are. We investigate the dependence of the LPR on the aspects of model training, particularly the composition of training data set, for a range of different problems from simple toy models to real chemical systems. We present strategies to systematically enhance the LPR, which can be used to improve the robustness, interpretability, and transferability of atomistic ML models.

Forces from Stochastic Density Functional Theory under Nonorthogonal Atom-Centered Basis Sets

We develop a formalism for calculating forces on the nuclei within the linear-scaling stochastic density functional theory (sDFT) in a nonorthogonal atom-centered basis set representation (Fabian et al. Wiley Interdiscip. Rev.: Comput. Mol. Sci.2019, 9, e1412, 10.1002/wcms.1412) and apply it to the Tryptophan Zipper 2 (Trp-zip2) peptide solvated in water. We use an embedded-fragment approach to reduce the statistical errors (fluctuation and systematic bias), where the entire peptide is the main fragment and the remaining 425 water molecules are grouped into small fragments. We analyze the magnitude of the statistical errors in the forces and find that the systematic bias is of the order of 0.065 eV/Å (∼1.2 × 10^–3E_h/a₀) when 120 stochastic orbitals are used, independently of system size. This magnitude of bias is sufficiently small to ensure that the bond lengths estimated by stochastic DFT (within a Langevin molecular dynamics simulation) will deviate by less than 1% from those predicted by a deterministic calculation.

Energy window stochastic density functional theory

Linear scaling density functional theory is important for understanding electronic structure properties of nanometer scale systems. Recently developed stochastic density functional theory can achieve linear or even sublinear scaling for various electronic properties without relying on the sparsity of the density matrix. The basic idea relies on projecting stochastic orbitals onto the occupied space by expanding the Fermi-Dirac operator and repeating this for N_χ stochastic orbitals. Often, a large number of stochastic orbitals are required to reduce the statistical fluctuations (which scale as N_χ ^-1/2) below a tolerable threshold. In this work, we introduce a new stochastic density functional theory that can efficiently reduce the statistical fluctuations for certain observable and can also be integrated with an embedded fragmentation scheme. The approach is based on dividing the occupied space into energy windows and projecting the stochastic orbitals with a single expansion onto all windows simultaneously. This allows for a significant reduction of the noise as illustrated for bulk silicon with a large supercell. We also provide theoretical analysis to rationalize why the noise can be reduced only for a certain class of ground state properties, such as the forces and electron density.

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.

Self-averaging stochastic Kohn-Sham density-functional theory

We formulate the Kohn-Sham density functional theory (KS-DFT) as a statistical theory in which the electron density is determined from an average of correlated stochastic densities in a trace formula. The key idea is that it is sufficient to converge the total energy per electron to within a predefined statistical error in order to obtain reliable estimates of the electronic band structure, the forces on nuclei, the density and its moments, etc. The fluctuations in the total energy per electron are guaranteed to decay to zero as the system size increases. This facilitates "self-averaging" which leads to the first ever report of sublinear scaling KS-DFT electronic structure. The approach sidesteps calculation of the density matrix and thus, is insensitive to its evasive sparseness, as demonstrated here for silicon nanocrystals. The formalism is not only appealing in terms of its promise to far push the limits of application of KS-DFT, but also represents a cognitive change in the way we think of electronic structure calculations as this stochastic theory seamlessly converges to the thermodynamic limit.

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 calculation of an n-type GaAs quantum dot

A linear scale method for calculating electronic properties of large and complex systems is introduced within a local density approximation. The method is based on the Chebyshev polynomial expansion and the time-dependent method, which is tested on the calculation of the electronic structure of a model n-type GaAs quantum dot.

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.

Dual-Basis Analytic Gradients. 1. Self-Consistent Field Theory

Analytic gradients of dual-basis Hartree-Fock and density functional theory energies have been derived and implemented, which provide the opportunity for capturing large basis-set gradient effects at reduced cost. Suggested pairings for gradient calculations are 6-31G/6-31G**, dual[-f,-d]/cc-pVTZ, and 6-311G*/6-311 + +G(3df,3pd). Equilibrium geometries are produced within 0.0005 A of large-basis results for the latter two pairings. Though a single, iterative SCF response equation must be solved (unlike standard SCF gradients), it may be obtained in the smaller basis set, and integral screening further reduces the cost for well-chosen subsets. Total nuclear force calculations exhibit up to 75% savings, relative to large-basis calculations.

Advances in methods and algorithms in a modern quantum chemistry program package

Advances in theory and algorithms for electronic structure calculations must be incorporated into program packages to enable them to become routinely used by the broader chemical community. This work reviews advances made over the past five years or so that constitute the major improvements contained in a new release of the Q-Chem quantum chemistry package, together with illustrative timings and applications. Specific developments discussed include fast methods for density functional theory calculations, linear scaling evaluation of energies, NMR chemical shifts and electric properties, fast auxiliary basis function methods for correlated energies and gradients, equation-of-motion coupled cluster methods for ground and excited states, geminal wavefunctions, embedding methods and techniques for exploring potential energy surfaces.

Spin-unrestricted linear-scaling electronic structure theory and its application to magnetic carbon-doped boron nitride nanotubes

We present an extension of the density-matrix-based linear-scaling electronic structure theory to incorporate spin degrees of freedom. When the spin multiplicity of the system can be predetermined, the generalization of the existing linear-scaling methods to spin-unrestricted cases is straightforward. However, without calculations it is hard to determine the spin multiplicity of some complex systems, such as many magnetic nanostuctures and some inorganic or bioinorganic molecules. Here we give a general prescription to obtain the spin-unrestricted ground state of open-shell systems. Our methods are implemented into the linear-scaling trace-correcting density-matrix purification algorithm. The numerical atomic-orbital basis, rather than the commonly adopted Gaussian basis functions, is used. The test systems include O2 molecule and magnetic carbon-doped boron nitride (BN)(5,5) and BN(7,6) nanotubes. Using the newly developed method, we find that the magnetic moments in carbon-doped BN nanotubes couple antiferromagnetically with each other. Our results suggest that the linear-scaling spin-unrestricted trace-correcting purification method is very powerful to treat large magnetic systems.

"Lagrange functions": a family of powerful basis sets for real-space order-N electronic structure calculations

Plane waves have unparalleled simplicity and systematic convergence by a single monotonic parameter, the energy cutoff, but they are limited to speriodic systems and require Fourier transforms that scale as N(2)logN, where N is the number of atoms. Real-space methods for order-N scaling are computationally complex and convergence depends on several variables. Here we introduce and demonstrate "Lagrange functions" as a family of analytical, complete, and orthonormal basis sets that are suitable for efficient, accurate, real-space, order-N electronic structure calculations. Convergence is controlled by a single monotonic parameter, the dimension of the basis set, and computational complexity is lower than that of plane waves.

Sparse matrix multiplications for linear scaling electronic structure calculations in an atom-centered basis set using multiatom blocks

A sparse matrix multiplication scheme with multiatom blocks is reported, a tool that can be very useful for developing linear-scaling methods with atom-centered basis functions. Compared to conventional element-by-element sparse matrix multiplication schemes, efficiency is gained by the use of the highly optimized basic linear algebra subroutines (BLAS). However, some sparsity is lost in the multiatom blocking scheme because these matrix blocks will in general contain negligible elements. As a result, an optimal block size that minimizes the CPU time by balancing these two effects is recovered. In calculations on linear alkanes, polyglycines, estane polymers, and water clusters the optimal block size is found to be between 40 and 100 basis functions, where about 55-75% of the machine peak performance was achieved on an IBM RS6000 workstation. In these calculations, the blocked sparse matrix multiplications can be 10 times faster than a standard element-by-element sparse matrix package.