DiagonalizeSupercell mutates the caller's dynamical matrix (in-place Re() at q = -q + G)

[mesonepigreco]

Summary

Phonons.DiagonalizeSupercell (and DiagonalizeSupercell_slow) silently modify the dynamical matrix they are called on. For every q = -q + G point the code enforces reality by writing the real part back into the object:

# Enforce reality to avoid complex polarization vectors
self.dynmats[iq] = re_part

The stripped imaginary part is tiny (~1e-8) but nonzero, so a byte-exact consumer of the matrix observes a changed input after the call.

Concrete failure

Breaks tdscha's distributed interpolated (atom-Fourier) backend:

1. tdscha builds the fine-mesh harmonic interpolation and a byte-exact fingerprint of the converged dynamical matrix before the master/worker MPI split.
2. sscha.Ensemble.update_weights(final_dyn, T) then calls final_dyn.DiagonalizeSupercell(...), which zeros the imaginary part of the non-Gamma blocks of the caller's matrix in place.
3. The Lanczos constructor validates the precomputed interpolation against the now-modified matrix and raises:

ValueError: the precomputed harmonic interpolation was not built from this
dynamical matrix and these interpolation settings; using it would contract
the ensemble in a mode basis that does not belong to it

Because the mutation is ~1e-8, any value-based (not byte-based) fingerprint hides it; only byte-exact digests expose the bug. Any other code that compares a dynamical matrix before/after a DiagonalizeSupercell/update_weights round trip is affected the same way.

Root cause

DyagDinQ(iq) reads self.dynmats[iq] internally, so the immediately following diagonalization must see the enforced-real matrix. The write-back was intended as a local temporary but is never undone, leaking a side effect into every caller.

Fix (patch attached as PR / in commit)

Keep the temporary write so DyagDinQ still sees the real matrix, and restore the caller's original matrix immediately after the diagonalization, in both the fast and slow implementations:

is_minus_q = False
dynmat_at_iq = None
if Methods.get_min_dist_into_cell(bg, q, -q) < 1e-6:
    is_minus_q = True
    dynmat_at_iq = self.dynmats[iq]
    re_part = np.real(dynmat_at_iq)
    ...
    self.dynmats[iq] = re_part   # temporary, seen by DyagDinQ below
...
if is_minus_q:
    # Undo the temporary reality enforcement: DiagonalizeSupercell must not
    # mutate the dynamical matrices it was given.
    self.dynmats[iq] = dynmat_at_iq

Numerics are unchanged: the enforced-real matrix is still diagonalized, and the returned frequencies/polarizations are bit-identical.

Verification

• Reproducer: snapshot dyn.dynmats before/after dyn.DiagonalizeSupercell(return_qmodes=True); before the fix the non-Gamma blocks differ by max |Δ| ≈ 3.4e-8, after the fix they are identical.
• tdscha Spectroscopy(EnsembleSource, backend="atom_fourier") on a 100 000-configuration hydrogen ensemble (fine mesh 8x6x8) now builds the engine and runs instead of raising the interpolation-mismatch error.

Environment

• cellconstructor 1.6.2 (installed) / master 1.7.0
• Python 3.14, numpy 2.4.6
• tdscha feature/symmetry-spectroscopy