We present an accurate and efficient real-space formulation of the Hellmann–Feynman stress tensor for $O(N)$ Kohn–Sham density functional theory (DFT). While applicable at any temperature, the formulation is most efficient at high temperature where the Fermi–Dirac distribution becomes smoother and the density matrix becomes correspondingly more localized. We first rewrite the orbital-dependent stress tensor for real-space DFT in terms of the density matrix, thereby making it amenable to $O(N)$ methods. We then describe its evaluation within the $O(N)$ infinite-cell Clenshaw–Curtis Spectral Quadrature (SQ) method, a technique that is applicable to metallic and insulating systems, is highly parallelizable, becomes increasingly efficient with increasing temperature, and provides results corresponding to the infinite crystal without the need of Brillouin zone integration. We demonstrate systematic convergence of the resulting formulation with respect to SQ parameters to exact diagonalization results and show convergence with respect to mesh size to the established plane wave results. We employ the new formulation to compute the viscosity of hydrogen at 10^{6} K from Kohn–Sham quantum molecular dynamics, where we find agreement with previous more approximate orbital-free density functional methods.

## I. INTRODUCTION

Kohn–Sham density functional theory (DFT)^{1,2} is among the most widely used first principles methods for understanding and predicting materials properties. The tremendous popularity of DFT can be attributed to its favorable accuracy-to-cost ratio and simplicity relative to other such *ab initio* theories. In DFT simulations for condensed matter systems, a fundamental quantity of interest, in addition to the ground-state energy and atomic forces, is the second order Hellmann–Feynman stress tensor. The components of this tensor represent derivatives of the energy density with respect to the six independent homogeneous strains that can be applied to the system, evaluated at the electronic ground state. The stress tensor has many applications, including the determination of equilibrium lattice parameters, equation of state (EOS), and the calculation of shear viscosity from quantum molecular dynamics (QMD) simulations.^{3–5}

The Hellmann–Feynman stress tensor in DFT, which has its origins in the work of Slater^{6} and Janak,^{7} has been developed in the context of pseudopotential plane wave calculations,^{8–11} the linearized augmented plane wave (LAPW) method,^{12} the projector augmented-wave (PAW) method,^{13} atom-centered orbital bases,^{14,15} the finite-element method,^{16} and the pseudopotential real-space finite-difference method.^{17} However, a common feature of all these formulations is their expression in terms of Kohn–Sham orbitals, the calculation of which is associated with computational cost and memory requirements that scale as $O(N3)$ and $O(N2)$ with respect to the number of atoms, respectively. This critical cubic scaling arises due to the orthonormality constraint on the Kohn–Sham orbitals, which also limits scalability in the context of high performance parallel computing, severely limiting the length and time scales that can be reached.

In order to overcome the critical $O(N3)$ bottleneck, a number of $O(N)$ approaches have been developed (e.g., Refs. 18–20 and references therein), which circumvent the calculation of the Kohn–Sham orbitals by proceeding instead through the density matrix to determine quantities of interest, achieving linear scaling by exploiting the exponential decay of the density matrix for insulating systems and metallic systems at finite temperature.^{21–25} These efforts have culminated in a number of mature codes;^{14,26–34} however, important challenges remain. These include limitations of underlying basis sets, large prefactors, the need for additional computational parameters, subtleties in determining sufficient numbers and/or centers of localized orbitals, the calculation of accurate atomic forces, and large-scale parallelization.^{19,35} In particular, to the best of our knowledge, the stress tensor has not been formulated and implemented in the context of $O(N)$ DFT, especially for calculations at high temperature.

Kohn–Sham calculations at high temperature occur in a range of applications, including the study of warm dense matter and dense plasmas in laser experiments and the interiors of giant planets and stars.^{36–41} However, these calculations present significant challenges due to the substantially larger number and lesser locality of orbitals that must be computed. Consequently, $O(N3)$ and local-orbital based $O(N)$ methods have very large prefactors, which makes QMD calculations for even small systems intractable. The recently developed Spectral Quadrature (SQ) method^{42–44} addresses scaling with the number of atoms as well as temperature while retaining systematic convergence to standard diagonalization results for metals and insulators alike. In particular, due to the increased locality of electronic interactions and enhanced smoothness of the Fermi–Dirac function, the cost of the $O(N)$ SQ approach decreases rapidly as temperature is increased. Furthermore, it is particularly well suited to massive parallelization since a majority of the communication is localized to nearby processors, the pattern of which remains fixed throughout the calculation.

