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 106 K from Kohn–Sham quantum molecular dynamics, where we find agreement with previous more approximate orbital-free density functional methods.

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 Slater6 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) method42–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 106 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.

Consider an orthorhombic unit cell Ω (Fig. 1), with nuclei positioned at R = {R1, R2, …, RN} and a total of Ne valence electrons. The lattice vectors corresponding to Ω are L1ê1, L2ê2, and L3ê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 as17 

σαβ=1|Ω|LF(Ψ,g,ϕ,RF)FαβG,α,β{1,2,3},
(1)

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, Ψ={ψ1,ψ2,,ψNs} is the collection of Kohn–Sham orbitals with occupations g={g1,g2,,gNs}, ϕ is the electrostatic potential,45,46 and G signifies the electronic ground state for the undeformed system.

FIG. 1.

Orthorhombic unit cell Ω (solid lines) with sides of length L1, L2, and L3 along the x1, x2, and x3 directions, respectively. The unit cell subsequent to the application of the infinitesimal homogeneous deformation is denoted by ΩF (dashed lines).

FIG. 1.

Orthorhombic unit cell Ω (solid lines) with sides of length L1, L2, and L3 along the x1, x2, and x3 directions, respectively. The unit cell subsequent to the application of the infinitesimal homogeneous deformation is denoted by ΩF (dashed lines).

Close modal

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 electrostatics47–49 can be written as17 

σαβ=1|Ω|σαβI+σαβII+σαβIII+σαβIV,
(2)

where σαβI, σαβII, σαβIII, and σαβIV are the contributions arising from the electronic kinetic energy Ts, exchange–correlation energy Exc, nonlocal pseudopotential energy Enl, and the electrostatic energy Eel, respectively,75 

σαβI=2n=1NsgnΩxαψn(x)xβψn(x)dx,
(3)
σαβII=δαβExc(ρ,ρ)ΩVxcρ(x),xρ(x)ρ(x)dxΩρ(x)εxcρ(x),xρ(x)xβρ(x)xαρ(x)dx,
(4)
σαβIII=δαβEnl(Ψ,g,R)4n=1NsgnJlmγJl×JΩχJlm(x,RJ)xRJβxαψn(x)dx×Ωχ̃Jlm(x,RJ)ψn(x)dx,
(5)
σαβIV=14πΩxαϕ(x,R)xβϕ(x,R)dx+IΩxαbI(x,RI)xRIβ×ϕ(x,R)12VI(x,RI)dx12IΩxαVI(x,RI)xRIβbI(x,RI)dx+12δαβΩb(x,R)ρ(x)ϕ(x,R)dxδαβEself(R)+σαβEc.
(6)

The electron density ρ can itself be expanded as

ρ(x)=2n=1Nsgnψn2(x),
(7)

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

14π2ϕ(x,R)=ρ(x,R)+b(x,R),
(8)

and the energy terms Exc and Enl take the form

Exc(ρ,ρ)=Ωεxcρ(x),xρ(x)ρ(x)dx,
(9)
Enl(Ψ,g,R)=2n=1NsgnJlmγJlΩχ̃Jlm(x,RJ)ψn(x)dx2.
(10)

In the above equations, xα is the αth component of the gradient vector x, δαβ 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, χ̃Jlm=JχJlm are the periodically mapped nonlocal projectors, with χJlm representing the projectors associated with the J′th atom and the index lm running over all azimuthal and magnetic quantum numbers, b = IbI is the total pseudocharge density of the nuclei, with bI being the pseudocharge density of the Ith nucleus that generates the potential VI, Eself=12IΩbI(x,RI)VI(x,RI)dx is the self-energy associated with the pseudocharge densities, σαβ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 Jth 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.

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) methods18–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

D(x,y)=n=1Nsgnψn(x)ψn(y),
(11)

whose diagonal entries are related to the electron density by

ρD(x)=2D(x,x).
(12)

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 σαβI and σαβIII have an explicit dependence on the orbitals. Therefore, we now rewrite them in terms of the density matrix,

