Symmetries and response

[antoine-levitt]

We need a consistent plan for using symmetries in response calculations. What follows is a discussion with Claude about the design, as a basis for discussion. Some of it is slop but I think the general diagnosis and proposed design sounds reasonable

Symmetry in linear response — a unified plan

Status: design proposal, for discussion. Nothing here is implemented yet.

Assumed baseline: this document assumes the time-reversal-symmetry PR (antiunitary SymOp,
θ tag, breaks_conjugation) and the electric-field dielectric PR (shift_kpoints,
compute_dielectric, compute_δHψ_dk) are both merged. It builds on top of them.

Audience. Someone who knows the DFTK codebase but finds the symmetry handling perpetually
confusing. So §1–§2 rederive the symmetry machinery from scratch with DFTK's exact conventions,
before §3 onward proposes the new abstractions. If you already dream in SymOps, skim to §3.

Goal. One mechanism, respecting the math, that lets every density-functional perturbation
theory (DFPT) response — dielectric (q=0 electric field), phonons (q≠0 atomic displacement),
strain, and future variants — exploit crystal symmetry to reduce the k-point sampling, without
each response reimplementing symmetry. Today the dielectric and phonon codes both punt with
symmetries=false; this plan removes that.

---

1. The symmetry of a periodic Hamiltonian, from scratch

1.1 What a symmetry is

A crystal has a finite group of space-group operations. Each is a pair (W, w) acting on
real-space (reduced/fractional) coordinates as

x  ↦  W x + w

with W an integer matrix (a point operation in the reduced basis) and w a translation, such
that the operation maps every atom onto an atom of the same species. In DFTK this is a SymOp
(src/SymOp.jl), which precomputes some redundant but convenient quantities:

struct SymOp{T}
    W::Mat3{Int}   # real-space:  x ↦ W x + w
    w::Vec3{T}
    S::Mat3{Int}   # = W'  (reciprocal-space point operation)
    τ::Vec3{T}     # = -W⁻¹ w
    θ::Int         # +1 unitary, -1 antiunitary (see §1.4)
end

A SymOp induces a unitary operator U on wavefunctions,

(U u)(x) = u(W x + w).                                    (θ = +1)

The reason S = Wᵀ and τ = -W⁻¹w appear is that in Fourier space this operator reads

(Û u)(G) = e^{-i G·τ} û(S⁻¹ G).                          (θ = +1)

(All these identities hold in reduced coordinates too, dropping factors of 2π; DFTK works in
reduced coordinates, so G·τ below is a plain dot product of the integer G with τ.)

1.2 Why symmetry lets us skip k-points

The self-consistent potential of a symmetric crystal is invariant under U, so U commutes
with the Hamiltonian: U H U⁻¹ = H. Because a Bloch state at k has crystal momentum k, and
U moves momentum k → S k, commutation means

H_{Sk}  is unitarily equivalent to  H_k,     with     u_{Sk} = U u_k.

So the eigenvalues at Sk equal those at k, and the eigenvectors are related by a known
operator U. We never have to diagonalize at Sk if we have already done k. This is the
entire point of symmetry: we sample only a set of irreducible k-points (the IBZ), and every
other k-point in the full BZ is reconstructed by applying some U.

Three nested groups matter (decreasing size):

1. symmetries of the lattice (ignores atoms),
2. symmetries of the crystal = model.symmetries (lattice ops that also map atoms to atoms),
3. symmetries that additionally preserve the discrete BZ mesh and FFT grid =
   basis.symmetries (a subset of model.symmetries).

basis.symmetries is the one used for k-point reduction. It is computed in the PlaneWaveBasis
constructor by filtering model.symmetries through symmetries_preserving_rgrid and
symmetries_preserving_kgrid (src/symmetry.jl, src/PlaneWaveBasis.jl).

1.3 Symmetrizing observables

An observable computed from IBZ data alone is not symmetric — it must be averaged over the group
to restore the full-BZ result. For a density (real, periodic) the average is, in Fourier space,

ρ_sym(G) = (1/|G_group|) Σ_S  e^{-i G·τ_S}  ρ(S⁻¹ G).

This is exactly accumulate_over_symmetries! (the sum) and symmetrize_ρ (the normalized
average) in src/symmetry.jl. Other observables symmetrize with their own representation:

• forces F_s (one vector per atom): average inv(Wᵀ)·F with atoms permuted by the op
  (symmetrize_forces).
• stresses σ (a rank-2 tensor): average W_cart σ W_cartᵀ (symmetrize_stresses).
• Hubbard occupations n_s (a matrix per atom/manifold): average D_l(W)ᵀ n D_l(W) with a
  Wigner-D matrix and atoms permuted (symmetrize_hubbard_n).

