Stress calculation in linear scaling DFT: convergence and dynamics

Abstract: We present the approach needed to calculate stress within density functional theory (DFT) using a localised orbital basis, both for exact diagonalisation and linear scaling approaches, and demonstrate our implementation within the large scale DFT code Conquest. For the linear scaling approach, we test the rate of convergence of stress with density matrix range, and compare it to the convergence of energy and forces for different materials with a range of band gaps. We show that excellent convergence is found for modest cutoffs, and show that large-scale isothermal-isobaric molecular dynamics is stable and accurate.

• The paper demonstrates that accurate stress calculation is attainable within linear scaling DFT using modest density matrix ranges, converging to within 0.1 GPa for most materials.
• It details the methodology for incorporating Hellmann-Feynman, Pulay, and additional grid-based contributions, clarifying the role of Lagrange multipliers in the stress tensor.
• The study validates stable and accurate NPT molecular dynamics simulations, enabling reliable ab initio simulations for systems with millions of atoms.

Stress Calculation in Linear Scaling Density Functional Theory: Convergence and Dynamics

Overview and Motivation

This work addresses the rigorous calculation of stress within density functional theory (DFT) for codes employing localised orbital basis sets, with a sharpened focus on both exact diagonalisation and linear scaling ($\mathcal{O}(N)$) approaches as implemented in the Conquest code. The authors systematically analyze the convergence properties of stress, energy, and forces as a function of the density matrix cutoff, with extensive benchmarks across materials with varying electronic structures and bonding motifs. The stability and accuracy of isothermal-isobaric ($NPT$) molecular dynamics under linear scaling DFT are also explicitly demonstrated.

Theoretical Formulation

Stress in DFT emerges as the derivative of the total energy with respect to strain; its evaluation within localized orbital frameworks requires precise treatment of both Hellmann-Feynman and Pulay contributions, alongside miscellaneous terms arising from grid-based integration and reciprocal space operations.

For $\mathcal{O}(N)$ DFT, the density matrix formalism is employed, with idempotency prescribed via the LNV method and electron number maintained using a Lagrange multiplier. Notably, the derivation clarifies that, while a crucial Lagrange-induced term in the force expression was previously (and justifiably) neglected for forces, it cannot be disregarded in the stress tensor, especially for modest density matrix cutoffs where deviations from perfect idempotency persist.

Pulay-like terms—as they relate to both overlap and basis function strain response—are carefully delineated. Supplementary stress terms address effects such as grid cell deformation and reciprocal lattice rescaling for Hartree potential evaluation, and explicit treatment is given for exchange-correlation (GGA) and non-linear partial core correction (PCC) modifications.

Numerical Analysis: Convergence of Stress, Energy, and Forces

Convergence tests use six archetype systems: elemental semiconductors (C, Si, Ge) and ionic crystals (MgO, LiCl, AlN), spanning insulating, semiconducting, and metallic regimes. For each, the authors provide an incisive comparison between the linear scaling and exact diagonalisation results, using a single-zeta PAO basis and a real-space grid of $0.2\,a_0$ spacing.

The main results show that the stress converges to within $0.1\,\mathrm{GPa}$ for a density matrix range of $20\,a_0$ across all systems except metallic Ge, where a slightly larger range is needed. Force errors are even smaller, sub-$5\times10^{-3}$ eV/Å for all practical cutoffs. Total energy convergence exhibits the anticipated exponential decay as a function of density matrix range, with faster convergence in wider gap materials.

Figure 1: Absolute value of the difference between linear scaling and fully converged exact diagonalisation calculations as a function of density matrix cutoff for (a) stress, (b) force, and (c) total energy; (d) stress for Si/Ge on a linear scale.

A notable physical insight is the rapid convergence of stress and force differences compared to total energy convergence, reinforcing the reliability of structural optimization and molecular dynamics even when the total energy itself converges more slowly.

Molecular Dynamics: Stability and Fidelity

The investigation is extended to $NPT$ molecular dynamics simulations at high pressure and ambient temperature, using bulk silicon as a prototype. The stability of the dynamics—i.e., conservation of the relevant extended Lagrangian quantity and maintenance of imposed thermodynamic variables—is maintained for density matrix cutoffs down to $12\,a_0$, with progressively better fidelity for larger cutoffs.

Crucially, even with a modest cutoff, agreement with exact diagonalisation trajectories is excellent over $NPT$0100 fs timescales provided the initial simulation cell pressure is matched. The authors elucidate that remaining trajectory divergence at longer scales stems primarily from initial pressure offsets rather than force inaccuracy, underlining the robustness of the method for large-scale DFT-based sampling.

Implications and Outlook

The work establishes that linear scaling DFT, as implemented here, affords reliable and accurate computation of stress—crucial for both cell optimization and $NPT$1 molecular dynamics—using manageable density matrix ranges. This enables routine $NPT$2-scaling simulations of systems with up to millions of atoms for both static and dynamic studies with full DFT accuracy.

Remaining technical challenges include the mitigation of issues related to overlap matrix inversion for large, more complete basis sets, essential for further improved accuracy and stability. Emerging approaches, such as alternative basis sets (blip, on-site support functions) and efficient extended Lagrangian methods requiring no self-consistency cycles, are highlighted as promising future directions for both efficiency and scalability.

Conclusion

This work provides a comprehensive analysis of stress calculation in linear scaling DFT, showcasing rapid and predictable convergence with modest density matrix ranges and validating stable, accurate molecular dynamics in the $NPT$3 ensemble. The results support the deployment of $NPT$4 DFT as a practical tool for large-scale atomistic simulations—including dynamical studies—without sacrificing formal accuracy or numerical stability. The trajectory towards ultra-large-scale, long time-scale ab initio molecular dynamics with quantum accuracy is thus open, contingent on continued algorithmic innovation in handling basis sets and computational efficiency.