In this work, we present an accurate and efficient real-space formulation of the stress tensor for $O(N)$ Kohn–Sham calculations. While applicable at any temperature, the formulation is most efficient at high temperature where the Fermi–Dirac distribution becomes smoother and the density matrix becomes correspondingly more localized. Starting with the recently derived expression for the stress tensor in real-space DFT,^{17} we develop a formulation in terms of the density matrix and then describe its evaluation within the infinite-cell Clenshaw–Curtis SQ method. The framework is applicable to metallic and insulating systems, is highly parallelizable, and becomes more efficient as the temperature is increased, and stresses corresponding to the infinite crystal can be computed without Brillouin zone integration. We demonstrate the systematic convergence of the stress tensor with respect to quadrature order as well as truncation radius—the two key parameters of the SQ method in addition to mesh size—to exact diagonalization results. In addition, we show convergence with respect to mesh size to the established plane wave results. Finally, we employ the new formulation to compute the viscosity of hydrogen at 10^{6} K from Kohn–Sham QMD.

The remainder of this manuscript is organized as follows: In Sec. II, we summarize the orbital-dependent stress tensor formulation for real-space DFT. Next, we develop a density matrix formulation for the stress tensor in Sec. III and describe its evaluation using the real-space SQ method in Sec. IV. Finally, we verify the proposed framework in Sec. V and provide concluding remarks in Sec. VI.

## II. STRESS TENSOR IN $O(N3)$ REAL-SPACE DFT

Consider an orthorhombic unit cell Ω (Fig. 1), with nuclei positioned at **R** = {**R**_{1}, **R**_{2}, …, **R**_{N}} and a total of *N*_{e} valence electrons. The lattice vectors corresponding to Ω are *L*_{1}*ê*_{1}, *L*_{2}*ê*_{2}, and *L*_{3}*ê*_{3}, where *ê*_{1}, *ê*_{2}, and *ê*_{3} are the lattice/Cartesian unit vectors. For an infinitesimal homogeneous deformation that maps Ω to Ω^{F}, as illustrated in Fig. 1, the Hellmann–Feynman stress tensor can be defined as^{17}

where |Ω| denotes the measure of the unit cell, the superscript (.)^{F} represents quantities after perturbation of the unit cell with a deformation gradient having components *F*_{αβ}, $L$ denotes the Lagrangian, $\Psi ={\psi 1,\psi 2,\u2026,\psi Ns}$ is the collection of Kohn–Sham orbitals with occupations $g={g1,g2,\u2026,gNs}$, *ϕ* is the electrostatic potential,^{45,46} and $G$ signifies the electronic ground state for the undeformed system.

Neglecting spin and Brillouin zone integration, the real-space stress tensor for the choice of a semilocal exchange–correlation functional and a local formulation of the electrostatics^{47–49} can be written as^{17}

where $\sigma \alpha \beta I$, $\sigma \alpha \beta II$, $\sigma \alpha \beta III$, and $\sigma \alpha \beta IV$ are the contributions arising from the electronic kinetic energy *T*_{s}, exchange–correlation energy *E*_{xc}, nonlocal pseudopotential energy *E*_{nl}, and the electrostatic energy *E*_{el}, respectively,^{75}

The electron density *ρ* can itself be expanded as

the electrostatic potential *ϕ* is the solution of the Poisson equation,

and the energy terms *E*_{xc} and *E*_{nl} take the form

