Add breaks_TRS term predicate and guard compute_dynmat

src/densities.jl

@@ -91,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

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
@@ -71,6 +95,12 @@ 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_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]

src/terms/terms.jl

@@ -50,6 +50,12 @@ include("Hamiltonian.jl")
 # is invalid)
 breaks_symmetries(::Any) = false
 
+# 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).
+breaks_time_reversal_symmetry(::Any) = false
+
 include("kinetic.jl")
 
 include("local.jl")
@@ -69,9 +75,11 @@ include("exact_exchange.jl")
 
 include("magnetic.jl")
 breaks_symmetries(::Magnetic) = true
+breaks_time_reversal_symmetry(::Magnetic) = true
 
 include("anyonic.jl")
 breaks_symmetries(::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