From dfd1e2676cc8c31ba51eca4d2f7b5b84ac70d094 Mon Sep 17 00:00:00 2001 From: Antoine Levitt Date: Thu, 7 May 2026 14:20:37 +0200 Subject: [PATCH 1/3] Add breaks_TRS mechanism and assert it in compute_dynmat MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Mirror the existing breaks_symmetries dispatch with a breaks_TRS predicate on term builders. Mark Magnetic and Anyonic as breaking TRS, and assert in compute_dynmat that no term in the model breaks TRS — the phonon response solver folds the +q and -q Sternheimer equations into one using TRS, so it is currently invalid otherwise (see #1310, Dal Corso 2019). Also document the derivation of the "ψ̄·δψ with factor 2" form used in compute_δρ from scratch (no χ₀), tracking momenta term by term and showing how TRS identifies the missing δψ̄·ψ contribution with the computed one. --- src/densities.jl | 31 +++++++++++++++++++++++++++++++ src/postprocess/phonon.jl | 7 +++++++ src/terms/terms.jl | 8 ++++++++ 3 files changed, 46 insertions(+) diff --git a/src/densities.jl b/src/densities.jl index 1a110e2d98..05e7cd1939 100644 --- a/src/densities.jl +++ b/src/densities.jl @@ -55,6 +55,37 @@ using an optional `occupation_threshold`. By default all occupation numbers are end # Variation in density corresponding to a variation in the orbitals and occupations. +# +# Density response and the "factor of 2" / dropped δψ̄·ψ term +# --------------------------------------------------------- +# A real perturbation at wavevector q has two complex Fourier components, +# δV(r) = δV_q e^{iq·r} + δV_q* e^{−iq·r}. Each Bloch state ψ_k acquires a response +# with two pieces: δψ_k^{(+)} at momentum k+q (sourced by δV_q) and δψ_k^{(−)} at +# k−q (sourced by δV_q*). Differentiating ρ = ∑_k f_k |ψ_k|² gives +# δρ = ∑_k f_k (ψ̄_k δψ_k + δψ̄_k ψ_k). +# Picking out the +q Fourier component of δρ, by counting momenta: +# • ψ̄_k δψ_k^{(+)} contributes at −k + (k+q) = +q ← what we compute below +# • ψ̄_k δψ_k^{(−)} contributes at −q (drops out of δρ_{+q}) +# • δψ̄_k^{(+)} ψ_k contributes at −q (drops out) +# • δψ̄_k^{(−)} ψ_k contributes at +q ← *not* computed directly +# So naively δρ_{+q} would need both δψ_k^{(+)} (Sternheimer at +q) and δψ_k^{(−)} +# (Sternheimer at −q). +# +# Time-reversal symmetry rescues us. Under TRS the Bloch state ψ_{−k} = ψ_k* is +# also occupied with the same energy, and conjugating its Sternheimer equation +# at +q, +# (H_{−k+q} − ε_{−k}) δu_{−k}^{(+)} = −δV_q u_{−k}, +# turns it (using H̄_{−k+q} = H_{k−q} and ū_{−k} = u_k) into +# (H_{k−q} − ε_k) (δu_{−k}^{(+)})* = −δV_q* u_k, +# which is exactly the Sternheimer equation defining δu_k^{(−)}. Hence +# δψ_k^{(−)} = (δψ_{−k}^{(+)})*, and after summing over k the missing term +# δψ̄_k^{(−)} ψ_k just reproduces ψ̄_k δψ_k^{(+)}. The two +q contributions are +# therefore equal — we compute one and multiply by 2 below. +# +# Without TRS this identification fails and the −q Sternheimer equation must be +# solved separately (see JuliaMolSim/DFTK.jl#1310 and Dal Corso, +# https://arxiv.org/abs/1906.11673). compute_dynmat asserts !breaks_TRS to keep +# this routine valid in the phonon use-case. @views @timing function compute_δρ(basis::PlaneWaveBasis{T}, ψ, δψ, occupation, δoccupation=zero.(occupation); occupation_threshold=zero(T), q=zero(Vec3{T})) where {T} diff --git a/src/postprocess/phonon.jl b/src/postprocess/phonon.jl index 40a14f3d2f..260f34f623 100644 --- a/src/postprocess/phonon.jl +++ b/src/postprocess/phonon.jl @@ -71,6 +71,13 @@ in reduced coordinates. @timing function compute_dynmat(basis::PlaneWaveBasis{T}, ψ, occupation; q=zero(Vec3{T}), ρ=nothing, ham=nothing, εF=nothing, eigenvalues=nothing, kwargs...) where {T} + # The phonon response solver assumes time-reversal symmetry: the trick used + # to compute δρ from a single Sternheimer equation at +q (instead of one at + # +q and one at -q) is only valid under TRS. See the discussion in + # JuliaMolSim/DFTK.jl#1310 and Dal Corso, https://arxiv.org/abs/1906.11673. + @assert !any(breaks_TRS, basis.model.term_types) ( + "compute_dynmat is currently only implemented for time-reversal-symmetric " + * "Hamiltonians; see JuliaMolSim/DFTK.jl#1310.") n_atoms = length(basis.model.positions) δρs = [zeros(complex(T), basis.fft_size..., basis.model.n_spin_components) for _ = 1:3, _ = 1:n_atoms] diff --git a/src/terms/terms.jl b/src/terms/terms.jl index 6490d5d0cc..d5be103c99 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -50,6 +50,12 @@ include("Hamiltonian.jl") # is invalid) breaks_symmetries(::Any) = false +# breaks_TRS on a term builder answers true if this term breaks +# time-reversal symmetry. Some algorithms (notably the phonon response +# solver, see compute_dynmat) rely on TRS to avoid solving Sternheimer +# equations at both +q and -q (cf. discussion in JuliaMolSim/DFTK.jl#1310). +breaks_TRS(::Any) = false + include("kinetic.jl") include("local.jl") @@ -69,9 +75,11 @@ include("exact_exchange.jl") include("magnetic.jl") breaks_symmetries(::Magnetic) = true +breaks_TRS(::Magnetic) = true include("anyonic.jl") breaks_symmetries(::Anyonic) = true +breaks_TRS(::Anyonic) = true # forces computes either nothing or an array forces[at][α] (by default no forces) compute_forces(::Term, ::AbstractBasis, ψ, occupation; kwargs...) = nothing From 0bdbb6f713290a596ccd13722e2c45f1948a8003 Mon Sep 17 00:00:00 2001 From: Antoine Levitt Date: Thu, 7 May 2026 20:42:07 +0200 Subject: [PATCH 2/3] Rewrite --- src/densities.jl | 33 ++------------------------------- src/postprocess/phonon.jl | 29 ++++++++++++++++++++++++++--- src/terms/terms.jl | 12 ++++++------ 3 files changed, 34 insertions(+), 40 deletions(-) diff --git a/src/densities.jl b/src/densities.jl index 05e7cd1939..affbd8657e 100644 --- a/src/densities.jl +++ b/src/densities.jl @@ -55,37 +55,6 @@ using an optional `occupation_threshold`. By default all occupation numbers are end # Variation in density corresponding to a variation in the orbitals and occupations. -# -# Density response and the "factor of 2" / dropped δψ̄·ψ term -# --------------------------------------------------------- -# A real perturbation at wavevector q has two complex Fourier components, -# δV(r) = δV_q e^{iq·r} + δV_q* e^{−iq·r}. Each Bloch state ψ_k acquires a response -# with two pieces: δψ_k^{(+)} at momentum k+q (sourced by δV_q) and δψ_k^{(−)} at -# k−q (sourced by δV_q*). Differentiating ρ = ∑_k f_k |ψ_k|² gives -# δρ = ∑_k f_k (ψ̄_k δψ_k + δψ̄_k ψ_k). -# Picking out the +q Fourier component of δρ, by counting momenta: -# • ψ̄_k δψ_k^{(+)} contributes at −k + (k+q) = +q ← what we compute below -# • ψ̄_k δψ_k^{(−)} contributes at −q (drops out of δρ_{+q}) -# • δψ̄_k^{(+)} ψ_k contributes at −q (drops out) -# • δψ̄_k^{(−)} ψ_k contributes at +q ← *not* computed directly -# So naively δρ_{+q} would need both δψ_k^{(+)} (Sternheimer at +q) and δψ_k^{(−)} -# (Sternheimer at −q). -# -# Time-reversal symmetry rescues us. Under TRS the Bloch state ψ_{−k} = ψ_k* is -# also occupied with the same energy, and conjugating its Sternheimer equation -# at +q, -# (H_{−k+q} − ε_{−k}) δu_{−k}^{(+)} = −δV_q u_{−k}, -# turns it (using H̄_{−k+q} = H_{k−q} and ū_{−k} = u_k) into -# (H_{k−q} − ε_k) (δu_{−k}^{(+)})* = −δV_q* u_k, -# which is exactly the Sternheimer equation defining δu_k^{(−)}. Hence -# δψ_k^{(−)} = (δψ_{−k}^{(+)})*, and after summing over k the missing term -# δψ̄_k^{(−)} ψ_k just reproduces ψ̄_k δψ_k^{(+)}. The two +q contributions are -# therefore equal — we compute one and multiply by 2 below. -# -# Without TRS this identification fails and the −q Sternheimer equation must be -# solved separately (see JuliaMolSim/DFTK.jl#1310 and Dal Corso, -# https://arxiv.org/abs/1906.11673). compute_dynmat asserts !breaks_TRS to keep -# this routine valid in the phonon use-case. @views @timing function compute_δρ(basis::PlaneWaveBasis{T}, ψ, δψ, occupation, δoccupation=zero.(occupation); occupation_threshold=zero(T), q=zero(Vec3{T})) where {T} @@ -122,6 +91,8 @@ end # use unnormalized plans for extra speed, normalize at the end ifft_normalization = basis.fft_grid.ifft_normalization + # For the (non-trivial) explanation of why we don't take real parts when q!=0, + # see comment in phonon.jl storage.δρ[:, :, :, kpt.spin] .+= ifft_normalization^2 .* real_qzero.( 2 .* occupation[ik][n] .* basis.kweights[ik] .* conj.(storage.ψnk_real) .* storage.δψnk_real diff --git a/src/postprocess/phonon.jl b/src/postprocess/phonon.jl index 260f34f623..589d3b2520 100644 --- a/src/postprocess/phonon.jl +++ b/src/postprocess/phonon.jl @@ -1,3 +1,27 @@ +# Calculation of phonons from DFPT. +# +# This implementation relies on time-reversal symmetry in the following way. +# A real perturbation at wavevector q has two complex Fourier components, +# δV(r) = δVq e^{iq·r} + δVq* e^{-iq·r}. Each Bloch state ψk acquires a response +# δψk = δψk+ + δψk- with two pieces: δψk+ at momentum k+q, and δψk- at momentum k-q. +# Differentiating ρ = ∑_k fk |ψk|² gives +# δρ = ∑_k fk (ψk* δψk + δψk* ψk). +# Because of the δψk*, everything is coupled: one cannot just +# separate the + and - parts, and would need two Sternheimer solves. +# More abstractly, the map δV -> δρ that one gets naively by +# δρ = ∑_k fk (ψk* δψk + δψk* ψk) is R-linear but not C-linear +# so it's not valid to close our eyes and take δVq e^{iq·r} as a perturbation +# (similar structure to Casida equations in TDDFT). + +# Under time-reversal symmetry (TRS) however, the contribution to δρ of +k and -k are linked: +# δψk+ = (δψ(-k)-)*. +# so +# δρ = 2 ∑_k fk (ψk* δψk) (this is the central equation) +# In this form, the map δV -> δρ becomes complex-linear and only one Sternheimer per kpoint is needed +# Without TRS this fails and one needs two Sternheimer (see +# JuliaMolSim/DFTK.jl#1310 and eg Dal Corso, +# https://arxiv.org/abs/1906.11673). + # Convert to Cartesian a dynamical matrix in reduced coordinates. function dynmat_red_to_cart(model::Model, dynmat) inv_lattice = model.inv_lattice @@ -75,9 +99,8 @@ in reduced coordinates. # to compute δρ from a single Sternheimer equation at +q (instead of one at # +q and one at -q) is only valid under TRS. See the discussion in # JuliaMolSim/DFTK.jl#1310 and Dal Corso, https://arxiv.org/abs/1906.11673. - @assert !any(breaks_TRS, basis.model.term_types) ( - "compute_dynmat is currently only implemented for time-reversal-symmetric " - * "Hamiltonians; see JuliaMolSim/DFTK.jl#1310.") + @assert !any(breaks_time_reversal_symmetry, basis.model.term_types) ( + "Phonons are currently only implemented in the presence of time-reversal-symmetry.") n_atoms = length(basis.model.positions) δρs = [zeros(complex(T), basis.fft_size..., basis.model.n_spin_components) for _ = 1:3, _ = 1:n_atoms] diff --git a/src/terms/terms.jl b/src/terms/terms.jl index d5be103c99..fc0a8ad9fe 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -51,10 +51,10 @@ include("Hamiltonian.jl") breaks_symmetries(::Any) = false # breaks_TRS on a term builder answers true if this term breaks -# time-reversal symmetry. Some algorithms (notably the phonon response -# solver, see compute_dynmat) rely on TRS to avoid solving Sternheimer -# equations at both +q and -q (cf. discussion in JuliaMolSim/DFTK.jl#1310). -breaks_TRS(::Any) = false +# time-reversal symmetry. Phonon computations rely implicitly on TRS +# to avoid solving Sternheimer equations at both +q and -q (cf. +# discussion in JuliaMolSim/DFTK.jl#1310). +breaks_time_reversal_symmetry(::Any) = false include("kinetic.jl") @@ -75,11 +75,11 @@ include("exact_exchange.jl") include("magnetic.jl") breaks_symmetries(::Magnetic) = true -breaks_TRS(::Magnetic) = true +breaks_time_reversal_symmetry(::Magnetic) = true include("anyonic.jl") breaks_symmetries(::Anyonic) = true -breaks_TRS(::Anyonic) = true +breaks_time_reversal_symmetry(::Anyonic) = true # forces computes either nothing or an array forces[at][α] (by default no forces) compute_forces(::Term, ::AbstractBasis, ψ, occupation; kwargs...) = nothing From c914dff4e83cc963bc6583f4264d3701d22aac5d Mon Sep 17 00:00:00 2001 From: Antoine Levitt Date: Fri, 8 May 2026 21:27:21 +0200 Subject: [PATCH 3/3] review --- src/terms/terms.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/terms/terms.jl b/src/terms/terms.jl index fc0a8ad9fe..f18cd6ac39 100644 --- a/src/terms/terms.jl +++ b/src/terms/terms.jl @@ -50,7 +50,7 @@ include("Hamiltonian.jl") # is invalid) breaks_symmetries(::Any) = false -# breaks_TRS on a term builder answers true if this term breaks +# breaks_time_reversal_symmetry on a term builder answers true if this term breaks # time-reversal symmetry. Phonon computations rely implicitly on TRS # to avoid solving Sternheimer equations at both +q and -q (cf. # discussion in JuliaMolSim/DFTK.jl#1310).