σαβI=2n=1NsgnΩxαψn(x)xβψn(x)dx=2n=1NsgnΩψn(x)xαxβψn(x)dx=2Ωyαyβn=1Nsgnψn(y)ψn(x)y=xdx=2ΩyαyβD(y,x)y=xdx
(13)

and

σαβIII=δαβ2n=1NsgnJlmγJlΩχ̃Jlm(x,RJ)ψn(x)dxΩχ̃Jlm(y,RJ)ψn(y)dy4n=1NsgnJlmγJlJΩχJlm(x,RJ)xRJβxαψn(x)dxΩχ̃Jlm(y,RJ)ψn(y)dy=δαβ2JlmγJlΩΩχ̃Jlm(x,RJ)n=1Nsgnψn(x)ψn(y)χ̃Jlm(y,RJ)dxdy4JlmγJlJΩΩχJlm(x,RJ)xRJβxαn=1Nsgnψn(x)ψn(y)χ̃Jlm(y,RJ)dxdy=δαβ2JlmγJlΩΩχ̃Jlm(x,RJ)D(x,y)χ̃Jlm(y,RJ)dxdy4JlmγJlJΩΩχJlm(x,RJ)xRJβxαD(x,y)χ̃Jlm(y,RJ)dxdy.
(14)

The second equality for σαβ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 σαβ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:

σαβ=1|Ω|2ΩyαyβD(y,x)y=xdx+δαβExc(ρD,ρD)ΩVxcρD(x),xρD(x)ρD(x)dx,ΩρD(x)εxcρD(x),xρD(x)xβρD(x)xαρD(x)dxδαβEnl(D,R)4JlmγJlJΩΩχJlm(x,RJ)xRJβxαD(x,y)χ̃Jlm(y,RJ)dxdy+14πΩxαϕ(x,R)xβϕ(x,R)dx+IΩxαbI(x,RI)xRIβϕ(x,R)12VI(x,RI)dx12IΩxαVI(x,RI)xRIβbI(x,RI)dx+12δαβΩb(x,R)ρD(x)ϕ(x,R)dxδαβEself(R)+σαβEc.
(15)

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 operator42,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) formalism42–44 without Brillouin zone integration.

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 Nd finite-difference nodes, the collection of which is denoted by KΩ. Then, we partition Ω into Np non-overlapping regions of equal size such that Ω=p=1NpΩp and KΩ=p=1NpKΩp, where Np is the total number of processors and KΩp denotes the collection of grid points associated with Ωp, the domain local to the pth processor. We define the nodal Hamiltonian HqRNc×Nc for node qKΩ as the restriction of the Hamiltonian to its region of influence—the cuboid of side 2Rcut centered at that point containing a total of Nc finite-difference nodes, with Rcut denoting the truncation radius of the density matrix beyond which the electronic interactions are ignored. Similarly, we define wqRNc×1, h,q1h,qRNc×Nc,2h,qRNc×Nc,3h,qRNc×Nc, Vnl,qIRNc×Nc, and Xh,qX1h,qRNc×Nc,X2h,qRNc×Nc,X3h,qRNc×Nc for node qKΩ as the restriction to its region of influence of the standard basis vector, gradient matrices, nonlocal pseudopotential matrix of the Ith 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., σαβI and σαβIII, can be written as follows:

σαβI2p=1NpqKΩpwqTαh,qβh,qj=0nplcqjTj(H^q)wq=2p=1NpqKΩpwqTαh,qβh,qj=0nplcqjtqj,
(16)
σαβIII2δαβp=1NpIDpcqKΩpwqTVnl,qIj=0nplcqjTj(H^q)wq4p=1NpIDpcqKΩpwqTVnl,qIXβh,qRI,βIαh,qj=0nplcqjTj(H^q)wq=2p=1NpIDpcqKΩpwqTVnl,qIδαβj=0nplcqjtqj+2Xβh,qRI,βIαh,qj=0nplcqjtqj,
(17)