In the above equations, $\u2207x\alpha $ is the *α*th component of the gradient vector $\u2207x$, *δ*_{αβ} is the Kronecker delta function, *λ*_{n} are the eigenvalues corresponding to *ψ*_{n}, *ε*_{xc} is the sum of the exchange and correlation energy per particle of a uniform electron gas, $Vxc$ is the exchange–correlation potential, $\chi \u0303Jlm=\u2211J\u2032\chi J\u2032lm$ are the periodically mapped nonlocal projectors, with $\chi J\u2032lm$ representing the projectors associated with the *J*′th atom and the index *lm* running over all azimuthal and magnetic quantum numbers, *b* = *∑*_{I}*b*_{I} is the total pseudocharge density of the nuclei, with *b*_{I} being the pseudocharge density of the *I*th nucleus that generates the potential *V*_{I}, $Eself=12\u2211I\u222b\Omega bI(x,RI)VI(x,RI)\u2009dx$ is the self-energy associated with the pseudocharge densities, $\sigma \alpha \beta Ec$ is the stress tensor contribution arising from the energy correction due to the overlapping pseudocharges,^{17} the summation index *J* runs over all atoms in Ω, the summation index *J*′ runs over the *J*th atom and its periodic images, and the summation index *I* runs over all atoms in $R3$. It is worth noting that all quantities in the above equations correspond to the electronic ground state, i.e., after the solution of the Kohn–Sham nonlinear eigenproblem.

## III. STRESS TENSOR IN $O(N)$ REAL-SPACEDFT

The Hellmann–Feynman stress tensor presented in Sec. II has been formulated in terms of the Kohn–Sham orbitals *ψ*_{n}. Although the evaluation of the stress tensor by itself scales as $O(N)$, calculation of the orbitals scales as $O(N3)$ with respect to system size.^{50} Since $O(N)$ methods^{18–20} bypass the calculation of the orbitals, computing the truncated density matrix (directly or indirectly) instead, the stress tensor needs to be suitably reformulated in terms of this quantity. To do so, we note that the density matrix is of the form

whose diagonal entries are related to the electron density by

Indeed, the density matrix has exponential decay for insulators as well as metals at finite electronic temperature,^{21,25} allowing for its truncation during computation.

Among the various stress tensor contributions in Eqs. (3)–(6), only $\sigma \alpha \beta I$ and $\sigma \alpha \beta III$ have an explicit dependence on the orbitals. Therefore, we now rewrite them in terms of the density matrix,

and

The second equality for $\sigma \alpha \beta I$ is obtained using integration by parts in conjunction with the divergence theorem, whereas the third equality is obtained by a rearrangement of terms. The second and the third equalities for $\sigma \alpha \beta III$ are obtained by rearrangement of terms.

Using the above relations and Eqs. (2), (4), and (6), we arrive at the following reformulation for the Hellmann–Feynman stress tensor in terms of the density matrix:

The evaluation of the stress tensor by itself scales as $O(N)$ with respect to system size by virtue of the truncated nature of the density matrix. The overall scaling also becomes $O(N)$ when a linear-scaling method (e.g., Refs. 18–20 and references therein) is used to calculate the truncated density matrix. In particular, approaches that rely on the expansion of the Hamiltonian-dependent Fermi operator^{42,51–53} are particularly attractive for high-temperature calculations. Note that, proceeding along the lines of Ref. 17, it is straightforward to generalize the above formulation to include non-orthogonal systems and Brillouin zone integration. However, we focus on the orthogonal case here for simplicity since non-orthogonal cells and Brillouin zone integration are typically not required in $O(N)$ calculations, where target systems tend to be large. In addition, the stresses corresponding to the infinite crystal can be computed in the Spectral Quadrature (SQ) formalism^{42–44} without Brillouin zone integration.

## IV. STRESS TENSOR IN $O(N)$ REAL-SPACE SPECTRAL QUADRATURE METHOD

We now describe the evaluation of the stress tensor within the framework of the $O(N)$ real-space Spectral Quadrature (SQ) method.^{42–44} In this approach, all quantities of interest (in discrete form) are expressed as bilinear forms or sums of bilinear forms, which are then approximated by spatially localized quadrature rules. The method does not assume the existence of a bandgap and so is applicable to metallic and insulating systems alike. In addition, the technique is particularly well suited to scalable high performance computing and becomes more efficient as the temperature is increased, by virtue of the electronic interactions becoming more localized and the representation of the Fermi–Dirac function becoming more compact.

It is clear from Eq. (15) that the off-diagonal components of the density matrix are required for the calculation of the stress tensor. In such a situation, Clenshaw–Curtis SQ is significantly more efficient than the Gauss SQ variant, motivating its selection here.^{43} In particular, we choose the *infinite-cell* version of Clenshaw–Curtis SQ, wherein results corresponding to the infinite crystal are obtained without recourse to Brillouin zone integration or large supercells. Specifically, rather than employing Bloch boundary conditions for the orbitals on the unit cell Ω, zero-Dirichlet (or equivalently periodic) boundary conditions are prescribed at infinity, and the relevant components of the density matrix for spatial points within Ω are calculated using the nearsightedness principle.^{54} Indeed, the infinite-cell approach reduces to the standard Γ-point calculation when the size of the truncation region is smaller than the domain size, a situation common in large-scale $O(N)$ DFT simulations, particularly for those at high temperature.