Key mental model. "Symmetrizing" is always the same operation: push the object through every
group element (spatially, and by whatever representation the object carries) and average. Densities
carry the trivial representation on top of the spatial action; forces carry the vector
representation; Hubbard occupations carry the Wigner-D representation. Keep this in mind — the whole
plan is "do this for the response density, with the representation the perturbation carries."

1.4 Time reversal as an antiunitary symmetry (the TRS PR)

Complex conjugation K (with Ku = ū) maps a Bloch state at k to one at -k with the same
energy — it is an antiunitary symmetry whenever no term breaks it (no external magnetic field,
no anyons). The TRS PR folds this into SymOp by tagging each op with θ ∈ {+1, -1}:

(U u)(x) = u(W x + w)          (θ=+1)          (U u)(x) = conj(u(W x + w))     (θ=-1)
(Û u)(G) = e^{-iG·τ} û(S⁻¹G)   (θ=+1)          (Û u)(G) = e^{+iG·τ} conj(û(-S⁻¹G))   (θ=-1)

The single unifying statement is that the BZ action of a symop is k ↦ θ·S·k:

transform_kpoint_coordinate(op, k) = op.θ * op.S * k

θ=+1 gives the usual Sk; θ=-1 gives -Sk, i.e. the conjugation partner. model.symmetries
is now the spatial group augmented with its θ=-1 partners when time reversal is unbroken
(default_symmetries → symmetry_operations(...; time_reversal=true)). The rule for everyone else
in the code is simple and mechanical:

Never write symop.S * k. Always write transform_kpoint_coordinate(symop, k).

For a real density, conjugation partners are redundant (ρ̂(-G) = conj(ρ̂(G)) = ρ̂(G) since ρ is
real), so symmetrize_ρ drops θ=-1 ops as an optimization. This optimization is specific to
real quantities — it will not carry over to the complex response densities of §4.

1.5 How a term breaks symmetry today

A term can break symmetry. The mechanism is two boolean, per-class traits on the term builder
(src/terms/terms.jl):

breaks_symmetries(::Magnetic)  = true    # breaks the spatial point group
breaks_conjugation(::Magnetic) = true    # breaks time reversal (θ = -1)

default_symmetries consults these: an external potential / magnetic field / anyonic term drops the
crystal to the trivial group (or drops just the θ=-1 ops). This is coarse — all or nothing per
class — but it is the right granularity for terms, which define the model group. The finer,
per-operation reasoning belongs to perturbations (§3), which carve a subgroup out of the model
group. Keep these two levels separate; conflating them is the usual source of confusion.

---

2. What linear response computes, from scratch

DFPT computes the first-order change of the electronic structure under a perturbation. The objects:

• a perturbing potential δV, at some wavevector q;
• the induced orbital response δψ (the change in the occupied Bloch states);
• the induced density response δρ;
• a contracted observable (dielectric tensor, dynamical matrix, ...).

The self-consistent (Dyson / Sternheimer) machinery to get δψ/δρ from δV already exists in
src/response/ (solve_ΩplusK_split, apply_χ0, apply_χ0_4P) and is symmetry-agnostic per
k-point. The only places symmetry enters the response are:

• compute_δρ ends with symmetrize_ρ(basis, δρ) (reconstruct full-BZ δρ from IBZ);
• apply_χ0 symmetrizes δV with symmetrize_ρ (inside solve_ΩplusK_split via
  DielectricAdjoint).

Both key off basis.symmetries. This is the leverage point of the whole plan: to make a
response exploit a particular symmetry group, we only need the basis to carry that group. No solver
changes.

2.1 The wavevector q

A perturbation at wavevector q has the Bloch form δV(r) = δV_q(r) e^{iq·r} + c.c., with δV_q
lattice-periodic. It couples a state at k to states at k+q. Two regimes:

• q = 0 (electric field, strain, zone-center phonon): k+q = k, everything lives on one basis.
• q ≠ 0 (phonon dispersion): the response at k lives at k+q; the code already threads this
  through transfer_blochwave_equivalent_to_actual and stores δρ_q as a complex periodic
  array (compute_δρ: real_qzero = identity when q≠0).

2.2 Why the naive symmetry reduction is wrong for a response

If you build the response on a symmetry-reduced basis today and let compute_δρ/apply_χ0
symmetrize δρ/δV over the full crystal group, you get the wrong answer — usually zero.

