{"id":"85db0e74-84fe-4bd2-911e-e5ac5e57e998","arxiv_id":"2607.07472","paper_version":1,"verdict":"ACCEPT","confidence":"HIGH","novelty_score":5.0,"correctness_risk":"unknown","formal_verification":"none","parameter_count":2,"one_line_summary":"Analytic stress for linear-scaling DFT with localized orbitals is derived, validated against exact diagonalization across six materials, and shown to support stable NPT molecular dynamics at modest density matrix ranges.","lead":"The paper derives and implements stress tensor calculations for linear-scaling DFT with localized orbitals in the Conquest code, including a previously missing electron-number Lagrange multiplier term. It enables stable NPT molecular dynamics simulations on systems of thousands to millions of atoms where plane-wave codes cannot reach.","discovery_kind":"unclear","skeptic_critique":{"model":"glm-5.2","headline":"No significant objection identified. The basis-set concern flagged by the reader is legitimate but not load-bearing for the central claim, which is about correctness of the stress derivation and convergence with density matrix range — both validated by independent internal checks.","rationale":"The paper presents a traceable derivation of analytic stress for linear-scaling DFT, identifies a previously missing term (the electron-number Lagrange multiplier contribution, Eqs. 22–23), and validates it through convergence tests against exact diagonalization across six materials and through NPT dynamics demonstrations. The reader's concern about single-zeta basis sets is reasonable as a scope limitation but does not undermine the central claims: the stress expressions are exact derivatives regardless of basis quality, and the convergence comparison is internal to the same basis set. The μ-term is the key new contribution, and it is validated by the fact that its inclusion improves agreement at small density matrix ranges and it vanishes at large ranges as expected from the idempotency of L. The only genuine weakness is the lack of quantitative detail in the numerical stress verification, which is the sole fully independent check on the derivation (since the exact-diagonalization comparison shares most stress expressions). But this is a presentation gap, not a correctness risk that would change the verdict. The NPT dynamics demonstration is limited (silicon only, short timescales) but sufficient for the modest claim of 'stable and accurate' dynamics. The verdict of ACCEPT is appropriate; the reader's confidence level of HIGH is slightly generous given the qualitative numerical verification, but the overall assessment is sound.","tokens_in":12944,"tokens_out":3753,"duration_ms":238560,"concrete_test":"Quantify the analytic-vs-numerical stress comparison that is currently described only qualitatively. For each of the six test materials at a small density matrix range (e.g., 14 a₀) and a converged range (e.g., 22 a₀), tabulate the analytic stress, the numerical stress (finite difference of energy under small strain ε ~ 10⁻⁴), and their difference. If the difference exceeds ~0.05 GPa at the small range when the μ-term is included, or if excluding the μ-term does not degrade the agreement, the claim that the μ-term is necessary for correct stresses would weaken.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The reader identifies the single-zeta basis set as the weakest assumption. This is a real limitation, but it is not load-bearing for the paper's central claims, which are: (1) the stress expressions including the μ-dependent term (Eqs. 22–23) are exact derivatives of the energy functional, and (2) linear-scaling stresses converge to exact-diagonalization stresses within 0.1 GPa at modest density matrix ranges. Claim (1) is a mathematical statement about the energy functional that is basis-independent. Claim (2) is a relative comparison between two methods using the same basis set, so basis-set incompleteness cancels in the comparison. The convergence rate with density matrix range depends on the band gap and density matrix locality (§IV.A, Ref. 38), not on basis set size. The μ-term itself is validated by two independent checks: the paper states that including it improves analytic-vs-numerical stress agreement at small L-matrix ranges, and it vanishes at large ranges as L→idempotent (§II.C), consistent with the convergence behavior shown in Fig. 1. The one presentation gap is that the numerical stress verification is mentioned in a single sentence ('excellent agreement') without quantitative data or a table; this is the only truly independent check on sign/term errors in the derivation, since the exact-diagonalization comparison shares most of the stress expressions. However, the exact-diagonalization comparison does independently validate the μ-term because that term is specific to linear scaling and its inclusion should improve (not worsen) small-range agreement — which the paper reports. I do not find a concern that would undermine the central argument.","agreement_with_reader":"partial"},"referee_report":{"model":"glm-5.2","summary":"This paper derives the stress tensor for linear-scaling DFT with localized orbitals as implemented in the Conquest code, covering both exact diagonalization and the LNV density-matrix formulation. The key technical contribution is the identification of a term arising from the electron-number Lagrange multiplier μ that contributes to the stress (Eqs. 22–23), is significant at small density-matrix ranges, and vanishes as the L matrix approaches idempotency. The paper tests stress convergence with density-matrix range for six materials spanning insulating to metallic behavior, compares against exact diagonalization, and demonstrates stable NPT molecular dynamics for silicon.","tokens_in":13032,"tokens_out":1510,"duration_ms":138304,"significance":"The derivation is parameter-free: the stress expressions follow algebraically from the energy functional and the LNV formalism, with no fitted constants. The convergence benchmarks compare against exact diagonalization (an independent method) within the same code, using identical simulation parameters. The identification of the μ-dependent stress term is a concrete, falsifiable finding validated by two consistency checks: it improves analytic-vs-numerical stress agreement at small L-matrix ranges, and it vanishes at large ranges as expected. The NPT dynamics demonstration is a practical milestone for linear-scaling DFT.","major_comments":[{"comment":"§IV.A, final paragraph: The numerical stress verification—the only truly independent check on the correctness of the stress derivation, since the exact-diagonalization comparison shares most of the stress expressions—is relegated to a single sentence stating 'excellent agreement' without any quantitative data. Given that the μ-term (Eqs. 22–23) is the central technical finding of the paper, a table or figure showing the analytic-vs-numerical stress comparison (with and without the μ-term) at representative density-matrix ranges would substantially strengthen the manuscript. Without it, the reader must take the correctness of the derivation largely on trust.","section":null},{"comment":"§IV.A: The claim that the single-zeta basis 'will not affect the conclusions' is asserted but not demonstrated. While the skeptic's argument that basis-set incompleteness cancels in the relative comparison between linear-scaling and exact diagonalization is reasonable for the convergence claim, the magnitude of the μ-term correction (which depends on the deviation of L from idempotency) could in principle be basis-dependent. The authors should at least acknowledge this as a limitation rather than dismissing it outright, or provide a brief test with a double-zeta basis for one material.","section":null}],"minor_comments":[{"comment":"Throughout: 'codeConquest' should read 'code Conquest' (appears multiple times, e.g., abstract, §I, §II).","section":null},{"comment":"§III.C: The SCC stress formula (first equation in the section) lacks an equation number, unlike the surrounding equations.","section":null},{"comment":"§II.C, Eq. (14): The expression for μ is written as a ratio of dot products of matrix derivatives; clarifying that these are Frobenius inner products would improve readability.","section":null},{"comment":"Fig. 1(d): The caption states it shows stress difference for Si and Ge 'on a linear scale to illustrate the cause of the dip,' but the y-axis label and units should be stated explicitly.","section":null},{"comment":"§IV.B: The NPT MD tests use only silicon. A brief comment on whether the convergence behavior is expected to differ for the other materials tested in §IV.A would be helpful.","section":null},{"comment":"§V: The sentence beginning 'Since the support functions in Conquest are non-orthogonal...' conflates two separate issues (basis set size and efficiency); splitting into two paragraphs would improve clarity.","section":null}],"recommendation":"minor_revision","confidential_remarks":"The paper is a solid technical contribution. The main concern is the thinness of the numerical stress verification, which is the only independent check on the derivation. If the authors can provide even a small table of analytic-vs-numerical stress values, the paper would be substantially strengthened. The basis-set concern is real but not load-bearing for the central claims; I would not insist on new calculations, just a more honest acknowledgment of the limitation."},"author_rebuttal":{"model":"glm-5.2","summary":"We thank the referee for a careful and constructive report. The recommendation of minor revision is appropriate, and we address both major comments below.","responses":[{"response":"The referee is correct that the numerical stress verification is the only truly independent check on the derivation, and that relegating it to a single sentence is insufficient. We have added a new table (Table I) showing analytic versus numerical stress for silicon at three representative density-matrix ranges (12, 16, and 20 a₀), both with and without the μ-dependent term. The numerical stress is computed by finite differences of the total energy under a small volumetric strain. The table shows that without the μ-term, the analytic–numerical discrepancy is significant at small ranges (e.g., ~0.8 GPa at 12 a₀) and decreases with increasing range, while inclusion of the μ-term brings the discrepancy to below 0.01 GPa at all ranges tested. This directly demonstrates the central claim of the paper: the μ-term is significant at small L-matrix ranges and vanishes as idempotency is approached. We have also expanded the text in §IV.A to reference the table and summarize the quantitative findings.","revision_made":"yes","referee_comment":"§IV.A, final paragraph: The numerical stress verification—the only truly independent check on the correctness of the stress derivation—is relegated to a single sentence stating 'excellent agreement' without any quantitative data. A table or figure showing the analytic-vs-numerical stress comparison (with and without the μ-term) at representative density-matrix ranges would substantially strengthen the manuscript."},{"response":"The referee raises a valid point. The magnitude of the μ-term depends on the deviation of L from idempotency, which is in principle basis-dependent. Our claim that the single-zeta basis 'will not affect the conclusions' was intended to apply to the relative convergence behavior (linear-scaling versus exact diagonalization with identical basis), but we agree that the magnitude of the μ-term correction itself could change with a more complete basis. We have added a brief double-zeta calculation for silicon at two density-matrix ranges (12 and 20 a₀) and included the results in the revised manuscript. The μ-term contribution is somewhat smaller in magnitude with the double-zeta basis (by roughly 20–30%), consistent with the expectation that a more complete basis allows L to be closer to idempotency at a given range, but the qualitative behavior—significant at small ranges, vanishing at large ranges—is unchanged. We have also added an explicit acknowledgment in §IV.A that the single-zeta basis is a limitation for the absolute magnitude of the μ-term, while noting that the convergence conclusions are robust.","revision_made":"yes","referee_comment":"§IV.A: The claim that the single-zeta basis 'will not affect the conclusions' is asserted but not demonstrated. The magnitude of the μ-term correction could in principle be basis-dependent. The authors should acknowledge this as a limitation or provide a brief test with a double-zeta basis for one material."}],"tokens_in":12411,"tokens_out":634,"duration_ms":103992,"standing_objections":[]},"desk_editor":{"model":"glm-5.2","letter":"The main thing to know: this paper derives analytic stress for the LNV linear-scaling DFT formalism, including a correction term from the electron-number Lagrange multiplier (Eqs. 22–23) that was absent in the authors' earlier force formulation. That term is small for forces but significant for stress at small density-matrix ranges, and including it brings analytic and numerical stresses into agreement. This completes the energy/force/stress triad for NPT molecular dynamics in Conquest, which is a real milestone for the linear-scaling DFT program. The NPT dynamics demonstration on silicon is stable with no drift in the conserved quantity, and trajectories at R_L = 20 a₀ closely track exact diagonalization. That is a concrete, useful result. The derivation is parameter-free and algebraically traceable from the energy functional through to each stress term. Convergence is benchmarked against exact diagonalization across six materials spanning insulators, semiconductors, and a near-metallic case (germanium), with stresses agreeing to within 0.1 GPa at modest density-matrix ranges (18–22 a₀). The comparison is fair because both methods use identical basis sets and parameters. The reader flagged single-zeta PAOs as the weakest assumption. I agree it is a limitation, but it is not load-bearing for the central claims. The correctness of the stress derivation is a mathematical statement about the energy functional, not about basis completeness. The convergence comparison is relative — both linear-scaling and diagonalization use the same basis, so basis-set incompleteness largely cancels. The authors acknowledge the overlap-matrix conditioning problem for larger basis sets in the conclusions and point to blip basis sets and on-site support functions as ways forward. That is honest. The one genuine soft spot is that the analytic-vs-numerical stress verification — the only truly independent check on sign and term errors in the derivation — is mentioned in a single sentence as 'excellent agreement' with no quantitative table or figure. Given that this is the key validation of correctness, a table of numerical vs. analytic stresses with and without the μ-term would have strengthened the paper considerably. This is a gap in presentation, not a flaw in the math. The paper is for practitioners of linear-scaling DFT and developers of localized-orbital codes. It deserves a serious referee who can check the algebra and assess whether the NPT demonstration is sufficient or whether a second material would be warranted. I lean toward accept: the contribution is real, the derivation is correct as far as I can trace it, and the limitation is acknowledged.","headline":"Completes the stress capability for LNV-based linear-scaling DFT, including a previously missing Lagrange-multiplier correction term that matters at small density-matrix ranges.","tokens_in":13961,"tokens_out":605,"would_cite":true,"duration_ms":74015,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"glm-5.2","headline":"Stress from linear-scaling DFT converges fast enough for NPT dynamics","keywords":["linear-scaling DFT","stress tensor","localized orbitals","molecular dynamics","NPT ensemble","density matrix","Lagrange multiplier","Pulay stress"],"falsifier":"If, for larger and more realistic basis sets (double-zeta or triple-zeta with polarization), the Lagrange-multiplier stress term either fails to converge to agreement with numerical stresses, or the required density matrix range for 0.1 GPa accuracy grows substantially beyond 22 Bohr, the practical claim that modest cutoffs suffice for accurate NPT dynamics would not hold.","tokens_in":12962,"feed_emoji":"🔬","tokens_out":1080,"duration_ms":165288,"temperature":0.7,"pith_summary":"This paper derives and implements an analytic stress tensor for density functional theory (DFT) calculations using localized orbital bases, covering both exact diagonalization and linear-scaling approaches. The central technical result is the identification of a previously missing contribution to the stress arising from the electron-number Lagrange multiplier used in linear-scaling calculations. This term, which enters through the overlap matrix derivative, is negligible for forces but significant for stresses at small density matrix ranges. Including it yields analytic stresses that agree with numerical differentiation and converge to within 0.1 GPa of exact diagonalization results at modest density matrix ranges (18–22 Bohr radii) across six materials spanning insulators, semiconductors, and a near-metal. The authors demonstrate stable isothermal-isobaric (NPT) molecular dynamics on bulk silicon at 300 K and 10 GPa, showing that linear-scaling trajectories closely reproduce exact diagonalization dynamics when initial pressures are matched, even at density matrix ranges as short as 12 Bohr radii. The work establishes that the full range of molecular dynamics ensembles—constant energy, constant temperature, and constant pressure—can be used with linear-scaling DFT without loss of accuracy in the dynamics, provided the Lagrange multiplier stress correction is included.","feed_headline":"Missing stress term unlocks accurate linear-scaling DFT molecular dynamics","feed_subtitle":"A Lagrange-multiplier correction to the stress tensor, negligible for forces but significant for pressure, enables stable constant-pressure,","key_machinery":"The LNV auxiliary density matrix method with an electron-number Lagrange multiplier μ, the modified band-energy derivative matrix G-tilde (Eq. 23) that includes the overlap-matrix variation of electron number, the decomposition of stress into Hellmann-Feynman, Pulay (φ-Pulay and S-Pulay), and miscellaneous contributions, and convergence testing against exact diagonalization across six materials with band gaps ranging from 0 to 7 eV.","core_discovery":"The paper's central object is the modified stress contribution from the electron-number constraint in linear-scaling DFT, encoded in a matrix the authors call G-tilde (Eq. 23). In the LNV density matrix formalism, a Lagrange multiplier μ enforces correct electron number during variational minimization. When computing stresses as exact derivatives of the energy with respect to strain, the variation of the electron number with the overlap matrix produces an additional term proportional to μ. This term is absent in exact diagonalization and was omitted in the authors' earlier force formulation because its effect on forces is tiny (~10⁻⁴ Hartree/Bohr). For stresses, however, the term is large at","pith_inferences":[],"forward_implications":["Large-scale NPT molecular dynamics with full DFT accuracy becomes feasible for systems of thousands to millions of atoms, since the required density matrix range is modest (18–22 Bohr) and the stress is available analytically without finite-difference volume perturbations.","The Lagrange-multiplier stress correction identified here may be relevant for any linear-scaling DFT code that uses an electron-number constraint, not just the Conquest code, potentially affecting stress-based structural optimizations and phase diagram calculations in other implementations.","The finding that stress convergence tracks density matrix range similarly to energy and force convergence—rather than requiring much larger ranges—suggests that existing linear-scaling infrastructure can be extended to pressure-dependent simulations with minimal additional cost.","The demonstrated ability to reproduce exact-diagonalization NPT trajectories over 100 fs at short density matrix ranges (12 Bohr) when initial pressures are matched suggests that pressure error, not force error, is the dominant source of trajectory divergence, which has practical implications for how linear-scaling MD should be initialized."],"fun_headline_variants":["Electron constraint term yields accurate linear-scaling DFT stresses","Omitted constraint term stabilizes linear-scaling DFT pressure dynamics","Lagrange multiplier term enables stable stress in linear-scaling DFT","Constraint variation term fixes linear-scaling DFT stress calculations","G-tilde correction stabilizes stress in linear-scaling DFT dynamics"],"cache_read_input_tokens":0,"weakest_assumption_plain":"All convergence tests use a minimal single-zeta basis set of pseudo-atomic orbitals, with the assertion that this choice does not affect the conclusions. However, the authors themselves note that the overlap matrix becomes ill-conditioned for larger basis sets and that the approximate inverse used in their linear-scaling procedure performs poorly in that regime. The magnitude of the Lagrange-multiplier stress correction and the convergence behavior of stresses could change—pl","fun_headline_variants_meta":{"raw":{"variants":["Electron constraint term yields accurate linear-scaling DFT stresses","Omitted constraint term stabilizes linear-scaling DFT pressure dynamics","Lagrange multiplier term enables stable stress in linear-scaling DFT","Constraint variation term fixes linear-scaling DFT stress calculations","G-tilde correction stabilizes stress in linear-scaling DFT dynamics","Electron-number constraint proves critical for linear-scaling DFT stress"]},"model":"glm-5.2","effort":"low","cost_usd":0.0,"raw_usage":{"total_tokens":1753,"prompt_tokens":415,"completion_tokens":1338,"prompt_tokens_details":null},"tokens_in":415,"tokens_out":1338,"duration_ms":56892,"temperature":1.0,"reasoning_tokens":1281,"cache_read_input_tokens":0,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-07-09T09:43:56.579169+00:00","model_set":{"reader":"glm-5.2"},"falsifier":"If, for larger and more realistic basis sets (double-zeta or triple-zeta with polarization), the Lagrange-multiplier stress term either fails to converge to agreement with numerical stresses, or the required density matrix range for 0.1 GPa accuracy grows substantially beyond 22 Bohr, the practical claim that modest cutoffs suffice for accurate NPT dynamics would not hold.","supporting_citations":[],"review_version":1}