Proceeding as in previous work,^{43,44} we discretize the domain Ω with a uniform grid containing *N*_{d} finite-difference nodes, the collection of which is denoted by *K*_{Ω}. Then, we partition Ω into *N*_{p} non-overlapping regions of equal size such that $\Omega =\u22c3p=1Np\Omega p$ and $K\Omega =\u22c3p=1NpK\Omega p$, where *N*_{p} is the total number of processors and $K\Omega p$ denotes the collection of grid points associated with Ω_{p}, the domain local to the *p*th processor. We define the nodal Hamiltonian $Hq\u2208RNc\xd7Nc$ for node *q* ∈ *K*_{Ω} as the restriction of the Hamiltonian to its *region of influence*—the cuboid of side 2*R*_{cut} centered at that point containing a total of *N*_{c} finite-difference nodes, with *R*_{cut} denoting the truncation radius of the density matrix beyond which the electronic interactions are ignored. Similarly, we define $wq\u2208RNc\xd71$, $\u2207h,q\u2261\u2009\u22071h,q\u2208RNc\xd7Nc,\u22072h,q\u2208RNc\xd7Nc,\u2009\u22073h,q\u2208RNc\xd7Nc\u2009$, $Vnl,qI\u2208RNc\xd7Nc$, and $Xh,q\u2261X1h,q\u2208RNc\xd7Nc,\u2009X2h,q\u2208RNc\xd7Nc,$ $X3h,q\u2208RNc\xd7Nc$ for node *q* ∈ *K*_{Ω} as the restriction to its region of influence of the standard basis vector, gradient matrices, nonlocal pseudopotential matrix of the *I*th atom, and the spatial location matrices of the grid points, respectively.

In the aforementioned framework, the stress tensor contributions that explicitly depend on the density matrix, i.e., $\sigma \alpha \beta I$ and $\sigma \alpha \beta III$, can be written as follows:

where *n*_{pl} is the order of the Clenshaw–Curtis quadrature, *T*_{j} denotes the Chebyshev polynomial of degree *j*, $H^q=(Hq\u2212\chi qI)/\zeta q$ is the scaled and shifted nodal Hamiltonian having spectrum in [−1, 1], with $I\u2208RNc\xd7Nc$ signifying the identity matrix, $\chi q=(\lambda qmax+\lambda qmin)/2$, and $\zeta q=(\lambda qmax\u2212\lambda qmin)/2$, where $\lambda qmax$ and $\lambda qmin$ denote the maximum and minimum eigenvalues of **H**_{q}, respectively, $Dp\u2032c$ is the set of all atoms in $R3$ whose nonlocal projectors overlap with the extended processor domain, i.e., processor domain extended by *R*_{cut} on each side, the summation with a prime indicates that the first term is halved, $cqj$ is the Chebyshev expansion coefficient of the Fermi–Dirac function,

where $\mu ^q=(\mu \u2212\chi q)/\zeta q$ is the scaled and shifted Fermi energy and $\sigma ^=\sigma /\zeta q$ is the scaled smearing, and $tqj\u2208RNc\xd71$ is the *q*th column of $Tj(H^q)$, determined using the three term recurrence relation for Chebyshev polynomials,

Thereafter, the expression for the stress tensor in Eq. (15) take the form

where |dΩ| denotes the volume associated with each grid point, $Dpb$ is the set of all atoms in $R3$ whose pseudocharges overlap with the processor domain Ω_{p}, *ρ*_{q} is the electron density at node *q*,

*∇*_{h} is the finite-difference approximation to the gradient with components ∇_{1h}, ∇_{2h}, and ∇_{3h}, and **x**_{βh,q} is the spatial coordinate of the node *q* in the *β* direction. Note that in certain instances, only the pressure is needed and not the entire stress tensor. In such cases, it is possible to compute the pressure by just taking the mean of the diagonal components of the stress tensor. However, a direct evaluation of the pressure—formulation presented in Appendix A—is slightly more efficient and does not involve second order derivatives of the density matrix, which can be beneficial in certain other real-space discretization schemes, such as finite elements^{45,55} and their mesh-free counterparts.^{46}