The reason is the crux of this whole document. Symmetrizing a density n over a group G' is only
a valid full-BZ reconstruction if n is invariant under G'. The self-consistent response
density δρ^α (say, to an electric field along α) is not invariant under the full group: a
rotation S maps the field direction α to another direction, so it maps δρ^α to δρ^{Sα}, a
different response. Averaging δρ^α over the full group averages together responses to different
perturbations, and for a field the vector average is zero. (ABINIT calls this "over-symmetrization
of the density response".)

So a response perturbation is compatible only with the subgroup that leaves it invariant, and
that subgroup — not the full group — is what its basis must carry. Deriving that subgroup, and the
representation by which the perturbation transforms under it, is the content of §3.

---

3. The unifying abstraction

3.1 A perturbation carries (q, components, representation)

Everything a response needs to know about symmetry is how the perturbation transforms under the
crystal group. That is captured by three things, and nothing else:

abstract type Perturbation end

perturbation_qpoint(p)                    -> Vec3          # q  (0 for field / strain / Γ-phonon)
perturbation_components(p)                                 # the index set {a}
perturbation_representation(p, model, S)  -> D_S           # how S mixes the components {a}

• q determines the little group — the ops that map the perturbation's wavevector to itself:

  little_group(model, q) = filter(model.symmetries) do S
      is_approx_integer(transform_kpoint_coordinate(S, q) - q)     # θ·S·q ≡ q  (mod reciprocal)
  end

  Note this uses transform_kpoint_coordinate, so it naturally includes θ=-1 ops: a θ=-1 op is
  in little_group(q) iff -S·q ≡ q, i.e. its spatial part maps q ↔ -q. That is the time-
  reversal link between +q and -q, expressed as group membership (see §5.2).

• components {a} index the family of perturbations solved together: Cartesian directions
  1:3 for a field, (atom s, direction α) for a phonon (3·n_atoms of them), Voigt indices
  1:6 for strain.

• representation D(S) is the matrix by which S ∈ little_group(q) mixes the components. This
  is where the physics of each perturbation lives, and the only place it lives:

  perturbation	q	components a	D(S)
  electric field (ε∞)	0	direction, 3	vector rep W_cart[S]
  atomic displacement (phonon)	q	(s,α), 3N	small rep: atom-permutation ⊗ W_cart, with q-phase
  strain (elastic/piezo)	0	Voigt, 6	symmetric rank-2 rep sym(W_cart ⊗ W_cart)

  For θ=-1 ops, D(S) is antilinear (it acts after complex conjugation).

Two more traits round out the interface:

compute_δHψ(p, basis, ψ, ρ; ω=0)   # build the RHS δH·ψ per component (see §3.3)

ω (frequency, default 0) is reserved for dynamical/TDDFT response; it shifts the Sternheimer
denominator and is symmetry-inert, but wiring it now avoids a signature break later.

3.2 The basis stays dumb: it honours a group

The basis must not know what an ElectricField is. It only needs to be constructible with a
reduced symmetry group. The single new primitive:

PlaneWaveBasis(model; Ecut, kgrid, symmetries = model.symmetries)

The symmetries kwarg is intersected with (validated against) the model group, then run through the
existing rgrid/kgrid filters. This is literally "this basis breaks some of the model's symmetries;
reduce the group, then build k-points as usual." It feeds the existing k-point-reduction code
unchanged, and it retires the current work-around of building a whole new Model(model; symmetries=…).

A generic response helper ties the two together (in src/response/):

# Build the basis on which to solve the response to `p`, carrying transferred ground-state data.
function response_basis(scfres, p::Perturbation)
    model = scfres.basis.model
    G     = little_group(model, perturbation_qpoint(p))          # or a subgroup, see §3.4
    basis = PlaneWaveBasis(model; scfres.basis.Ecut, scfres.basis.kgrid, symmetries = G)
    full  = unfold_bz(scfres)                                    # ψ on the full BZ
    ψ     = transfer_blochwave(full.ψ, full.basis, basis)        # onto IBZ(G)
    # eigenvalues/occupations are symmetry-invariant: map them by k-coordinate
    (; basis, ψ, eigenvalues=…, occupation=…, unfold_map=…)
end

This is perturbation-agnostic: it consults only perturbation_qpoint (and, for the subgroup
variant, perturbation_representation). The dielectric and phonon drivers both call it.

3.3 Two orthogonal axes — do not conflate them

A recurring temptation is to unify "how the RHS δH·ψ is built" with "how symmetry is reduced".
They are independent and must stay so:

• Axis A — building the RHS δHψ. Per-perturbation, already solved, left alone:

  • field: the d/dk velocity operator (compute_δHψ_dk / shift_kpoints), because E·r is not a
    valid periodic term (modern theory of polarization);
  • phonon: ∂/∂(atomic position) of the local + nonlocal terms (compute_δHψ_αs);
  • strain: ∂/∂(lattice) (à la compute_stresses_cart).

  These are three genuinely different constructions; forcing them into one abstraction ("a
  perturbation is an infinitesimal term") fails precisely on the field, which has no term. Leave
  them as three functions.

• Axis B — the symmetry group and representation. This is (q, D(S)), entirely independent of
  how the RHS was built. This is what the plan unifies.

3.4 One knob: IBZ(little group) with a covariant symmetrization

Given the little group L(q) and the representation D(S), the response is computed on IBZ(L(q))
and the density response is reconstructed by a covariant symmetrization (§4). The components
couple only through that symmetrization step, because χ0 and the Hartree/xc kernel K are
block-diagonal in the component index a. Concretely the response driver is:

for each component a:  δρ^a = χ0[ δV^a_ext + K δρ^a ]      (per-component, reuse existing solvers)
covariantly symmetrize the family {δρ^a} over L(q)         (the one new coupling step)
iterate to self-consistency

There is one mechanism, parameterized by (q, {a}, D(S)):

• field: q=0, L(0) = full crystal group, D(S) = W_cart (3 coupled directions) →
  smallest possible k-set (IBZ(full group)).
• phonon: q, L(q) = little group, D(S) = displacement small rep (3N coupled) → IBZ(L(q)).

For the field there is also a cheaper "solve one direction at a time on its stabilizer subgroup"
option (scalar, no component coupling; this is what an earlier dielectric plan called "option B").
It does not generalize to phonons (a single atomic displacement has essentially no residual
symmetry), so this plan makes the covariant path primary and treats the scalar per-direction
path as an optional optimization, not the foundation.

---

4. The covariant response symmetrizer (the one genuinely new kernel)

Everything above reduces to a single new operation: symmetrize a family of complex response
densities {δρ_q^a} over L(q) with the representation D(S). This is not a new kernel from
scratch — it is accumulate_over_symmetries! generalized on three axes, with the phase factor
unchanged. Precisely (derivation in Appendix A):

(S ⋆_q δρ_q)^a (G)  =  cis2pi(-θ · G·τ) · Σ_b D(S)_{ab} · Θ_S[ δρ_q^b( θ·S⁻¹·(G − G_S) ) ]

    θ    = S.θ
    G_S  = transform_kpoint_coordinate(S, q) − q          (integer reciprocal vector; ≠0 only at zone boundary)
    Θ_S  = identity        if θ = +1
         = complex-conj    if θ = -1

and the full symmetrization is δρ_q ← (1/|L(q)|) Σ_{S∈L(q)} (S ⋆_q δρ_q).

Compared with the tested accumulate_over_symmetries!
(ρaccu(G) = Σ_S cis2pi(-G·τ) · ρin(S⁻¹G)), exactly three things change and the phase does not:

1. index θ·S⁻¹·(G − G_S) instead of S⁻¹·G — the θ (±q) linkage plus the zone-boundary
   G_S shift;
2. conjugate the fetched coefficient when θ = -1;
3. mix components by D(S) (antilinear when θ=-1), a small dense matmul outside the G-loop.

Two consequences worth stating loudly because they trip people up:

• Antiunitary ops are retained. symmetrize_ρ drops them (§1.4) because ground-state ρ is real.
  δρ_q is complex, and the θ=-1 ops are exactly what carry the -q information. So the
  response symmetrizer is a separate code path from symmetrize_ρ (or symmetrize_ρ gains
  q, D, and a keep_antiunitary=true flag).

• No new q-dependent phase in the spatial kernel. The density-symmetrizer phase stays
  cis2pi(-θ·G·τ), independent of q. Every q-dependent phase (the non-symmorphic e^{2πi q·w},
  the phonon small-representation atom-permutation phase) lives inside D(S), computed per
  perturbation with find_symmetry_preimage + W_cart + the q-translation phase (Gonze–Lee). This
  clean separation means the spatial kernel is a 3-line change to tested code, and the representation
  subtlety is confined to perturbation_representation, unit-testable in isolation.

At q=0, θ=+1, single component, D=1, this collapses exactly to
accumulate_over_symmetries! — the consistency anchor.

---

5. How this subsumes the current special cases

5.1 Dielectric (electric field, q=0)

Drop the symmetries=false requirement. Run the covariant driver with ElectricField
(q=0, D(S)=W_cart, 3 coupled directions) on IBZ(full group). The RHS is the existing
compute_δHψ_dk. Contract to ε∞[β,α]. The current per-direction scalar path becomes an optional
optimization.

5.2 Phonons (atomic displacement, q≠0) — and TRS for free

For each q, build response_basis(scfres, PhononPerturbation(q)) with symmetries = L(q), run the
covariant driver with the displacement small rep, contract to the dynamical matrix. This reduces
the k-point set per q, which is the headline win over the current symmetries=false code.

The phonon code's hand-rolled time-reversal trick ("one Sternheimer at +q instead of ±q",
src/postprocess/phonon.jl) becomes automatic: when time reversal is unbroken, L(q) contains
θ=-1 ops with -S·q ≡ q, and covariantly symmetrizing over L(q) reconstructs the -q
contribution from the +q solve via the Θ_S = conj branch of §4. A perturbation that breaks
conjugation (a magnon) simply has no θ=-1 ops in its little group, so more k-points survive and the
genuine two-solve cost reappears — through the same machinery, no special path. breaks_conjugation
(term level) and "does the perturbation keep the θ=-1 ops" (perturbation level) are the same idea
at the two levels of §1.5.

5.3 The RHS transfer at q≠0 (the one real implementation gotcha)

The TRS PR made find_equivalent_kpt error when a k-point is present only as a symmetry image
(it points you at unfold_bz). On an IBZ(L(q)) basis, the k+q point generally is only present
as an image, so transfer_blochwave_equivalent_to_actual must be extended to resolve k+q through
apply_symop (now θ-aware). This is bounded, well-identified work, and it is the same obstacle for
symmetry-reduced dielectric and phonons.

---

6. Higher and mixed derivatives (why the first-order object is the thing to lock down)

Every higher/mixed derivative is a contraction of first-order responses (the 2n+1 theorem,
Baroni RMP 2001 §IV). None of them needs a new response primitive:

• Born effective charges / IR ∂²E/∂u∂E: contract the displacement response with the field
  response (both at q=0, same IBZ(full group)). Output Z*_{sα,β} symmetrized by the
  (displacement ⊗ vector) representation.
• Electron–phonon g_{mn}(k,q) = ⟨ψ_{m,k+q}|δH^{u,qν}|ψ_{n,k}⟩: a matrix element of the phonon
  RHS. Same L(q).
• Anharmonic force constants ∂³E/∂u³ at (q₁,q₂,q₃), Σqᵢ=0: contract first-order responses
  at each qᵢ, on a common mesh.
• Raman / nonlinear optics ∂³E/∂u∂E², ∂³E/∂E³: contractions of first-order field/displacement
  responses.

Requirement banked now: the first-order response must be storable and unfoldable to the full
BZ (keep δψ, δρ, and the IBZ→full k-mapping), because cross-derivatives at different q live
on different reduced bases and must be brought to a common mesh (unfold_bz/transfer_blochwave).
The output tensor is then covariantly symmetrized once at the end, exactly like symmetrize_stresses.

---

7. Extension points (design for them, implement later)

The abstraction must be stated generally so these never force a redesign, even if v1 implements
only the scalar-density, unitary-spatial case.

• Non-collinear magnetism. SymOp will grow an SU(2) spinor matrix (the TRS PR's own TODO); the
  density becomes a spin matrix and Θ_S generalizes from conj to the spinor action. The
  symmetrizer's "physical representation" axis (currently trivial for a scalar density, Wigner-D for
  Hubbard) is where this plugs in. Nothing in the perturbation interface changes.
• DFT+U. Orthogonal to symmetry: +U adds a Hubbard block to the kernel K (needs the term's
  apply_kernel). The response occupation δn symmetrizes with the same Wigner-D representation
  symmetrize_hubbard_n already uses — another instance of "symmetrize a quantity by its
  representation".
• Meta-GGA. Adds δτ (kinetic-energy density), a scalar that symmetrizes like δρ. Fine once
  the τ-response is wired into the kernel (a term-level concern, not a symmetry one).
• Dynamical / TDDFT response χ(ω). ω shifts the Sternheimer denominator and is symmetry-inert;
  reserved in compute_δHψ/the driver signature from day one.
• Finite electric field (Berry-phase, nonlinear): not a Sternheimer response, so it never
  touches the driver — but it does break symmetry to the field stabilizer, which is handled purely
  by the §3.2 basis primitive. This is positive evidence that "basis carries a group" and
  "perturbation drives a response" must remain separate concepts.

---

8. Prior art (we are not inventing an architecture)

• ABINIT indexes every DFPT perturbation as (ipert, idir, qpt) — displacements, d/dk,
  electric field, strain — and "select[s] the k-point set for each perturbation using the symmetries
  that leave it invariant", explicitly warning about over-symmetrization of the density response.
  This is exactly §3.
• Quantum ESPRESSO (PHonon) treats the field as three independent directions, phonons with the
  small group of q, and builds Born charges / dielectric as contractions — exactly §5–§6.
• Baroni, de Gironcoli, Dal Corso, Giannozzi, RMP 73, 515 (2001): monochromatic perturbation at
  q (§II), the d/dk electric-field trick (§II.C), the 2n+1 theorem (§IV). This is the layering
  of §2–§6.
• Gonze & Lee, PRB 55, 10355 (1997): the phonon small-representation phases that live in D(S).

Conforming to this proven structure is the point: a wrong abstraction here costs a rewrite.

---

9. Proposed layering (summary)

Layer	What	Knows physics?
basis.symmetries (exists)	single source of truth for reduction + symmetrization	no
PlaneWaveBasis(…; symmetries=G) (new kwarg)	build a basis on a reduced group	no
breaks_symmetries / breaks_conjugation (exist)	coarse per-class term traits → model group	minimal
Perturbation + little_group + response_basis (new, src/response/)	perturbation → subgroup → basis+data	one method each
covariant response symmetrizer (new, generalize accumulate_over_symmetries!)	reconstruct full-BZ δρ_q from IBZ	no (takes D(S), q)
response driver (generalize solve_ΩplusK_split usage)	per-component χ0/K + covariant symmetrize	no
higher-derivative contractions (new)	Born/Raman/e-ph/anharmonic from stored first-order responses	yes

Invariants to bank (so nothing forces a later change):

1. Never S*k; always transform_kpoint_coordinate. Little groups and transfers are θ-aware.
2. ±q/TRS is θ, not a code path. Phonon single-solve = θ=-1 ops in L(q); magnon two-solve =
   their absence. One machinery.
3. Model group = coarse per-class term booleans; response subgroup = per-op perturbation geometry.
   Do not generalize term traits to per-op.
4. The covariant response symmetrizer keeps antiunitary ops (complex δρ_q), diverging from
   symmetrize_ρ's real-ρ shortcut; adds D(S) and the G_S shift; phase unchanged.
5. Non-collinear = SymOp gains an SU(2) field (already planned); the symmetrizer applies it. Not a
   response-layer invention.
6. First-order response is the storable/unfoldable atom; higher derivatives are 2n+1 contractions.

---

10. Suggested build order

1. Basis primitive. PlaneWaveBasis(model; symmetries=G) + validation. Retire the
   Model(model; symmetries=…) work-around. (Small, self-contained.)
2. Perturbation interface + ElectricField. q, components, D(S)=W_cart, and reuse
   compute_δHψ_dk for the RHS.
3. Covariant symmetrizer (§4), unitary + antiunitary, with D(S); the q≠0/G_S path guarded.
   Anchor test: at q=0, θ=+1, D=1, it must equal accumulate_over_symmetries! bit-for-bit.
4. response_basis + generic driver. Wire the dielectric onto it (§5.1); confirm it reproduces
   the symmetries=false ε∞ on silicon.
5. k+q transfer through symmetry images (§5.3): extend transfer_blochwave_equivalent_to_actual.
6. PhononPerturbation + displacement small rep (D(S), the Gonze–Lee phases). Refactor
   compute_dynmat onto response_basis; drop its symmetries=false and its hand-rolled TRS trick.
7. Higher derivatives (separate PRs): Born charges first (cheap, q=0).

11. Acceptance tests (the correctness gates)

• T1 — anchor. Covariant symmetrizer at q=0, θ=+1, D=1 reproduces accumulate_over_symmetries!.
• T2 — field. q=0, D(S)=W_cart covariant path reproduces the symmetries=false dielectric
  ε∞ on silicon; tensor isotropic.
• T3 — equivariance (the direct lemma check). Random δV_q: χ0 δV_q on the full basis equals
  χ0 on IBZ(L(q)) + covariant symmetrize, for
  (a) interior q, G_S=0, on Si (has inversion ⇒ θ=-1 ops present);
  (b) GaAs (no inversion ⇒ tests θ handling with fewer antiunitary ops);
  (c) a zone-boundary q (G_S ≠ 0).
• T4 — non-symmorphic. A structure with w≠0 to exercise the D(S) q-phase.
• T5 — end-to-end phonon. Dynamical matrix at a q with nontrivial L(q), symmetry-reduced vs
  symmetries=false, agree to solver tolerance.

T3(c) and T4 are the ones to watch: they exercise the only two ingredients not already covered by
tested q=0 code (the G_S shift; the D(S) q-phase).

---

Appendix A. Derivation of the covariant symmetrizer (§4)

Claim. With δρ_q the full-BZ response density (insulator, δf=0) stored as in compute_δρ,

δρ_q^full = (1/|L(q)|) Σ_{S∈L(q)} S ⋆_q δρ_q^IBZ,

with S ⋆_q as in §4.

Conventions (from code). Build everything from the tested periodic-part law of apply_symop:

û_{θSk}(G) = cis2pi(-θ·G·τ) · Θ_S[ û_k(θ·S⁻¹·G) ].

δρ_q is stored (insulator) as, in real space,

δρ_q(r) = Σ_k w_k Σ_n 2 f_nk · conj(u_nk(r)) · δu_{n,k+q}(r),

whose Fourier coefficient is the cross-correlation
δρ_q(G) = Σ_k w_k Σ_n 2 f_nk Σ_{G''} conj(û_nk(G'')) · δû_{n,k+q}(G''+G).

Step 1 — χ0 is equivariant under L(q). The non-interacting response solves, per occupied
(n,k),

(H_{k+q} − ε_nk) δu_{n,k+q} = −P_c^{k+q} (δV_q u_nk).

U_S intertwines H_k with H_{θSk} (ground-state symmetry, both θ); ε is real; P_c is
spectrally defined, hence commutes with U_S. For S∈L(q) the image of k+q is
θS k + θS q = θS k + q — again a +q problem. Applying U_S:

• θ=+1: linearly maps the (n,k) solve to the (n, θSk) solve of the transformed perturbation;
• θ=-1: antilinearly — conjugates the equation. As H−ε has real spectrum and δV_q → conj at
  the mirrored momentum, the conjugated +q solve at the image reproduces the -q data.

Hence χ0 ∘ U_S = U_S ∘ χ0 on q-perturbations, i.e. χ0 commutes with S ⋆_q. Since the
Hartree/xc kernel K also commutes with symmetry and is block-diagonal in the component index, the
full self-consistent response is equivariant too, and the components couple only through the
symmetrization.

Step 2 — transformation of the stored δρ_q. Relabel the full-BZ k-sum by k → θSk and
substitute the û law into the cross-correlation. Phases: conj(û_nk) contributes
cis2pi(+θ G''·τ), δû_{n,k+q} contributes cis2pi(-θ(G''+G)·τ); the G'' pieces cancel, leaving
cis2pi(-θ·G·τ) — q-independent, identical to the q=0 kernel. The index acquires θ·S⁻¹, with
Θ_S = conj for θ=-1 (since conj(conj(û))·conj(δû) = û·conj(δû), the conjugate of the original
correlation at the mirrored G). The zone-boundary shift G_S enters via the k+q sphere's wrap
into [-1/2,1/2). This is exactly S ⋆_q of §4. All q-dependent phases (non-symmorphic τ+q,
atom permutation) are carried by D(S) in the component mixing, not the spatial kernel.

Step 3 — assembly. Partition the full BZ into L(q)-orbits; equivariance turns each orbit's sum
into S ⋆_q applied to the IBZ representative's contribution; averaging over L(q) with the usual
weight bookkeeping (unchanged from the scalar case) gives the claim. At q=0, θ=+1, one component,
D=1, G_S=0, it reduces to accumulate_over_symmetries!. ∎

Assumptions / failure modes (each must be respected or guarded):

• Insulator (δf=0). The metallic Fermi-surface term is real and scalar-symmetrizing; add later,
  gate on a gap for now.
• Degeneracy/gauge. δρ_q sums over degenerate n, so it is gauge-invariant; safe.
• P_c commutation. Holds because the occupied/virtual split is symmetry-invariant.
• Zone-boundary q (G_S≠0). Essential; the generalized lowpass_for_symmetry! must drop
  wrapped-out G (where index_G_vectors returns nothing), exactly as at q=0.
• The real-ρ antiunitary-drop optimization must NOT be reused for the complex δρ_q.
• Non-collinear needs the SymOp SU(2) field; Θ_S generalizes accordingly.

[jaemolihm]

I just skimmed but couldn't resist to reply 😂

I think (q, components, representation) is the right abstraction. This is what I implemented in QE to make the DFPT kernel fully agnostic to the perturbation type.
(The MODULE that stores the representation is https://gitlab.com/QEF/q-e/-/blob/develop/LR_Modules/lrcom.f90#L113-146)

For those who are interested, I attached the slides where I explain this refactoring in a QE developers meeting
JaeMoLihm_QE_dev_meeting.pdf

An alternative way to implement symmetry in DFPT is the Abinit way, where "each perturbation has its own list of irreducible k-points" (https://docs.abinit.org/faq/dfpt). There, the number of perturbations can be reduced to the irreducible ones (e.g. for E field in cubic materials you need to solve for x only, not x,y,z as in QE way), but this limits the k point reduction to the symmetry of that perturbation (i.e. symmetries that keep the perturbation unchaged). One should do some math, but I think this is slightly less efficient than the QE way if one has only a few perturbations.

[antoine-levitt]

Thanks, super helpful as always!! I wish I had /model jaemolihm :-) Yes, I had hoped to avoid treating several perturbations at the same time but it looks like it's cleaner that way... How do you mix the deltarhos in this case? Do you just mix them independently, or put them together in a big array?

[Technici4n]

What I had in mind was a partial unfold for each perturbation, i.e. for each perturbation (each DFPT rhs), select a subset of the symmetries and create a new basis that only has these symmetries (i.e. a "partial unfold_bz"). Then do the DFPT and post-processing in that basis. The only non-trivial thing to do is to implement the partial unfold logic. Is there anything more to it than this? What is the purpose of combining multiple perturbations in the same DFPT call (i.e. why is lr_npert not always 1 in QE for example)?

(As a side note, I find it hard to identify the signal in the issue text :( )

[antoine-levitt]

Yes there should be a partial unfold_bz, but there's more than that.

What is the purpose of combining multiple perturbations in the same DFPT call

Symmetries couple different perturbations together: (S ⋆_q δρ_q)^a (G) = cis2pi(-θ · G·τ) · Σ_b D(S)_{ab} · Θ_S[ δρ_q^b( θ·S⁻¹·(G − G_S) ) ]. Also we sort of already half-bake the +q/-q coupling by TRS in a way that never felt right to me (for instance it probably wouldn't work for the computation of magnons), I was kind of hoping to generalize this mechanism in a way that feels less hackish (my main motivation for this is to clean up by finding the right abstractions, speed is a bonus).

(As a side note, I find it hard to identify the signal in the issue text :( )

Well of course it's AI slop, but I think there's actual signal there, and I don't understand this well enough yet to come with something better. Honestly, I find reading this marginally less annoying than actually reading physics papers :-p there's plenty of special-case derivations, but it's very hard to find a source discussing this in a general way (but maybe I haven't looked enough). I thought posting this publicly would be a good basis for discussion with interested parties. (and it worked, it summoned @jaemolihm!)

[jaemolihm]

😂
The partial unfold is the abinit way I mentioned above. When treating multiple perturbations that form irrep together, you don't need any unfolding at all. (Or the calculation for different perturbations is implicitly doing the partial unfold.)

For mixing, we concatenate ans mix in a large array, but in principle since drho from one perturbation determines the rest via symmetry, I think mixing one and reconstructing the rest with symmetry could be simpler.

[Technici4n]

for E field in cubic materials you need to solve for x only, not x,y,z as in QE way

Ah so QE will do all 3 perturbations at once (so "3x the work") and keep the IBZ from the full symmetries (or little group for phonons). Whereas ABINIT can do only 1 perturbation, but then has to do it with a larger IBZ. Is QE really able to for example compute electric field perturbations without removing any symmetry? I think I almost understand...

I never thought about the QE way, I always assumed one would do it the ABINIT way. 🤯

(BTW, thanks a lot for chiming in, Jae-Mo! I echo Antoine in saying that your input is extremely useful :)).

[Technici4n]

I understand it now after some Claude prompting and more thinking. It seems that the QE way is the right thing to do for standard properties where you want the full set anyway. Actually, I never understood this part of the 2009 QE paper before but now it makes sense:

Notice how it has to mix the density responses across different perturbations 😉

It seems that this issue is suggesting a rather elaborate interface and that we could get away with something much simpler. I would suggest starting a prototype with the elastic constants since we already have a helper for those and it doesn't need to involve q-points.

Maybe one thing that will be annoying with the QE approach is the symmetrization for hybrid DFPT. But I suppose we'll cross that bridge once we get there.

[antoine-levitt]

Hm I just remembered #1251 where I argued that symmetrization is not that important for phonons. So maybe that + abinit way is the simplest... still have to think about the +q/-q business.

[mfherbst]

So maybe that

@antoine-levitt To what precisely are your referring to with 'that' in your comment ?

[antoine-levitt]

#1251: computation on the full BZ, clever selection of the perturbations, and reconstruction of the full set from that.