where npl is the order of the Clenshaw–Curtis quadrature, Tj denotes the Chebyshev polynomial of degree j, H^q=(HqχqI)/ζq is the scaled and shifted nodal Hamiltonian having spectrum in [−1, 1], with IRNc×Nc signifying the identity matrix, χq=(λqmax+λqmin)/2, and ζq=(λqmaxλqmin)/2, where λqmax and λqmin denote the maximum and minimum eigenvalues of Hq, respectively, Dpc is the set of all atoms in R3 whose nonlocal projectors overlap with the extended processor domain, i.e., processor domain extended by Rcut 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,

cqj=2π11g(r,μ^q,σ^q)Tj(r)1r2dr,
(18)

where μ^q=(μχq)/ζq is the scaled and shifted Fermi energy and σ^=σ/ζq is the scaled smearing, and tqjRNc×1 is the qth column of Tj(H^q), determined using the three term recurrence relation for Chebyshev polynomials,

tqi+1=2H^qtqitqi1,i=1,2,,npltq1=H^qwq,tq0=wq.
(19)

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

σαβ=1|Ω||dΩ|p=1NpqKΩp2|dΩ|wqTαh,qβh,qj=0nplcqjtqj+δαβϵxc(ρq)ρqVxc(ρq)ρqρqεxcρ,hρβhρqαhρq2p=1NpIDpcqKΩpwqTVnl,qIδαβj=0nplcqjtqj+2Xβh,qRI,βIαh,qj=0nplcqjtqj+|dΩ|p=1NpqKΩp14παhϕq×βhϕq+12δαβbqρqϕq+|dΩ|p=1NpIDpbqKΩpαhbIqxβh,qRI,βϕq12VI,q12αhVIq×xβh,qRI,βbI,q12bI,qVI,q+σαβEc,
(20)

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,