## V. RESULTS AND DISCUSSION

We have implemented the proposed formulation for the stress tensor in the SQDFT code.^{44} We now verify the accuracy and convergence of the developed framework and demonstrate its practical utility by calculating the viscosity of hydrogen under extreme conditions from quantum molecular dynamics (QMD). In all simulations, we use the Gauss SQ method to evaluate the electron density during the self-consistent field (SCF) iteration^{42,56} and the Clenshaw–Curtis SQ method to evaluate the Hellmann–Feynman atomic forces.^{43,44} In addition, we employ a twelfth-order accurate finite-difference discretization, the periodic Pulay method^{57} for acceleration of the self-consistent field (SCF) iteration, and the Alternating Anderson–Richardson (AAR) linear solver^{58,59} for calculation of the electrostatic potential, as well as for the application of the real-space Kerker preconditioner.^{60} Note that although we focus on high temperature calculations in the following, the proposed approach is equally applicable to calculations at ambient temperature, as shown in Appendix B.

### A. Accuracy and convergence

We verify the accuracy and convergence of the proposed formulation and implementation by considering three representative systems with atoms randomly perturbed by up to 15% of the equilibrium interatomic distance: (i) a 32-atom cell of face-centered cubic aluminum (Al) with lattice constant 7.78 bohrs, smearing *σ* = 4 eV (i.e., 46 418 K), and local-density approximation (LDA) exchange–correlation,^{2} (ii) a 64-atom cell of lithium hydride (LiH) with lattice constant 7.37 bohrs, smearing *σ* = 4 eV (i.e., 46418 K), and LDA exchange-correlation, and (iii) a 64-atom cell of diamond cubic carbon (C) with lattice constant 4.61 bohrs, smearing *σ* = 21.5 eV (i.e., 250 000 K), and generalized gradient approximation (GGA) exchange-correlation.^{61} For the aluminum and lithium hydride systems, we use Troullier–Martins pseudopotentials,^{62} whereas for the carbon system, we use an Optimized Norm-Conserving Vanderbilt (ONCV) pseudopotential.^{63}

First, we verify the convergence of the stress tensor with respect to the two additional parameters introduced by the SQ method, i.e., order of quadrature *n*_{pl} and truncation radius *R*_{cut}. We choose mesh-sizes of *h* = 0.65 bohr, 0.46 bohr, and 0.23 bohr for the aluminum, lithium hydride, and carbon systems, respectively. The values of *R*_{cut} and *n*_{pl} are chosen to be large enough so as to not influence the convergence behavior of *n*_{pl} and *R*_{cut}, respectively. We present the results so obtained in Fig. 2, with the reference corresponding to diagonalization-based values obtained by SPARC^{64,65} at the same mesh-size and with a 4 × 4 × 4 Monkhorst–Pack grid for Brillouin zone integration. The proposed formulation shows exponential convergence in the stress tensor with respect to both *n*_{pl} and *R*_{cut}. In particular, {*n*_{pl}, *R*_{cut}} ∼ {55, 6}, {*n*_{pl}, *R*_{cut}} ∼ {70, 5.6}, and {*n*_{pl}, *R*_{cut}} ∼ {30, 2.3} are sufficient to obtain 1% accuracy in the stresstensor—errors typical in production simulations—for the aluminum, lithium hydride, and carbon systems, respectively. Note that the convergence with *R*_{cut} is dictated by the smearing,^{66} with values of $O(30)$ bohr required for ambient-like conditions and $O(1\u22122)$ bohr for high temperature calculations with smearing $O(40\u221280)$ eV. The convergence rate with *n*_{pl} is determined by the smearing (linear dependence) and the mesh-size (inverse square root dependence).^{42} A detailed discussion on the convergence of the SQ method with respect to these parameters can be found in Refs. 42 and 66.

