# Fileset

[Mujahed_2026_J._Phys.__Condens._Matter_38_425501.pdf](https://mdr.nims.go.jp/filesets/c59d4b2f-7a66-426d-8120-4be004610b22/download)

## Creator

[Shereif Y Mujahed](https://orcid.org/0000-0002-7121-073X), [Tsuyoshi Miyazaki](https://orcid.org/0000-0003-3534-4404), [David R Bowler](https://orcid.org/0000-0001-7853-1520)

## Rights

[Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/)

## Other metadata

[Stress calculation in linear scaling DFT: convergence and dynamics](https://mdr.nims.go.jp/datasets/d5e60d4d-9699-40bc-b0c6-aad544e6ca55)

## Fulltext

Journal of Physics:Condensed Matter      PAPER • OPEN ACCESSStress calculation in linear scaling DFT:convergence and dynamicsTo cite this article: Shereif Y Mujahed et al 2026 J. Phys.: Condens. Matter 38 425501 View the article online for updates and enhancements.You may also likeIntroductory Remarks: Linear ScalingMethodsD R Bowler, J-L Fattebert, M J Gillan et al.-A localized-basis scheme for moleculardynamicsZ Lin and J Harris-Linear-scaling first-principles moleculardynamics of complex biological systemswith the Conquest codeTakao Otsuka, Makoto Taiji, David R.Bowler et al.-This content was downloaded from IP address 144.213.253.16 on 21/10/2026 at 02:28https://doi.org/10.1088/1361-648X/aeae97/article/10.1088/0953-8984/20/29/290301/article/10.1088/0953-8984/20/29/290301/article/10.1088/0953-8984/5/8/008/article/10.1088/0953-8984/5/8/008/article/10.7567/JJAP.55.1102B1/article/10.7567/JJAP.55.1102B1/article/10.7567/JJAP.55.1102B1/article/10.7567/JJAP.55.1102B1/article/10.7567/JJAP.55.1102B1https://pagead2.googlesyndication.com/pcs/click?xai=AKAOjsuGjp4Dnw50dnc0DJY3wG2bEd04hFs0DxLAz7aWi85ZvdIGkLs54t9Cy5-DtX3Tmk0B0gbeubdIL4TQn60vNLqqEfd7zZJOE8_oxuPWkkumjAz8y9OI6osdWLfxepb8sfqpg-Zf2EH7I-dSH6elF3x-NuaHxwdDKgyflBiCTkhrPfNA2YWsWZJT2ohevrCQa4lyO0xYSsIfd1JlwlBq04axulVHI_Mur8oic457JlI_QFYtYICYTOSukrsCJhY2vux2OqqNzK563d1tBHRDH_AQsk-LbKOTnFJP1u4ulKF-gO6EJBwkAdYdK2RQswC5rjPtoyW5eUzjZhVWY7az1q0FRxLtF2Dp283dP9sEHeEocI90wFth_EqNcNJkkXvWp661bDc&sig=Cg0ArKJSzNPnK2z0vCUg&fbs_aeid=%5Bgw_fbsaeid%5D&adurl=https://www.i-tech.si/quantum-lugus/%3Futm_source%3Dphysics-world%26utm_medium%3Dhalf-page%26utm_campaign%3Dquantum-lugusJ. Phys.: Condens. Matter 38 (2026) 425501 https://doi.org/10.1088/1361-648X/aeae97OPEN ACCESSRECEIVED8 July 2026REVISED24 September 2026ACCEPTED FOR PUBLICATION30 September 2026PUBLISHED20 October 2026Original content fromthis work may be usedunder the terms of theCreative CommonsAttribution 4.0 licence.Any further distributionof this work mustmaintain attribution tothe author(s) and the titleof the work, journalcitation and DOI.PAPERStress calculation in linear scaling DFT: convergence anddynamicsShereif Y Mujahed1,2, Tsuyoshi Miyazaki3 and David R Bowler1,2,3,∗1 Dept. Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, United Kingdom2 London Centre for Nanotechnology, University College London, 17–19 Gordon Street, London, WC1H 0AH, United Kingdom3 Research Centre for Materials Nanoarchitectonics (WPI-MANA), National Institute for Materials Science (NIMS), 1–1 Namiki,Tsukuba, Ibaraki 305-0044, Japan∗ Author to whom any correspondence should be addressed.E-mail: david.bowler@ucl.ac.ukKeywords: stress, linear scaling, molecular dynamicsAbstractWe present the approach needed to calculate stress within density functional theory (DFT) using alocalized orbital basis, both for exact diagonalisation and linear scaling approaches, and demon-strate our implementation within the large scale DFT code CONQUEST. For the linear scalingapproach, we test the rate of convergence of stress with density matrix range, and compare it tothe convergence of energy and forces for different materials with a range of band gaps. We showthat excellent convergence is found for modest cutoffs (below 0.1 GPa error for range of 20 a0), andshow that large-scale isothermal-isobaric molecular dynamics is stable and accurate, showing neg-ligible drift over picoseconds of dynamics.1. IntroductionDensity functional theory (DFT) has become the method of choice for electronic structure calculations[1]. In condensed matter research, plane wave basis sets have often been the basis of choice, which hastended to limit the size of all but the largest simulations to hundreds of atoms [2]; the adoption of localorbital basis sets has enabled larger scale simulations, both with exact diagonalisation and with linearscaling methods to find the ground state [3–7]. The CONQUEST code [7] has demonstrated simulations ofseveral thousand atoms with exact diagonalisation [8], and millions of atoms with linear scaling [9, 10].Simulations of condensed and soft matter systems require that the forces on the ions can be effi-ciently calculated as exact differentials of the energy; in many cases, the stresses on the simulation cellare also required. We have discussed the calculation of forces in CONQUEST for both exact diagonalisa-tion and linear scaling before [11], and noted that the motion of basis functions with the ions requiresevaluation of Pulay terms [12] which are readily found analytically. The evaluation of stresses is standardin plane-wave DFT [1] though for plane-wave basis set there is a missing term which arises from thechange of basis with simulation cell size (a Pulay stress), and cannot be calculated analytically, leading toa need to increase basis set size for stress evaluation. The calculation of stress for local orbital basis setshas been shown before [13, 14] but the calculation of stress with linear scaling, and the convergence ofthe stress, has only been investigated in the context of finite electronic temperature simulations [15]. Wereported preliminary results [7] on stress with linear scaling methods, but gave only an overview of theapproach; below, we give full details, and note an important correction to our preliminary results. Wealso demonstrate stable NPT molecular dynamics.The paper is organized as follows. We recap the expressions required to evaluate the energy both inexact diagonalisation and linear scaling modes; we note an important constraint that is applied to main-tain electron number in linear scaling calculations. We then examine the approach needed to calculatestresses with localized orbitals and the linear scaling approach in CONQUEST. We conclude by presentingvarious tests of the implementation; first, testing the convergence of energy, forces and stresses with den-sity matrix range, and then demonstrating stable molecular dynamics in different ensembles.© 2026 The Author(s). Published by IOP Publishing Ltdhttps://doi.org/10.1088/1361-648X/aeae97https://crossmark.crossref.org/dialog/?doi=10.1088/1361-648X/aeae97&domain=pdf&date_stamp=2026-10-20https://creativecommons.org/licenses/by/4.0/https://creativecommons.org/licenses/by/4.0/https://orcid.org/0000-0002-7121-073Xhttps://orcid.org/0000-0003-3534-4404https://orcid.org/0000-0001-7853-1520mailto:david.bowler@ucl.ac.ukJ. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et al2. EnergiesFor a system with total energy E we define the forces on the atoms and the stresses on the simulationcell in the usual way, with Fi =−∇iE for atom i, and σλµ = (1/Ω)∂E/∂ϵλµ, for the strain ϵλµ (with λand µ each representing one of x, y or z) and Ω the simulation cell volume. This definition of stress hasunits of pressure; an alternative definition, often known as the internal stress, omits the prefactor of 1/Ωto give units of energy [16, 17].In order to implement accurate forces and stresses, we have to consider the details of how energiesare calculated, and ensure that the forces and stresses are exact derivatives of the energy. In CONQUEST, incommon with other local orbital pseudopotential codes such as Siesta [14] and OpenMX [18], we use areformulation of the electrostatic part of the total energy based on neutral atom charges. By associatinga pseudo-atomic density (PAD) n(a)i (r) with each atom [19], we can rewrite the Hartree energy of theelectrons and the ionic core–core energy, avoiding any long range terms, leaving only a screened core-core interation (SCC). We combine the local part of the pseudopotential and the Hartree potential fromthe PAD to give a neutral atom potential: VNA(r) =∑iVPS,Loc,i(r)+V(a)Ha,i(r). We define the total PADfor the system as nPAD(r) =∑i n(a)i (r), and a charge density difference between the ground state chargedensity and the PAD as δn(r) = n(r)− nPAD(r). This density difference gives a potential VδHa and anenergy EδHa, and we define the Harris–Foulkes [20, 21] and Kohn–Sham (KS) [22] energies below inequations (2) and (4), with a Hamiltonian written as:Ĥ =(−12∇2 + V̂NA + V̂PS,NL +VδHa +VXC)(1)where V̂PS,NL is the non-local part of the pseudopotential.CONQUEST uses a basis set of localized functions, known as support functions [23], to represent theKS eigenstates. Most commonly these functions are represented by pseudo-atomic orbitals (PAOs) whichare generated by solving the Schrödinger equation for isolated atoms, using pseudopotentials [24]. Theeigenstates are then expressed as: |ψn⟩=∑iα ciαn |ϕiα⟩, where iα indicates the support functions on atomi, with the support functions indexed by α. Using these basis functions, the Schrödinger equation canbe re-written as a matrix equation, Hiαjβcjβn = ϵnSiαjβcjβn , noting that the support functions are non-orthogonal, with overlap matrix Siαjβ , and the matrices defined in the usual way: Hiαjβ = ⟨ϕiα|Ĥ|ϕjβ⟩and Siα,jβ = ⟨ϕiα|ϕjβ⟩. The individual terms in the energy can be defined either in terms of the eigen-values and eigenstates, or in terms of the Hamiltonian and density matrices; the latter approach, coupledwith sparse matrices, enables a linear scaling formulation of DFT to be constructed [25, 26].2.1. Energy in terms of eigenvalues and eigenstatesWe present first the Harris–Foulkes (HF) and KS energies used when solving by diagonalisation; belowwe turn to the linear scaling formulation. We also consider the inclusion of non-linear partial core cor-rections (PCC) [27] which include an additional charge density in the core region of each atom, nCi (r),in the calculation of the exchange-correlation energy, and for simplicity we define ñ(r) = n(r)+ nC(r),where nC(r) =∑i nCi (r).EHF = EBand +∆EHa +∆EXC + ESCC (2)EHF =∑n2fnϵn +[−12ˆdrδVHa (r)δn(r)−ˆdrδVHa (r)nPAD (r)]+ˆdr(ñ(r)ϵXC [ñ(r)]− n(r)µXC [ñ(r)])+12∑i,jZiZj∣∣Ri −Rj∣∣ −ˆdrV(a)Ha,j (r)n(a)i (r) (3)EKS = EKE + ENA + EPS,NL + EδHa + EXC + ESCC (4)2J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alEKS =−12∑n2fn⟨ψn∣∣∇2∣∣ψn⟩+ˆdrVNA (r)n(r)+∑n2fn⟨ψn|V̂PS,NL|ψn⟩+12ˆdrδn(r)δVHa (r)+ˆdrñ(r)ϵXC [ñ(r)]+12∑i,jZiZj∣∣Ri −Rj∣∣ −ˆdrV(a)Ha,j (r)n(a)i (r) (5)where ϵn are the eigenvalues indexed with subscript n, fn are the occupancies of the eigenstates and thenon-local pseudopotential is defined as V̂PS,NL =∑i∑γ |χiγ⟩Aγ⟨χiγ |, summing over atoms i and projec-tors γ in the usual way [28].When operating without self-consistency [29], we use the fixed superposition of atomic densitiesdefined above as the density, so that n(r) = nPAD(r), which fits well with the neutral-atom approach. Inthis case, only the HF energy is reliable, and ∆EHa = 0 as δn(r) = 0, while ∆EXC is found as usual, butwith ñ(r) = nPAD(r)+ nC(r).2.2. Energy in terms of the density matrixFor linear scaling operation, all calculations are performed in terms of the density matrix and thecharge density [9, 30]. The density matrix in real space is defined in terms of the support functionsand a matrix Kiαjβ as ρ(r,r ′) =∑iαjβ ϕiα(r)Kiαjβϕjβ(r ′), and the connection to eigenvector coeffi-cients can be seen in the formal definition of Kiαjβ =∑n fnciα⋆n cjβn . The charge density is then writtenas n(r) =∑iαjβ ϕiα(r)Kiαjβϕjβ(r) and the energies are given by:EHF =∑n2Tr [KH] +[−12ˆdrδVHa (r)δn(r)−ˆdrδVHa (r)nPAD (r)]+ˆdr(ñ(r)ϵXC [ñ(r)]− n(r)µXC [ñ(r)])+12∑i,jZiZj∣∣Ri −Rj∣∣ −ˆdrV(a)Ha,j (r)n(a)i (r) (6)EKS =−12∑iαjβ2Kjβiα⟨ϕiα∣∣∇2∣∣ϕjβ⟩+ˆ drVNA (r)n(r)+∑iαjβ2Kjβiα⟨ϕiα|V̂PS,NL|ϕjβ⟩+12ˆdrδn(r)δVHa (r)+ˆdrñ(r)ϵXC [ñ(r)]+12∑i,jZiZj∣∣Ri −Rj∣∣ −ˆdrV(a)Ha,j (r)n(a)i (r) (7)For relatively small systems (up to a few thousand atoms) the ground state is found by direct diago-nalisation of the Hamiltonian [7]. To go beyond systems of this size, linear scaling methods are needed[25], as described in the next section.2.3. Linear scaling and electron numberThe ground state density matrix of semiconducting and insulating systems (i.e. systems with a gap) canbe found with computational effort that scales linearly with system size [25, 26] by working with local-ized functions as basis functions, and restricting the range of the density matrix (a controlled approxi-mation). Within CONQUEST, a variational minimization of the energy with respect to the density matrixelements is performed, while maintaining two properties of the system: idempotency of the densitymatrix; and correct total electron number.Idempotency is imposed approximately, using the LNV approach [31, 32], where the density matrixK is written in terms of an auxiliary density matrix, L, as K= 3LSL− 2LSLSL and the elements of theauxiliary density matrix, Liαjβ , are varied in the search for the ground state. The initial density matrix isconstructed using the Palser-Manolopoulos generalization of the McWeeny transform [33, 34].The electron number is found as N(K)e = 2Tr[KS], where the superscript (K) is used to distinguish thecalculated electron number from the desired electron number, Ne. To optimize the energy with respectto the elements of the auxiliary density matrix, Liαjβ , while maintaining electron number, we use aLagrange multiplier, µ, to impose the condition [23]. As above, we have Eband = 2Tr [KH], and we definea modified band energy which includes the Lagrange multiplier:Ẽband = Eband −µ(N(K)e −Ne)(8)= 2Tr [KH]−µ(2Tr [KS]−Ne) (9)3J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alNaturally, when the density matrix gives the target electron number, this energy is identically equal tothe band energy. It is sometimes convenient to define H ′ =H−µS and Ẽband = 2Tr [KH ′] +µNe toinclude the Lagrange multiplier in the Hamiltonian [23], but this can obscure a contribution to theforces discussed below, and we do not follow this convention.To find the optimum ground state, we need the derivative of this energy with respect to the auxiliarydensity matrix elements:∂Ẽband∂L=∂Eband∂L−µ∂N(K)e∂L∂Ẽband∂L= [6(HLS+ SLH)− 4(HLSLS− SLHLS− SLSLH)]−µ [12(SLS− SLSLS)] . (10)Here we have defined the change in electron number calculated from the K matrix with respect tothe elements of the auxiliary density matrix as ∂N(K)e /∂L. We can then set the value of the Lagrangemultiplier so that, to first order, the electron number is preserved while the energy is minimized. Thisis imposed by ensuring that the variation in energy is orthogonal to the variation in electron number:∂Ẽband∂L· ∂N(K)e∂L= 0(∂Eband∂L−µ∂N(K)e∂L)· ∂N(K)e∂L= 0⇒(∂Eband∂L· ∂N(K)e∂L)/(∂N(K)e∂L· ∂N(K)e∂L)= µ. (11)Of course, the Lagrange multiplier, µ, is the chemical potential for the electrons.At the electronic ground state, where ∂Ẽband/∂L= 0, we can find the contribution to the ionic forcesand the stresses from the band energy (and associated Lagrange multiplier term) by considering varia-tions in the energy with respect to the Hamiltonian and overlap matrices [11]:Ẽband = Eband −µ(N(K)e −Ne)= 2Tr [KH]−µ(2Tr [KS]−Ne) (12)δẼband = δEband −µδN(K)e (13)δEband = 2Tr [(3LδSL− 2LδSLSL− 2LSLδSL)H] + 2Tr [(3LSL− 2LSLSL)δH]= 2Tr [KδH]− 2Tr [GδS] (14)G= 2LHLSL+ 2LSLHL− 3LHL (15)δN(K)e = 2Tr [(3LδSL− 2LδSLSL− 2LSLδSL)S] + 2Tr [(3LSL− 2LSLSL)δS]= 2Tr [(6LSL− 6LSLSL)δS] (16)δẼband = 2Tr [KδH]− 2Tr [GδS]− 2µTr [(6LSL− 6LSLSL)δS] = 2Tr [KδH]− 2Tr[G̃δS](17)G̃= 2LHLSL+ 2LSLHL− 3LHL+µ(6LSL− 6LSLSL) . (18)The sign of the matrix G is chosen for compatibility with the equivalent matrix required for force andstress calculations with exact diagonalisation:G iαjβ =∑nfnϵnciα⋆n cjβn . (19)It is important to note that in O(N) calculations there is a contribution to the forces and stressesfrom the variation of the electron number with the overlap matrix which arises from the Lagrange mul-tiplier: the final term in equations (13) and (17). This term makes an almost negligible change to theforces, typically around 1× 10−4 Ha/a0, and was neglected in our original formulation of forces [11];the effect of the term is only visible as a small variation in the energy conservation of NVE moleculardynamics runs [35], and makes no difference to optimized structures. However, the term forms a signif-icant part of the stresses, and must be included for correct agreement between numerical and analyticstresses at small L matrix ranges. As the L matrix range increases, the contribution tends to zero, sincethe L matrix tends towards perfect idempotency, with LSL= L; the small fluctuations in energy overtime in NVE MD runs also reduce with increasing density matrix range, and are eliminated when thisextra contribution to the force is included.4J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et al3. Formulating stress with localized orbitals and linear scalingWe have separated the contributions to the stress into different types, in a similar way to the presenta-tion of forces [11]. Hellmann–Feynman stresses arise from explicit dependence of the potentials on theatom positions via the strain, while Pulay stresses come from the change in the basis when atoms move(since the basis is atom-centered). There are then further miscellaneous stresses which are important butdo not fit neatly into any category. We note that the calculation of the stress requires negligible time, asit is included with the force evaluation (this is standard for local orbital codes, in contrast to plane wavebasis sets).3.1. Hellmann–Feynman stressThe Hellmann–Feynman stress is straightforward, and equivalent to other codes, coming from the pseu-dopotentials; there is an extra contribution from the screened Hartree potential in the neutral atom for-mulation.σHFλµ =∑iˆdrn(r)∂VNAi (r)∂RiλRiµ − 2∑i∑abKab⟨ϕa | (∇iλχb)Riµ⟩Ai⟨χi | ϕb⟩−∑iˆdrδVHa (r)∂nPADi∂RiλRiµ. (20)To understand why we classify the final term as a Hellmann–Feynman term, we note that the screenedHartree energy is written as EδHa =12´drδVHa(n(r)−∑i nPADi (r)), so there is a contribution to theforce and stress from the change of nPADi (r) with atomic position. The equivalence to a Hellmann–Feynman force comes because the integral can be recast as δn(r) interacting with VPADHa (r), the potentialdue to the atomic density.3.2. Pulay stressesFollowing the forces [11], we find that there are two contributions to the Pulay stress, coming fromchanges in the Hamiltonian and overlap matrices; we term these ϕ-Pulay and S-Pulay stresses, respec-tively. In equation (21) below, the first three terms form the ϕ-Pulay stress (contributions from the localpotentials, i.e. neutral atom, Hartree and XC potentials, in the first term; the kinetic energy in the sec-ond term; and the non-local potential in the third term) while the fourth contribution forms the S-Pulaystress.σPulayλµ =−4∑i∑abKab⟨ϕa | V̂NA + δV̂Ha + V̂XC | (∇iλϕb)Riµ⟩+12∑i∑abKab∂Tab∂RiλRiµ − 2∑i∑abkKab⟨(∇iλϕa) | χk⟩RiµAk⟨χk | ϕb⟩− 4∑i∑abG̃ab⟨ϕa | (∇iλϕb)⟩Riµ. (21)It is important to note that, for linear scaling, the matrix G̃ab is the same matrix found in equation (18),including the effect of the variation of the electron number, while for exact diagonalisation it issimply Gab.3.3. Miscellaneous stressesThe remaining stresses do not fit easily into the two previous categories, and are given here in no partic-ular order. The ion-ion stress within the neutral atom potential is written as follows:σSCCλµ =∑i12∑j[ZiZj| Ri −Rj |3(Riλ −Rjλ)δλµ −ˆdrVHa,PADj (r)∂nPADi∂Riλ](Riµ −Rjµ). (22)Various terms in the total energy are calculated using a numerical integration on a real-space grid(the screened Hartree, exchange-correlation and neutral atom energies). There is a stress associated withthe change of the Jacobian (i.e. the grid point volume):σJacobianλµ = δλµ(EδHa + EXC + ENA). (23)5J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alSimilarly, since the Hartree potential is found using fast Fourier transforms, there is a stress which arisesfrom changes of the reciprocal lattice vectors:σHaλµ =ˆdrδn(r)FT[2GλGµG 4δ̃n](24)where δ̃n indicates the Fourier transform of δn, and FT[] indicates the inverse Fourier transform.For a GGA exchange-correlation functional, which depends on the gradient of the electron density,there is an extra contribution to the stress which can be written as [36, 37]:σXC,GGAλµ =ˆdr∂f∂∇λn(r)∇µn(r) . (25)When using non-linear PCC, there is a stress associated with the change of the core densities asatoms move, which can be written as:σPCCλµ =∑iˆdr∇iλnCi (r)VXC[n(r)+ nC (r)]Riµ. (26)It might be argued that this is similar to a Hellmann–Feynman stress, but the distinction is notimportant.Finally, when running without self-consistency, there are stresses equivalent to the forces we havedescribed elsewhere [29]:σNon−SCFλµ =−ˆdr∆VHa (r)∑i∂RiλnPADi Riµ−ˆdr∆n(r)µ ′XC[nPAD + nC]∑i∂Riλ(nPADi + nCi)Riµ. (27)There is also an extra Jacobian stress for the exchange and correlation:σXCJacobian,Non−SCF =ˆdr(nout − nPAD)VXC[nPAD + nC]. (28)The Hartree stress arising from the change in reciprocal lattice vectors is identically zero for non-self-consistent calculations.4. Tests and resultsIn this section we present comprehensive tests of the convergence of pressure, force and total energywith density matrix range for a variety of materials with differing properties, and then test moleculardynamics in the NPT ensemble. In all cases we compare results found with linear scaling to those foundwith converged exact diagonalisation. The convergence of total energy with respect to density matrixrange is related to the gap of the material [38]; while there are no formal results on the convergence offorces or stresses with density matrix range, energy differences converge faster than total energies, andmost linear scaling calculations rely on the fast convergence of energy differences, even when total ener-gies are not converged. It is important to note that the same is true for many parameters in standardDFT calculations, e.g. k-point convergence, or plane wave cutoff in plane wave codes. For simplicity, weuse a simple single zeta basis set of PAOs [24] for all calculations; the only significant effect that a largerbasis set will have on the results will be if the gap is changed by the basis. Any basis set superpositionerrors will have a negligible effect on stress calculations in condensed phases.All tests used Conquest v1.5 [39] which is publicly available from GitHub, and includes the stressimplementation described in this paper.4.1. Convergence testsWe consider six materials which provide a variety of types of bonding and electronic structure. Theelemental semiconductors carbon, silicon and germanium are all tetrahedrally bonded, and have vary-ing gaps (the values with the single zeta basis are not accurate, but nevertheless span an appropriaterange: 5.9 eV for carbon; 2.1 eV for silicon; and 0.0 eV for germanium at the Γ point), covering behav-iors from insulating, through semiconducting to metallic. We also include ionic materials, all of which6J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alFigure 1. Absolute value of the difference between linear scaling and fully converged exact diagonalisation calculations as a func-tion of density matrix cutoff (in a0) for: (a) stress (in GPa); (b) force (in Ha/a0); and (c) total energy (in Ha). (d) stress differencefor Si and Ge, as in (a), on a linear scale to illustrate the cause of the dip in the stress difference curve.have large gaps: MgO and LiCl, which are both cubic, with gaps of 5.3 eV and 7.0 eV; and AlN, whichis hexagonal (wurtzite), and polar along the c-axis, with a gap of 3.4 eV. This variety of materials cov-ers differences in bonding and electronic structure, and provides a stringent test of the approach. ForMgO and LiCl we use the local density approximation with Troullier-Martins pseudopotentials [40],while for the other materials we use PBE [41] with Hamann pseudopotentials [28] from the PseudoDojodatabase [42].In our comparisons, we use the same simulation cells, with all parameters identical except for theground state search. In all cases, we used a perfect crystal but with one atom slightly displaced along thec-axis of the cell to give non-zero forces, with a lattice constant reduced by 1% from the relaxed value(for exact diagonalisation). We use a real space integration grid with a spacing of 0.2 a0, which givesgood convergence. All diagonalisation calculations used a Monkhorst–Pack grid which fully convergedthe total energy.We consider three quantities: the total energy; the force; and the stress along the c-axis. In all practi-cal circumstances we are interested in energy differences rather than total energies, so the stress and forcegraphs are much more important. The results are shown in figure 1, showing the stress in (a), the forcein (b) and the total energy in (c). The unit cells of all materials apart from AlN are cubic; we have plot-ted the stress along the c axis in AlN, which is the direction of the displacement of the ion, though theconvergence is similar in all directions.Considering first the stress, we see that the discrepancy is rather small, being below 1 GPa for allmaterials at all density matrix ranges (except carbon with the exceptionally small density matrix range of12 a0) and below 0.1GPa for very reasonable density matrix ranges: 18 a0 for the four insulators (C, LiCl,MgO and AlN), 20 a0 for silicon, and 22 a0 for germanium. The dip in the curve seen around 25 a0 forSi and Ge is due to the non-monotonic convergence of the stress with density matrix range, as shown infigure 1(d); there is no fundamental reason why the stress should converge monotonically.The force convergence is even better than the stress, with the discrepancies below 1×10−4 Ha/a0(≃ 5× 10−3 eV Å−1) for all density matrix ranges and all materials (except AlN along a direction inthe basal plane where this level requires a range of 18 a0). This indicates that structural optimization forall elements considered here will be reliable regardless of the restriction placed on the density matrix.The final curve, showing convergence of total energy, displays the expected trends of exponential conver-gence, with larger gaps converging faster. Absolute energies are much less important than energy differ-ences, which converge much faster (as seen in the stresses and forces). We note that comparison of the7J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alTable 1. Bulk modulus in GPa calculated by fitting Birch–Murnaghan curves toenergy (B(E)0 ) and pressure (B(P)0 ) for the six materials, for linear scaling(RL = 18a0 for C, LiCl, MgO, AlN, RL = 20a0 for Si and RL = 22a0 for Ge) andexact diagonalisation.C Si Ge LiCl MgO AlNO(N) B(E)0 391.8 72.5 46.0 39.8 187.3 182.1Diag B(E)0 391.8 71.3 46.9 39.6 185.5 181.7O(N) B(P)0 395.6 71.9 46.2 39.3 185.6 180.9Diag B(P)0 395.5 70.5 47.9 39.1 185.4 180.4Figure 2. Characterization of the effect of density matrix range on NPT molecular dynamics. (a) Deviation of conserved quan-tity; (b) temperature and (c) pressure variation. In all cases results for full diagonalisation are compared to three density matrixranges (12, 16 and 20 a0).analytic stresses presented here with numerical stresses found from the energy change resulting from asmall volume change are in excellent agreement.As a further test of the stress calculation, we evaluate the bulk modulus for all materials by fittingto seven point curves in both the energy (as is standard for DFT calculations) and the pressure on thesimulation cell (which tests the stress explicitly). For the linear scaling calculations, we use the same den-sity matrix range identified above where the stress error is less than 0.1 GPa (18a0 for C, LiCl, MgO andAlN; 20 a0 for Si and 22 a0 for Ge). We found that rather than using a fixed density matrix range foreach point on the binding energy curve, better convergence is obtained by scaling the density matrixrange with simulation cell size [43]. The results are given in table 1, and show excellent agreement.4.2. Molecular dynamics testsWe now turn to a practical application of the calculated stress: isobaric-isothermal molecular dynam-ics. We have already demonstrated NVE molecular dynamics with linear scaling DFT on systems up to32 678 atoms [35], using the extended Lagrangian Born–Oppenheimer (XL-BOMD) approach [44]. Inthat work, we showed that the total energy was well conserved, with small fluctuations (no more than1 meV atom−1 even for a small density matrix range of 14 a0) which reduced with increasing densitymatrix range. We have also demonstrated dynamics under the NVT ensemble [45], again showing stabledynamics with good conservation. Here, we show that stable linear scaling molecular dynamics is pos-sible for the NPT ensemble, and characterize the behavior with density matrix range. We choose a tem-perature of 300 K and a pressure of 10 GPa for a sample of bulk silicon. We are particularly interested inthe effects of the density matrix range on fidelity of the dynamics: both the stability and conservation ofthe simulations, and how closely the full diagonalisation dynamics can be followed.In figure 2 we consider the case where we used identical simulation settings for all tests, includ-ing the same starting seed for the pseudo-random number generator and starting simulation cell size.We note that the initial pressures in these cases differ with density matrix range (as seen in figure 1),with a difference to exact diagonalisation of less than 0.1 GPa at RL = 20a0. As the density matrix rangeincreases, the linear scaling dynamics converge towards the exact diagonalisation dynamics; in all cases,we see that the conserved quantity does not drift, and the pressure and temperature are maintained attheir set point.As would be expected, the smaller density matrix ranges show more differences in the detaileddynamics to the exact diagonalisation results, while the RL = 20 a0 results follow the exact diagonali-sation results closely. In particular, the pressure fluctuations are out of phase with the exact diagonal-isation results, and the conserved quantity follows a somewhat different path (though with the sameunderlying frequency of oscillation). The key question to understand in terms of fidelity of dynamicsis whether this is purely due to the difference in initial pressure, or also due to differences in forces;8J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alFigure 3. Comparison of NPT molecular dynamics starting from the same initial pressure for diagonalisation and linear scalingcalculations. (a) Deviation of conserved quantity; (b) temperature and (c) pressure variation. In all cases, results for full diago-nalisation are compared to two density matrix ranges (12 and 16 a0).from the results in figure 1(b) we would expect it to be almost entirely from pressure. We investigateddynamics for short density matrix ranges (RL = 12 and 16 a0) with simulation cell sizes adjusted to givethe same initial pressure as in the exact diagonalisation simulation (to within 0.1GPa). The results areshown in figure 3.We see that for the first hundred femtoseconds of dynamics the short density matrix range simu-lations follow the exact diagonalisation trajectories extremely closely. The increase in difference overtime is an inevitable result of slight variations in dynamics. The result of these simulations is extremelyencouraging for large scale, long-time DFT dynamics: we see that high quality simulations can beachieved using very modest density matrix ranges. Comparing diagonalisation and O(N) calculations[46] we find a drift of less than 0.001 eV ps−1 atom−1 for simulations spanning multiple picoseconds.5. ConclusionsWe have shown how stress can be calculated in the context of the local orbital DFT code CONQUEST, bothin terms of exact diagonalisation and linear scaling operation. Tests of the convergence of the stress cal-culated with linear scaling to the exact diagonalisation result show that relatively modest density matrixranges are required for excellent agreement, with all materials converged to within 0.1 GPa with a rangeof 20 a0. Tests of NPT dynamics show no drift in the conserved quantity, and stable dynamics overmeaningful timescales. Moreover, adjusting the starting simulation cell size so that the initial pressuresof both linear scaling and exact diagonalisation are the same allows us to reproduce the exact dynamicsover 100 fs for very small ranges, as low as 12 a0 for silicon.Our results demonstrate that the full range of molecular dynamics ensembles can be used with lin-ear scaling DFT methods with no loss of accuracy in the dynamics of the system, even with modestdensity matrix ranges. This opens the prospect of very large scale dynamics calculations with full DFTaccuracy. The remaining questions for linear scaling MD simulations concern basis set size and over-all efficiency. Since the support functions in CONQUEST are non-orthogonal, we require the inverse ofthe overlap matrix for efficient calculations [32, 47], but with linear scaling time. At present, we findan approximate, short range inverse matrix using Hotelling’s method, though this approach performspoorly for larger PAO basis sets as the overlap matrix becomes increasingly ill-conditioned. The lackof an inverse overlap matrix potentially restricts the basis size that can be used, though the use of thealternative blip basis sets [23, 30] and the on-site support function approach outlined recently [7] offerways forward here. The question of efficiency is key to long timescale MD runs. We already use theextended Lagrangian approach to molecular dynamics [44] to ensure energy conservation, though thiscan require more iterations to find the ground state than simple reuse of the density matrix; exten-sions of the extended Lagrangian approach that require no self-consistency cycles [48] may help withefficiency. Once these questions have been addressed, the route to very large scale, long DFT moleculardynamics runs is open.AcknowledgmentsThis research was partly carried out using the high performance computing cluster in the LondonCenter for Nanotechnology. The authors acknowledge the use of the UCL Kathleen High Performance9J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et alComputing Facility (KathleenUCL), and associated support services, in the completion of this work. Weare grateful for computational support from the UK Materials and Molecular Modeling Hub, which ispartially funded by EPSRC (EP/T022213/1, EP/W032260/1 and EP/P020194/1), for which access wasobtained via the UKCP consortium and funded by EPSRC grant ref EP/X035891/1.Data availability statementThe data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.5281/zenodo.21262374 [49].Author contributionsShereif Y Mujahed  0000-0002-7121-073XFormal analysis (equal), Methodology (equal), Software (equal), Writing – original draft (equal)Tsuyoshi Miyazaki  0000-0003-3534-4404Conceptualization (equal), Supervision (equal), Writing – review & editing (equal)David R Bowler  0000-0001-7853-1520Conceptualization (equal), Formal analysis (equal), Investigation (equal), Methodology (equal),Software (equal), Supervision (equal), Writing – original draft (equal), Writing – review &editing (equal)References[1] Martin R M 2020 Electronic Structure: Basic Theory and Practical Methods 2nd edn (Cambridge University Press) (https://doi.org/10.1017/9781108555586)[2] Carnimeo I et al 2023 Quantum espresso: One further step toward the exascale J. Chem. Theory Comput. 19 6992–7006[3] Ratcliff L E et al 2020 Flexibilities of wavelets as a computational basis set for large-scale electronic structure calculations J. Chem.Phys. 152 194110[4] García A et al 2020 Siesta: Recent developments and applications J. Chem. Phys. 152 204108[5] Prentice J C A et al 2020 The onetep linear-scaling density functional theory program J. Chem. Phys. 152 174111[6] Kühne T D et al 2020 Cp2k: an electronic structure and molecular dynamics software package - quickstep: Efficient and accurateelectronic structure calculations J. Chem. Phys. 152 194103[7] Nakata A et al 2020 Large scale and linear scaling DFT with the CONQUEST code J. Chem. Phys. 152 164112[8] Nakata A, Bowler D R and Miyazaki T 2022 Large-scale DFT methods for calculations of materials with complex structuresJ. Phys. Soc. Jpn. 91 091011[9] Bowler D R and Miyazaki T 2010 Calculations on millions of atoms with density functional theory: linear scaling shows its poten-tial J. Phys.: Condens. Matter 22 074207[10] Arita M, Arapan S, Bowler D R and Miyazaki T 2014 Large-scale DFT simulations with a linear-scaling DFT code CONQUEST onK-computer J. Adv. Simul. Sci. Eng. 1 87–97[11] Miyazaki T, Bowler D R, Choudhury R and Gillan M J 2004 Atomic force algorithms in density functional theory electronic-structure techniques based on local orbitals J. Chem. Phys. 121 6186–94[12] Pulay P 1969Mol. Phys. 17 197[13] Feibelman P J 1991 Pulay-type formula for surface stress in a local-density-functional, linear combination of atomic orbitals,electronic-structure calculation Phys. Rev. B 44 3916–25[14] Soler J M, Artacho E, Gale J D, García A, Junquera J, Ordejón P and Sánchez-Portal D 2002 The SIESTA method for ab initioorder-N materials simulation J. Phys.: Condens. Matter 14 2745–70[15] Sharma A, Hamel S, Bethkenhagen M, Pask J E and Suryanarayana P 2020 Real-space formulation of the stress tensor for o(n)density functional theory: application to high temperature calculations J. Chem. Phys. 153 034112[16] Nielsen O H and Martin R M 1983 First-principles calculation of stress Phys. Rev. Lett. 50 697–700[17] Nielsen O H and Martin R M 1985 Quantum-mechanical theory of stress and force Phys. Rev. B 32 3780–91[18] Ozaki T and Kino H 2005 Efficient projector expansion for the ab initio lcao method Phys. Rev. B 72 045121[19] We have a pseudo-atomic density, since we use a pseudopotential and solve for pseudo-atomic orbitals (PAOs)[20] Harris J 1985 Phys. Rev. B 31 1770[21] Foulkes W and Haydock R 1989 Phys. Rev. B 39 12520[22] Kohn W and Sham L J 1965 Phys. Rev. 140 A1133[23] Hernández E, Gillan M J and Goringe C M 1996 Linear-scaling density-functional-theory technique: The density-matrixapproach Phys. Rev. B 53 7147–57[24] Bowler D R, Baker J S, Poulton J T L, Mujahed S Y, Lin J, Yadav S, Raza Z and Miyazaki T 2019 Highly accurate local basis sets forlarge-scale DFT calculations in conquest Jpn. J. Appl. Phys. 58 100503[25] Bowler D R and Miyazaki T 2012O(N) methods in electronic structure calculations Rep. Prog. Phys. 75 036503[26] Goedecker S 1999 Linear scaling electronic structure methods Rev. Mod. Phys. 71 1085[27] Louie S G, Froyen S and Cohen M L 1982 Phys. Rev. B 26 1738[28] Hamann D R 2013 Optimized norm-conserving vanderbilt pseudopotentials Phys. Rev. B 88 08511710https://doi.org/10.5281/zenodo.21262374https://doi.org/10.5281/zenodo.21262374https://orcid.org/0000-0002-7121-073Xhttps://orcid.org/0000-0002-7121-073Xhttps://orcid.org/0000-0003-3534-4404https://orcid.org/0000-0003-3534-4404https://orcid.org/0000-0001-7853-1520https://orcid.org/0000-0001-7853-1520https://doi.org/10.1017/9781108555586https://doi.org/10.1017/9781108555586https://doi.org/10.1021/acs.jctc.3c00249https://doi.org/10.1021/acs.jctc.3c00249https://doi.org/10.1063/5.0004792https://doi.org/10.1063/5.0004792https://doi.org/10.1063/5.0005077https://doi.org/10.1063/5.0005077https://doi.org/10.1063/5.0004445https://doi.org/10.1063/5.0004445https://doi.org/10.1063/5.0007045https://doi.org/10.1063/5.0007045https://doi.org/10.1063/5.0005074https://doi.org/10.1063/5.0005074https://doi.org/10.7566/JPSJ.91.091011https://doi.org/10.7566/JPSJ.91.091011https://doi.org/10.1088/0953-8984/22/7/074207https://doi.org/10.1088/0953-8984/22/7/074207https://doi.org/10.15748/jasse.1.87https://doi.org/10.15748/jasse.1.87https://doi.org/10.1063/1.1787832https://doi.org/10.1063/1.1787832https://doi.org/10.1080/00268976900100941https://doi.org/10.1080/00268976900100941https://doi.org/10.1103/PhysRevB.44.3916https://doi.org/10.1103/PhysRevB.44.3916https://doi.org/10.1088/0953-8984/14/11/302https://doi.org/10.1088/0953-8984/14/11/302https://doi.org/10.1063/5.0016783https://doi.org/10.1063/5.0016783https://doi.org/10.1103/PhysRevLett.50.697https://doi.org/10.1103/PhysRevLett.50.697https://doi.org/10.1103/PhysRevB.32.3780https://doi.org/10.1103/PhysRevB.32.3780https://doi.org/10.1103/PhysRevB.72.045121https://doi.org/10.1103/PhysRevB.72.045121https://doi.org/10.1103/PhysRevB.31.1770https://doi.org/10.1103/PhysRevB.31.1770https://doi.org/10.1103/PhysRevB.39.12520https://doi.org/10.1103/PhysRevB.39.12520https://doi.org/10.1103/PhysRev.140.A1133https://doi.org/10.1103/PhysRev.140.A1133https://doi.org/10.1103/PhysRevB.53.7147https://doi.org/10.1103/PhysRevB.53.7147https://doi.org/10.7567/1347-4065/ab45afhttps://doi.org/10.7567/1347-4065/ab45afhttps://doi.org/10.1088/0034-4885/75/3/036503https://doi.org/10.1088/0034-4885/75/3/036503https://doi.org/10.1103/RevModPhys.71.1085https://doi.org/10.1103/RevModPhys.71.1085https://doi.org/10.1103/PhysRevB.26.1738https://doi.org/10.1103/PhysRevB.26.1738https://doi.org/10.1103/PhysRevB.88.085117https://doi.org/10.1103/PhysRevB.88.085117J. Phys.: Condens. Matter 38 (2026) 425501 S Y Mujahed et al[29] Torralba A S, Bowler D R, Miyazaki T and Gillan M J 2009 Non-self-consistent Density-Functional Theory Exchange-CorrelationForces for GGA Functionals J. Chem. Theory Comput. 5 1499[30] Bowler D R, Miyazaki T and Gillan M J 2002 Recent progress in linear scaling ab initio electronic structure techniques J. Phys.:Condens. Matter 14 2781–98[31] Li X-P, Nunes R W and Vanderbilt D 1993 Density-matrix electronic-structure method with linear system-size scaling Phys. Rev. B47 10891–4[32] Nunes R W and Vanderbilt D 1994 Generalization of the density-matrix method to a nonorthogonal basis Phys. Rev. B 50 17611–4[33] Bowler D R and Gillan M J 1999 Density matrices in O(N) electronic structure calculations: theory and applications Comput.Phys. Commun. 120 95–108[34] Palser A H R and Manolopoulos D E 1998 Canonical purification of the density matrix in electronic-structure theory Phys. Rev. B58 12704–11[35] Arita M, Bowler D R and Miyazaki T 2014 Stable and efficient linear scaling first-principles molecular dynamics for 10,000+atoms J. Chem. Theory Comput. 10 5419[36] Corso A D and Resta R 1994 Density-functional theory of macroscopic stress: Gradient-corrected calculations for crystalline sePhys. Rev. B 50 4327–31[37] Balbás L C, Martins J L and Soler J M 2001 Evaluation of exchange-correlation energy, potential and stress Phys. Rev. B 64 165110[38] Prodan E and Kohn W 2005 Nearsightedness of electronic matter Proc. Natl Acad. Sci. USA 102 11635[39] Conquest v1.5 2026 (available at: https://github.com/OrderN/CONQUEST-release) (24 September 2026)[40] Troullier N and Martins J L 1991 Efficient pseudopotentials for plane-wave calculations Phys. Rev. B 43 1993–2006[41] Perdew J P, Burke K and Ernzerhof M 1996 Generalized gradient approximation made simple Phys. Rev. Lett. 77 3865–8[42] van Setten M J, Giantomassi M, Bousquet E, Verstraete M J, Hamann D R, Gonze X and Rignanese G M 2018 The pseudodojo:training and grading a 85 element optimized norm-conserving pseudopotential table Comput. Phys. Commun. 226 39–54[43] This ensures that the number of non-zero density matrix elements is constant across all calculations, in close analogy to similarcalculations performed with a plane wave basis set, which generally use a constant number of plane waves[44] Niklasson A M N 2008 Phys. Rev. Lett. 100 123004[45] Hirakawa T, Suzuki T, Bowler D R and Miyazaki T 2017 Canonical-ensemble extended Lagrangian Born–Oppenheimer moleculardynamics for the linear scaling density functional theory J. Phys.: Condens. Matter 29 405901[46] To simplify this comparison, we exclude the thermostat energy so that any difference is solely due to the electronic structure solver[47] Bowler D R, Bush I J and Gillan M J 2000 Practical methods for ab initio calculations on thousands of atoms Int. J. Quant. Chem.77 831–42[48] Souvatzis P and Niklasson A M N 2014 First principles molecular dynamics without self-consistent field optimization J. Chem.Phys. 140 044117[49] Mujahed S Y, Miyazaki T and Bowler D R 2026 Zenodo (https://doi.org/10.5281/zenodo.21262374)11https://doi.org/10.1021/ct8005425https://doi.org/10.1021/ct8005425https://doi.org/10.1088/0953-8984/14/11/303https://doi.org/10.1088/0953-8984/14/11/303https://doi.org/10.1103/PhysRevB.47.10891https://doi.org/10.1103/PhysRevB.47.10891https://doi.org/10.1103/PhysRevB.50.17611https://doi.org/10.1103/PhysRevB.50.17611https://doi.org/10.1016/S0010-4655(99)00221-0https://doi.org/10.1016/S0010-4655(99)00221-0https://doi.org/10.1103/PhysRevB.58.12704https://doi.org/10.1103/PhysRevB.58.12704https://doi.org/10.1021/ct500847yhttps://doi.org/10.1021/ct500847yhttps://doi.org/10.1103/PhysRevB.50.4327https://doi.org/10.1103/PhysRevB.50.4327https://doi.org/10.1103/PhysRevB.64.165110https://doi.org/10.1103/PhysRevB.64.165110https://doi.org/10.1073/pnas.0505436102https://doi.org/10.1073/pnas.0505436102https://github.com/OrderN/CONQUEST-releasehttps://doi.org/10.1103/PhysRevB.43.1993https://doi.org/10.1103/PhysRevB.43.1993https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1103/PhysRevLett.77.3865https://doi.org/10.1016/j.cpc.2018.01.012https://doi.org/10.1016/j.cpc.2018.01.012https://doi.org/10.1103/PhysRevLett.100.123004https://doi.org/10.1103/PhysRevLett.100.123004https://doi.org/10.1088/1361-648X/aa810dhttps://doi.org/10.1088/1361-648X/aa810dhttps://doi.org/10.1002/(SICI)1097-461X(2000)77:53.0.CO;2-Ghttps://doi.org/10.1002/(SICI)1097-461X(2000)77:53.0.CO;2-Ghttps://doi.org/10.1063/1.4862907https://doi.org/10.1063/1.4862907https://doi.org/10.5281/zenodo.21262374 Stress calculation in linear scaling DFT: convergence and dynamics 1. Introduction 2. Energies 2.1. Energy in terms of eigenvalues and eigenstates 2.2. Energy in terms of the density matrix 2.3. Linear scaling and electron number 3. Formulating stress with localized orbitals and linear scaling 3.1. Hellmann–Feynman stress 3.2. Pulay stresses 3.3. Miscellaneous stresses 4. Tests and results 4.1. Convergence tests 4.2. Molecular dynamics tests 5. Conclusions References