Add breaks_TRS term predicate and guard compute_dynmat

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}

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]

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