Next, we verify convergence of the stress tensor with respect to spatial discretization using highly converged diagonalization-based results from the plane wave code ABINIT^{67} as reference. Specifically, we utilize {*n*_{pl}, *R*_{cut}} = {150, 12}, {200, 12}, and {180, 6} for the aluminum, lithium hydride, and carbon systems, respectively. These parameters are sufficient to put the associated errors well below the mesh errors of interest, as evident from the results in Fig. 2. In ABINIT, we employ a plane wave cutoff of 60 Ha and a 4 × 4 × 4 Monkhorst–Pack grid for Brillouin zone integration, which translates to stresses that are converged to within 0.01%. We present the results so obtained in Fig. 3, from which it is clear that there is systematic convergence in the stress tensor computed using the proposed approach. In particular, stresses accurate to within 0.1% are readily obtained.

### B. High temperature QMD: Viscosity of hydrogen at 10^{6} K

Hydrogen is the most abundant element in the universe and its properties have important implications in astrophysics and planetary science as well as inertial confinement fusion (ICF) experiments where isotopes of deuterium and tritium serve as the fuel. In particular, material properties such as the mass diffusivity, thermal conductivity, and viscosity of warm dense hydrogen impact the onset of convection or turbulence in a hydrodynamic system, which, in turn, impacts the distribution of energy in a star, planet, or ICF capsule. These important quantities are very difficult to measure under extreme conditions of pressure and temperature. However, they can be accurately calculated from first principles QMD simulations.

In this work, we calculate the viscosity of hydrogen at a temperature of 10^{6} K and density of 2 g/cm^{3}. To do so, we consider a canonical ensemble of 64 hydrogen atoms and employ the LDA exchange–correlation functional, a local ONCV pseudopotential suitable for the target temperature, and the isokinetic Gaussian thermostat.^{68} In order to extract a statistical uncertainty of the final result, we average over QMD simulations corresponding to 10 different initial conditions for the atom positions and velocities, where each simulation has been run for more than 25 000 steps with a time step of 0.04 fs. Since macroscopic dynamics properties can be written as the time-integral of a microscopic time correlation function using Green–Kubo (GK) relations,^{69} we calculate the viscosity using the relation

where *k*_{B} is the Boltzmann constant, *T* is the temperature, ⟨·⟩ denotes the ensemble average, and *s*_{i} are the independent components of the deviatoric (i.e., traceless) stress tensor, i.e., *σ*_{12}, *σ*_{23}, *σ*_{31}, (*σ*_{11} − *σ*_{22})/2, and (*σ*_{22} − *σ*_{33})/2.^{70}

We present the results so obtained in Fig. 4. It is clear that the mean of the ensemble average for the independent components of the deviatoric stress tensor decays in around 50 fs, resulting in a viscosity of 164 ± 39 mPa s. This full Kohn–Sham DFT result is consistent with the recent OFDFT calculations of Sjostrom and Daligault^{71} (Fig. 5) while being free of the kinetic energy functional approximation inherent in OFDFT.

## VI. CONCLUDING REMARKS

In this work, we have presented an accurate and efficient real-space formulation for computation of the Hellmann–Feynman stress tensor in $O(N)$ Kohn–Sham DFT calculations. While applicable at any temperature, the formulation is most efficient at high temperature where the Fermi–Dirac distribution becomes smoother and the density matrix becomes correspondingly more localized. After rewriting the orbital-dependent real-space stress tensor in terms of the density matrix, we developed an approach for evaluating it using the infinite-cell Clenshaw–Curtis variant of the Spectral Quadrature (SQ) method. Notably, the developed framework is applicable to metallic and insulating systems, is highly parallelizable, becomes more efficient as the temperature is increased, and can be used to compute stresses corresponding to the infinite crystal without the need of Brillouin zone integration. We demonstrated that the resulting formulation converges systematically with respect to both polynomial order and localization radius to the exact diagonalization result and with respect to mesh size to the established plane wave results. Finally, we employed the new formulation to compute the viscosity of hydrogen at 10^{6} K from Kohn–Sham quantum molecular dynamics, where we found agreement with the previous more approximate orbital-free density functional methods.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

This work was performed, in part, under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Support from the Advanced Simulation and Computing/Physics and Engineering Models program at LLNL is gratefully acknowledged. P.S. also acknowledges the support of the U.S. National Science Foundation (NSF) under Grant No. 1663244. Time on the Quartz supercomputer was provided by the Computing Grand Challenge program at LLNL. We thank Donald Hamann for use of and assistance with a development version of the ONCVPSP code.