ρq=2|dΩ|j=0nplcqjρqj,ρqj=wqTtqj,
(21)

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 elements45,55 and their mesh-free counterparts.46 

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) iteration42,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 method57 for acceleration of the self-consistent field (SCF) iteration, and the Alternating Anderson–Richardson (AAR) linear solver58,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.

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 npl and truncation radius Rcut. 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 Rcut and npl are chosen to be large enough so as to not influence the convergence behavior of npl and Rcut, respectively. We present the results so obtained in Fig. 2, with the reference corresponding to diagonalization-based values obtained by SPARC64,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 npl and Rcut. In particular, {npl, Rcut} ∼ {55, 6}, {npl, Rcut} ∼ {70, 5.6}, and {npl, Rcut} ∼ {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 Rcut is dictated by the smearing,66 with values of O(30) bohr required for ambient-like conditions and O(12) bohr for high temperature calculations with smearing O(4080) eV. The convergence rate with npl 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.

FIG. 2.

Convergence of the stress tensor with respect to (a) quadrature order npl and (b) truncation radius Rcut for aluminum, lithium hydride, and carbon systems. The error is defined to be the maximum difference in any component from the corresponding results obtained by SPARC.

FIG. 2.

Convergence of the stress tensor with respect to (a) quadrature order npl and (b) truncation radius Rcut for aluminum, lithium hydride, and carbon systems. The error is defined to be the maximum difference in any component from the corresponding results obtained by SPARC.

Close modal

Next, we verify convergence of the stress tensor with respect to spatial discretization using highly converged diagonalization-based results from the plane wave code ABINIT67 as reference. Specifically, we utilize {npl, Rcut} = {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.

FIG. 3.

Convergence of the stress tensor with mesh size for aluminum, lithium hydride, and carbon systems. The error is defined to be the maximum difference in any component from the corresponding results obtained by ABINIT.

FIG. 3.

Convergence of the stress tensor with mesh size for aluminum, lithium hydride, and carbon systems. The error is defined to be the maximum difference in any component from the corresponding results obtained by ABINIT.

Close modal

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 106 K and density of 2 g/cm3. 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

η(t)=|Ω|kBT0t15i=15si(t)si(0)dt,
(22)

where kB is the Boltzmann constant, T is the temperature, ⟨·⟩ denotes the ensemble average, and si 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 Daligault71 (Fig. 5) while being free of the kinetic energy functional approximation inherent in OFDFT.

FIG. 4.

Ensemble average (top) and viscosity (bottom) for hydrogen at 106 K and density 2 g/cm3.

FIG. 4.

Ensemble average (top) and viscosity (bottom) for hydrogen at 106 K and density 2 g/cm3.

Close modal
FIG. 5.

Comparison between OFDFT71 and SQDFT for the viscosity of hydrogen at a density of 2 g/cm3.

FIG. 5.

Comparison between OFDFT71 and SQDFT for the viscosity of hydrogen at a density of 2 g/cm3.

Close modal

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 106 K from Kohn–Sham quantum molecular dynamics, where we found agreement with the previous more approximate orbital-free density functional methods.

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

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.

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 elements45,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) method43 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

P=13σ11+σ22+σ33=13|Ω|4ΩDH(x,x)dx+3Exc(ρD,ρD)ΩVxcρD(x),xρD(x)+εxcρD(x),xρD(x)xρD(x)xρD(x)×ρD(x)dxEnl(D,R)4JlmγJlJΩΩχJlm(x,RJ)xRJxD(x,y)χ̃Jlm(y,RJ)dxdy+14πΩxϕ(x,R)2dx+IΩxbI(x,RI)xRIϕ(x,R)12VI(x,RI)dx12IΩxVI(x,RI)xRIbI(x,RI)dx+12ΩρD(x)+3b(x,R)ϕ(x,R)dx3Eself(R)+i=13σiiEc,
(A1)

where the Hamiltonian

H=122+Vxc+ϕ+Vnl
(A2)

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

P=13|Ω||dΩ|p=1NpqKΩp4|dΩ|j=0nplχqcqj+ζqdqjρqj+3ϵxc(ρq)ρqVxc(ρq)ρqρqεxcρ,hρhρqhρq2p=1NpIDpcqKΩpwqTVnl,qIj=0nplcqjtqj+2Xh,qRIIh,qj=0nplcqjtqj+|dΩ|p=1NpqKΩp14πhϕq2+12ρq+3bqϕq+|dΩ|p=1NpIDpbqKΩphbIqxh,qRIϕq12VI,q12hVIq.xh,qRIbI,q32bI,qVI,q+i=13σiiEc,
(A3)

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 npl and truncation radius Rcut in Figs. 6(a) and 6(b), respectively, with the values from the real-space code SPARC64,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 ABINIT67 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 

FIG. 6.

Convergence of the pressure with respect to (a) quadrature order npl, (b) truncation radius Rcut, and (c) mesh size for aluminum, lithium hydride, and carbon systems. The error in (a) and (b) is defined to be the difference from the corresponding results obtained by SPARC, whereas in (c), it is the difference from the corresponding results obtained by ABINIT.

FIG. 6.

Convergence of the pressure with respect to (a) quadrature order npl, (b) truncation radius Rcut, and (c) mesh size for aluminum, lithium hydride, and carbon systems. The error in (a) and (b) is defined to be the difference from the corresponding results obtained by SPARC, whereas in (c), it is the difference from the corresponding results obtained by ABINIT.

Close modal

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 npl and Rcut 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 Rcut = 28 bohrs and npl = 650 for the convergence with npl and Rcut, respectively. As observed previously for larger smearing values, there is exponential convergence with both parameters. Note that larger values of npl and Rcut are required due to the reduced smoothness of the Fermi–Dirac function and the decreased locality of the electronic interactions, respectively.

FIG. 7.

Convergence of the stress tensor with respect to (a) quadrature order npl and (b) truncation radius Rcut for aluminum under ambient conditions. The error in the stress tensor is defined to be the maximum difference in any component from the corresponding results obtained by SPARC.

FIG. 7.

Convergence of the stress tensor with respect to (a) quadrature order npl and (b) truncation radius Rcut for aluminum under ambient conditions. The error in the stress tensor is defined to be the maximum difference in any component from the corresponding results obtained by SPARC.

Close modal

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 npl and Rcut, 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.

1.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
2.
W.
Kohn
and
L. J.
Sham
,
Phys. Rev.
140
,
A1133
(
1965
).
3.
G. A.
de Wijs
,
G.
Kresse
,
L.
Vočadlo
,
D.
Dobson
,
D.
Alfè
,
M. J.
Gillan
, and
G. D.
Price
,
Nature
392
,
805
(
1998
).
4.
N.
Jakse
and
A.
Pasturel
,
Sci. Rep.
3
,
3135
(
2013
).
5.
S.
Zhang
,
A.
Lazicki
,
B.
Militzer
,
L. H.
Yang
,
K.
Caspersen
,
J. A.
Gaffney
,
M. W.
Däne
,
J. E.
Pask
,
W. R.
Johnson
,
A.
Sharma
 et al,
Phys. Rev. B
99
,
165103
(
2019
).
6.
J. C.
Slater
,
J. Chem. Phys.
57
,
2389
(
1972
).
7.
J. F.
Janak
,
Phys. Rev. B
9
,
3985
(
1974
).
8.
M. T.
Yin
,
Phys. Rev. B
27
,
7769
(
1983
).
9.
O. H.
Nielsen
and
R. M.
Martin
,
Phys. Rev. B
32
,
3780
(
1985
).
10.
O. H.
Nielsen
and
R. M.
Martin
,
Phys. Rev. B
32
,
3792
(
1985
).
11.
A. D.
Corso
and
R.
Resta
,
Phys. Rev. B
50
,
4327
(
1994
).
12.
T.
Thonhauser
,
C.
Ambrosch-Draxl
, and
D. J.
Singh
,
Solid State Commun.
124
,
275
(
2002
).
13.
M.
Torrent
,
F.
Jollet
,
F.
Bottin
,
G.
Zérah
, and
X.
Gonze
,
Comput. Mater. Sci.
42
,
337
(
2008
).
14.
J. M.
Soler
,
E.
Artacho
,
J. D.
Gale
,
A.
García
,
J.
Junquera
,
P.
Ordejón
, and
D.
Sánchez-Portal
,
J. Phys.: Condens. Matter
14
,
2745
(
2002
).
15.
F.
Knuth
,
C.
Carbogno
,
V.
Atalla
,
V.
Blum
, and
M.
Scheffler
,
Comput. Phys. Commun.
190
,
33
(
2015
).
16.
P.
Motamarri
and
V.
Gavini
,
Phys. Rev. B
97
,
165132
(
2018
).
17.
A.
Sharma
and
P.
Suryanarayana
,
J. Chem. Phys.
149
,
194104
(
2018
).
18.
S.
Goedecker
,
Rev. Mod. Phys.
71
,
1085
(
1999
).
19.
D. R.
Bowler
and
T.
Miyazaki
,
Rep. Prog. Phys.
75
,
036503
(
2012
).
20.
J.
Aarons
,
M.
Sarwar
,
D.
Thompsett
, and
C.-K.
Skylaris
,
J. Chem. Phys.
145
,
220901
(
2016
).
21.
S.
Goedecker
,
Phys. Rev. B
58
,
3501
(
1998
).
22.
S.
Ismail-Beigi
and
T. A.
Arias
,
Phys. Rev. Lett.
82
,
2127
(
1999
).
23.
X.
Zhang
and
D. A.
Drabold
,
Phys. Rev. B
63
,
233109
(
2001
).
24.
S. N.
Taraskin
,
P. A.
Fry
,
X.
Zhang
,
D. A.
Drabold
, and
S. R.
Elliott
,
Phys. Rev. B
66
,
233101
(
2002
).
25.
M.
Benzi
,
P.
Boito
, and
N.
Razouk
,
SIAM Rev.
55
,
3
(
2013
).
26.
M. J.
Gillan
,
D. R.
Bowler
,
A. S.
Torralba
, and
T.
Miyazaki
,
Comput. Phys. Commun.
177
,
14
(
2007
).
27.
C. K.
Skylaris
,
P. D.
Haynes
,
A. A.
Mostofi
, and
M. C.
Payne
,
J. Chem. Phys.
122
,
084119
(
2005
).
28.
E.
Tsuchida
,
J. Phys. Soc. Jpn.
76
,
034708
(
2007
).
29.
D.
Osei-Kuffuor
and
J.-L.
Fattebert
,
Phys. Rev. Lett.
112
,
046401
(
2014
).
30.
S.
Mohr
,
L. E.
Ratcliff
,
P.
Boulanger
,
L.
Genovese
,
D.
Caliste
,
T.
Deutsch
, and
S.
Goedecker
,
J. Chem. Phys.
140
,
204110
(
2014
).
31.
See www.openmx-square.org for OpenMX; accessed on 2020-06-05.
32.
N.
Bock
,
M.
Challacombe
,
C. K.
Gan
,
G.
Henkelman
,
K.
Nemeth
,
A. M. N.
Niklasson
,
A.
Odell
,
E.
Schwegler
,
C. J.
Tymczak
, and
V.
Weber
, FreeON,
Los Alamos National Laboratory (LA-CC 01-2; LA-CC-04-086), Copyright University of California
,
2014
.
33.
J.
Aarons
and
C.-K.
Skylaris
,
J. Chem. Phys.
148
,
074107
(
2018
).
34.
S.
Mohr
,
M.
Eixarch
,
M.
Amsler
,
M. J.
Mantsinen
, and
L.
Genovese
,
Nucl. Mater. Energy
15
,
64
(
2018
).
35.
A.
Ruiz-Serrano
,
N. D. M.
Hine
, and
C.-K.
Skylaris
,
J. Chem. Phys.
136
,
234101
(
2012
).
36.
Frontiers and Challenges in Warm Dense Matter
, Lecture Notes in Computational Science and Engineering, edited by
F.
Graziani
,
M. P.
Desjarlais
,
R.
Redmer
, and
S. B.
Tricky
(
Springer
,
2014
).
37.
F. R.
Graziani
,
V. S.
Batista
,
L. X.
Benedict
,
J. I.
Castor
,
H.
Chen
,
S. N.
Chen
,
C. A.
Fichtl
,
J. N.
Glosli
,
P. E.
Grabowski
,
A. T.
Graf
 et al,
High Energy Density Phys.
8
,
105
(
2012
).
38.
P.
Renaudin
,
C.
Blancard
,
J.
Clérouin
,
G.
Faussurier
,
P.
Noiret
, and
V.
Recoules
,
Phys. Rev. Lett.
91
,
075002
(
2003
).
39.
M.
Dharma-Wardana
,
Phys. Rev. E
73
,
036401
(
2006
).
40.
R.
Ernstorfer
,
M.
Harb
,
C. T.
Hebeisen
,
G.
Sciaini
,
T.
Dartigalongue
, and
R. J. D.
Miller
,
Science
323
,
1033
(
2009
).
41.
T. G.
White
,
S.
Richardson
,
B. J. B.
Crowley
,
L. K.
Pattison
,
J. W. O.
Harris
, and
G.
Gregori
,
Phys. Rev. Lett.
111
,
175002
(
2013
).
42.
P.
Suryanarayana
,
Chem. Phys. Lett.
584
,
182
(
2013
).
43.
P. P.
Pratapa
,
P.
Suryanarayana
, and
J. E.
Pask
,
Comput. Phys. Commun.
200
,
96
(
2016
).
44.
P.
Suryanarayana
,
P. P.
Pratapa
,
A.
Sharma
, and
J. E.
Pask
,
Comput. Phys. Commun.
224
,
288
(
2018
).
45.
P.
Suryanarayana
,
V.
Gavini
,
T.
Blesgen
,
K.
Bhattacharya
, and
M.
Ortiz
,
J. Mech. Phys. Solids
58
,
256
(
2010
).
46.
P.
Suryanarayana
,
K.
Bhattacharya
, and
M.
Ortiz
,
J. Comput. Phys.
230
,
5226
(
2011
).
47.
J. E.
Pask
and
P. A.
Sterne
,
Phys. Rev. B
71
,
113101
(
2005
).
48.
P.
Suryanarayana
and
D.
Phanish
,
J. Comput. Phys.
275
,
524
(
2014
).
49.
S.
Ghosh
and
P.
Suryanarayana
,
J. Comput. Phys.
307
,
634
(
2016
).
50.
R.
Martin
,
Electronic Structure: Basic Theory and Practical Methods
(
Cambridge University Press
,
2004
).
51.
S.
Goedecker
and
L.
Colombo
,
Phys. Rev. Lett.
73
,
122
(
1994
).
52.
A. M. N.
Niklasson
,
Phys. Rev. B
68
,
233104
(
2003
).
53.
A. M. N.
Niklasson
,
M. J.
Cawkwell
,
E. H.
Rubensson
, and
E.
Rudberg
,
Phys. Rev. E
92
,
063301
(
2015
).
54.
E.
Prodan
and
W.
Kohn
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
11635
(
2005
).
55.
J. E.
Pask
and
P. A.
Sterne
,
Modell. Simul. Mater. Sci. Eng.
13
,
R71
(
2005
).
56.
P.
Suryanarayana
,
K.
Bhattacharya
, and
M.
Ortiz
,
J. Mech. Phys. Solids
61
,
38
(
2013
).
57.
A. S.
Banerjee
,
P.
Suryanarayana
, and
J. E.
Pask
,
Chem. Phys. Lett.
647
,
31
(
2016
).
58.
P. P.
Pratapa
,
P.
Suryanarayana
, and
J. E.
Pask
,
J. Comput. Phys.
306
,
43
(
2016
).
59.
P.
Suryanarayana
,
P. P.
Pratapa
, and
J. E.
Pask
,
Comput. Phys. Commun.
234
,
278
(
2019
).
60.
S.
Kumar
,
Q.
Xu
, and
P.
Suryanarayana
,
Chem. Phys. Lett.
739
,
136983
(
2020
).
61.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
62.
N.
Troullier
and
J. L.
Martins
,
Phys. Rev. B
43
,
1993
(
1991
).
63.
D. R.
Hamann
,
Phys. Rev. B
88
,
085117
(
2013
).
64.
Q.
Xu
,
A.
Sharma
,
B.
Comer
,
H.
Huang
,
E.
Chow
,
A. J.
Medford
,
J. E.
Pask
, and
P.
Suryanarayana
, arXiv:2005.10431 (
2020
).
65.
S.
Ghosh
and
P.
Suryanarayana
,
Comput. Phys. Commun.
216
,
109
(
2017
).
66.
P.
Suryanarayana
,
Chem. Phys. Lett.
679
,
146
(
2017
).
67.
X.
Gonze
,
J.-M.
Beuken
,
R.
Caracas
,
F.
Detraux
,
M.
Fuchs
,
G.-M.
Rignanese
,
L.
Sindic
,
M.
Verstraete
,
G.
Zerah
,
F.
Jollet
 et al,
Comput. Phys. Commun.
25
,
478
(
2002
).
68.
P.
Minary
,
G. J.
Martyna
, and
M. E.
Tuckerman
,
J. Chem. Phys.
118
,
2510
(
2003
).
69.
J.
Hansen
and
I.
McDonald
,
Theory of Simple Liquids: With Applications to Soft Matter
(
Elsevier Science
,
2013
).
70.
D.
Alfè
and
M. J.
Gillan
,
Phys. Rev. Lett.
81
,
5161
(
1998
).
71.
T.
Sjostrom
and
J.
Daligault
,
Phys. Rev. E
92
,
063304
(
2015
).
72.
S.
Ghosh
and
P.
Suryanarayana
,
Comput. Phys. Commun.
212
,
189
(
2017
).
73.
A.
Sharma
and
P.
Suryanarayana
,
Chem. Phys. Lett.
700
,
156
(
2018
).
74.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
75.

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