From 0124e9af95bbb06805bc8eeb78d7a09b56c52672 Mon Sep 17 00:00:00 2001 From: Antoine Levitt Date: Thu, 16 Jul 2026 11:43:59 +0200 Subject: [PATCH] Fix Float32 precision loss in lda_c_pw under AD MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In `energy_per_particle_c_pw` the two `(1 ± ζ)` bases of equation (9) are already typed `T`, but the remaining base was the integer literal `2`. With a `Dual` exponent an integer base is differentiated via `log(2)` in Float64, so a `Dual{Tag,Float32}` density came back as `Dual{Tag,Float64}`. Only `energy_density` was affected: `potential_terms` pins the type via `arithmetic_type` + `similar(ρ, TT, …)`. Downstream in DFTK this made the XC energy the only Float64-valued term in an otherwise Float32 forward-mode evaluation, which broke Float32 stress computations. The change is value-neutral for non-`Dual` types (`2^x === T(2)^x` bit-identical for Float32 and Float64); for a plain `Float32`, `promote(2, ::Float32)` is already `Float32`, which is why this only showed up under AD. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_017an8Gie5b1ULeSjhep3Sjb --- src/functionals/lda.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/functionals/lda.jl b/src/functionals/lda.jl index 04a4632..f5de403 100644 --- a/src/functionals/lda.jl +++ b/src/functionals/lda.jl @@ -83,7 +83,7 @@ function energy_per_particle_c_pw(ρ::T; A, α₁, β₁, β₂, β₃, β₄, f end # equation (9) - f(ζ) = ((1 + ζ)^(4 / T(3)) + (1 - ζ)^(4 / T(3)) - 2) / (2^(4 / T(3)) - 2) # == 0 for non-spin-polarised + f(ζ) = ((1 + ζ)^(4 / T(3)) + (1 - ζ)^(4 / T(3)) - 2) / (T(2)^(4 / T(3)) - 2) # == 0 for non-spin-polarised ε_0(rₛ) = G(sqrt(rₛ), A[1], α₁[1], β₁[1], β₂[1], β₃[1], β₄[1]) # ε_c(rₛ, 0) ε_1(rₛ) = G(sqrt(rₛ), A[2], α₁[2], β₁[2], β₂[2], β₃[2], β₄[2]) # ε_c(rₛ, 1)