### APPENDIX A: PRESSURE FORMULATION, IMPLEMENTATION, AND VERIFICATION

In certain applications such as the calculation of equations of state,^{5} only the pressure is required and not the entire stress tensor. Although the pressure can be computed from diagonal elements of the stress tensor, here, we present a suitable reformulation that can not only improve the efficiency of the calculation, albeit slightly, but also make it more amenable to other real-space discretizations such as finite elements^{45,55} and their mesh-free counterparts.^{46} We now present the reformulation for the pressure, its evaluation within the infinite-cell Clenshaw–Curtis Spectral Quadrature (SQ) method^{43} as implemented in the SQDFT code,^{44} and its verification through selected examples.

The expression for the pressure in terms of the density matrix can be written as

where the Hamiltonian

is introduced to eliminate the kinetic energy contribution that involves explicit second order derivatives of the density matrix.

In the context of the infinite-cell Clenshaw–Curtis method described in Sec. IV, the expression for the pressure takes the form

where the derivation for the first term can be found in previous work.^{43} In particular, rather than performing a Chebyshev polynomial expansion for the density matrix, it is done for the product of the Hamiltonian with the density matrix. This bypasses the need to calculate second-order derivatives of the density matrix, thereby making the evaluation of the pressure slightly more efficient as well as more amenable to other real-space discretization schemes.

We now verify the accuracy and convergence of the proposed formulation for the pressure using the same examples, as in Sec. V A, where a detailed description of the tests can also be found. First, we check convergence of the pressure with respect to quadrature order *n*_{pl} and truncation radius *R*_{cut} in Figs. 6(a) and 6(b), respectively, with the values from the real-space code SPARC^{64,72} again used as reference. It is clear that there is exponential convergence in the pressure with respect to both parameters, the curves being very similar to those obtained for the stress tensor. Next, we check convergence with spatial discretization in Fig. 6(c), the results from the plane wave code ABINIT^{67} again used as reference. We observe that there is systematic convergence in the pressure similar to the stress tensor. Note that the convergence is faster for the pressure, likely a consequence of the stress tensor’s off-diagonal components containing mixed derivatives, which are found to generally cause slower convergence with discretization.^{73}

### APPENDIX B: VERIFICATION OF STRESS TENSOR FORMULATION AT AMBIENT CONDITIONS

We now demonstrate the accuracy of the proposed stress tensor formulation and implementation for Kohn–Sham calculations at ambient temperature. Specifically, we consider a randomly perturbed 4-atom aluminum system and choose a smearing of *σ* = 0.27 eV, typical of the values employed for metallic systems under ambient conditions to facilitate self-consistent convergence.^{67,74} In Fig. 7, we verify the convergence of the stress tensor with respect to *n*_{pl} and *R*_{cut} for mesh size *h* = 0.65 bohr, with the reference diagonalization results obtained by SPARC at the same mesh-size and a 8 × 8 × 8 Monkhorst–Pack grid for Brillouin zone integration. We employ *R*_{cut} = 28 bohrs and *n*_{pl} = 650 for the convergence with *n*_{pl} and *R*_{cut}, respectively. As observed previously for larger smearing values, there is exponential convergence with both parameters. Note that larger values of *n*_{pl} and *R*_{cut} are required due to the reduced smoothness of the Fermi–Dirac function and the decreased locality of the electronic interactions, respectively.

It is worth noting that for the system chosen here, Brillouin zone integration is critical for the accurate calculation of the stress tensor. For example, there is 29% error in values of the stress when a Γ-point calculation is performed using ABINIT/SPARC. This verifies the accuracy of the infinite-cell SQ method in calculating the stress tensor corresponding to the infinite crystal, without the need for Brillouin zone integration or large supercells. Also note that we do not show convergence with mesh size here, since given large enough *n*_{pl} and *R*_{cut}, the results are independent of temperature/smearing, and therefore, the same curves as in Fig. 3 are obtained. Overall, it can be concluded that the proposed stress tensor formulation and implementation are not restricted by the choice of smearing/temperature.

## REFERENCES

*FreeON*,

Additional contributions arising from each of the energy terms vanish when taken together, by virtue of the Hellmann–Feynman theorem.