QSpace: fix the Schur fill between repeated irreps, and make the q-space path linear in memory

[SorBalda]

QSpace: fix the Schur fill between repeated irreps, and make the q-space path linear in memory

Three commits, reviewable independently. The first is a correctness fix and does
not depend on anything else.

Part of a three-repository series:

• CellConstructor: Phonons: q-only supercell diagonalization, in-place sort, LO-TO and AdjustQStar fixes CellConstructor#126 (provides DiagonalizeSupercell(q_only=True))
• python-sscha: Ensemble: opt-in linear-memory q-space mode, plus four fixes on the Fourier path python-sscha#428 (provides Ensemble(qspace_light=True))
• tdscha: this PR

This PR is safe on its own: with unpatched companion packages it detects
their absence and falls back with an explicit warning, rather than failing. But
the memory gain only materializes when all three are in, since the two
primitives live in the other repositories.

---

1. QSpaceHessian: the Schur fill between repeated irreps is wrong

Filling G on a degenerate block used the scalar shortcut c*I for the cross
block between two blocks as well. By Schur's lemma the diagonal block of a
d-dimensional irrep copy is c*I in any orthonormal basis of that copy, but
the cross block between two copies of the same irrep is c*U, with a unitary
intertwiner U that is not the identity: the blocks come from degeneracy in
frequency, and eigh picks an arbitrary basis inside each degenerate subspace,
so the two bases are unrelated. The shortcut is exact only where the coupling
vanishes, i.e. between distinct irreps.

The resulting error is symmetric by construction, so it survives the final
symmetrization and leaves no diagnostic.

_adaptive_schur_fill measures the coupling on the already-solved
representative columns, groups the coupled blocks with union-find, and solves
those groups column by column exactly. Distinct irreps keep the scalar path, so
nothing changes where the shortcut was valid, and a false positive only costs
solves. Cost: every block of a coupled group is solved in full, so a group of
k blocks of dimension d costs k(d-1) extra column solves — the earlier
wording said d-1 per group, which was wrong.

Measured on a cubic LaAlO3 supercell (5 atoms/cell, 27 q points, 405 modes,
200 configurations), against the exact reference obtained with
use_mode_symmetry=False, which solves every non-acoustic band as its own
column and therefore makes no Schur assumption at all. All three Hessians come
from the same ensemble, at tol=1e-10:

| | max |dPhi| | relative |
|---|---|---|
| scalar shortcut (current) | 3.3e-02 | 8.4e-02 |
| adaptive fill (this PR) | 1.5e-09 | 4.0e-09 |

The shortcut is off by 8% on this system; the adaptive fill matches the exact
solution to nine digits. The detection fires at 3 of the 27 q points. The size
of the error depends on the ensemble it is measured on -- the same comparison
at 200 vs 20 configurations gives 8% and 20% -- but its order of magnitude does
not, and it is always many orders above the adaptive fill.

This also means that turning mode symmetry on currently degrades the result,
which matches what we had been seeing empirically before finding the cause.

tests/test_qspace/test_schur_fill.py drives _adaptive_schur_fill with an
exact solver, so it needs neither Julia nor ensemble data.

What is still missing, and where we would like guidance. There is no
physics-level regression for this in CI. The LaAlO3 comparison above is a
manual three-way run against an ensemble that is not in the test data, and the
unit tests, while exact, are synthetic. A small committed ensemble exercising a
repeated irrep would make it permanently checkable — we are happy to add one if
you can point at a suitable system.

2. tdscha_qspace.jl: preallocate the symmetry-rotation buffers

The four q-space kernels allocated x_rot and y_rot inside the
(configuration, symmetry) loop, producing two fresh ComplexF64 vectors of
length n_q*n_bands per iteration for the GC. They are overwritten every
iteration and never outlive it, so they are hoisted out and filled in place with
mul!.

This makes the buffers shared across the loop: parallelizing it would require
making them thread-local. There is no @threads in this file, and a note sits
next to each buffer.

3. QSpaceLanczos: linear-memory path, with explicit capability fallbacks

The q-space Lanczos allocated the dense (3N, 3N) supercell polarization matrix
even though it only ever uses the per-q eigenmodes. It now asks for
DiagonalizeSupercell(q_only=True) and builds the ensemble with
Ensemble(qspace_light=True).

Both live in the companion packages, so their absence is detected rather than
assumed, by two signature probes:

• no q_only → fall back to return_qmodes with a warning. Same numbers, but
  the dense matrix is allocated, so the path is no longer linear.
• no Ensemble(qspace_light=) → the distributed loader builds a standard
  ensemble. This must be checked before building the light ensemble, because
  Ensemble.__setattr__ calls DiagonalizeSupercell(q_only=True) as soon as the
  flag is set: otherwise the failure happens there, before this class exists.

__JULIA_EXT__ cannot carry either decision: it only means "juliacall is
importable", which says nothing about the runtime or about the companion APIs.
load_distributed_tdscha also checks up front, before any I/O, that the ensemble
exposes the q-space cache API — otherwise the failure came only after load_bin
and a real-space update_weights had already been paid for.

Other fixes in this commit:

• N_eff is kept as a float. int() truncated the weighted count of each
  MPI slice, while Julia normalizes by the exact sum(rho).
• units="hartree" raises NotImplementedError on the q-space path instead of
  silently rescaling: the Fourier caches are built with Ry/Angstrom factors and
  Ensemble.convert_units does not update them.
• the symmetrization refuses to build a wrong matrix when a symmetry operation
  does not map the atoms onto themselves within tolerance; point-group
  operations that do not preserve an anisotropic q grid are skipped with a
  warning (they form a subgroup, so averaging over the rest is still a valid
  projector) instead of crashing.
• the ensemble's cache coherence marker is no longer written from here.

---

Testing. tests/test_qspace/ passes (26 passed, 7 skipped). The q-space and
real-space free-energy Hessians agree to 3.6e-5 cm-1 on LaAlO3 3x3x3 (405
modes, 200 configurations, same ensemble). test_neff_cast.py moves from work/ into tests/test_qspace/.

[mesonepigreco]

Review pass 2 — code hygiene / dependency contract (follow-up to the Schur-fill review).

The Schur analysis in commit 3 is solid (verified independently: the old shortcut is off by ~20% on a repeated-irrep pair, the adaptive fill is exact). This review is about the integration code in commit 3 (QSpaceLanczos.py), which currently works around the companion packages with runtime introspection instead of declaring the contract up front. Three points, all inline:

1. Pin the versions in requirements.txt / pyproject.toml and delete the inspect.signature capability probes (_cc_has_q_only, _ensemble_supports_light, line 67-88) — probing the internal API of CellConstructor/python-sscha is fragile and makes the fallback paths dead code that CI can never exercise. The return_qmodes=True dense fallback in __init__ and the standard-ensemble fallback in load_distributed_hessian should go away with them.
2. Top-module imports only: import inspect as _inspect is mid-file and import sscha.Ensemble is inside a function (line 81).
3. No hasattr (lines 1872-1873): the refresh_qspace_caches_from_real_space / _refresh_qspace_caches_from_real_space probe is a third, mutually inconsistent way of detecting the same capability; with the pinned version it is simply removable.

Both requirements.txt and pyproject.toml currently list bare cellconstructor / python-sscha (no minimum), which is what made these probes feel necessary — fix it at the source.

[mesonepigreco]

Review pass 3 — same import hygiene rule applied to the new test (test_neff_cast.py): top-module imports only. Two inline comments on run_mpi and _load_dyns.

[mesonepigreco]

Good, there are some edge cases for the new Shurr implementation that can be improved (see comments within the file changes, not sure it is good). Then there are multiple things that needs to be changed: the code use inspect which is very hacky. This should be removed and replaced with a hard constraint on cellconstructor and python-sscha versions (1.7, i.e., the next version). Also, imports must be module level, while tests are using it function level, which also feels hacky.

[SorBalda]

Thanks — the single-block gap is real and I had missed it. Pushed in e74a90c, together with the import cleanup.

Reproduced your leakage numbers exactly: 0.0000 for a clean copy, 0.7000 for α⊕α, 1.396 for E+A (bound |c_E−c_A|/2 = 1.400). On a lone reducible block the old code gives err = 7.60e-01, the fix 9.16e-16.

Two changes to the patch, both needed. I split the detection into two passes instead of doing it in one loop:

1. Order dependence. The relaxed dimension shortcut reads full_solve, which the single-loop version populates as it goes: coup_self(B) is only computed when the outer loop reaches i=j, so a pair (A, B) with different dimensions is skipped if B has not been visited yet. Reversing the block order changes the result. Regression test test_detection_is_order_independent fails on the single-loop version and passes with two passes.
2. Accumulated scale. scale = max(scale, ||xB||) inside the j loop never decreases, so the threshold is monotonic within a fixed i. Since ||column of G|| ~ 1/w², one soft-mode spectator sitting between two coupled blocks raises the threshold by orders of magnitude for every pair examined after it, and can only lose couplings. I recompute it per pair. I could not build a case where this actually bites in practice, so treat it as a structural fix rather than a measured bug.

Three regression tests added (test_schur_fill.py); two fail on the previous code.

Cost formula: you are right, and the PR body is fixed — k(d−1) for a group of k blocks, not d−1. The commit message already said "all their columns".

One limit worth recording, which neither criterion covers: if eigh happens to return an already symmetry-adapted basis, the intra-block leakage is exactly 0.0 while the two Schur constants still differ, and the fill is wrong (~78% on a synthetic case). Detecting it needs the columns the shortcut exists to avoid solving. It is documented in the docstring. A cheap upper bound would be to force the full solve for len(block) > 3 (no crystallographic irrep exceeds d=3), but on a Γ-supercell that is nearly all blocks — on LaAlO3 3×3×3 as a single cell, 390/405 modes — while on a primitive cell with an irreducible q grid it is free. Happy to add it behind a flag if you prefer correctness over the shortcut there.

On pinning the versions. I'd like to, but it isn't writable yet: no released version exposes the API. PyPI cellconstructor 1.6.2 and python-sscha 1.6.1, and both master branches, contain zero occurrences of q_only / qspace_light — the features only exist in CellConstructor #126 and python-sscha #428, both still open. A pin to an unreleased version would make tdscha uninstallable.

There is also a detail that cuts the other way: python-testsuite.yml (lines 48-66) git clones both companions from master and installs them directly, so it never reads requirements.txt. In CI _CC_HAS_Q_ONLY is therefore False, which means the fallback is currently the only branch CI can take — and it is numerically exact (max|Δw_q| = 0.0 between q_only=True and return_qmodes=True on a (72,4) grid). Dropping it today would give exactly the TypeError: DiagonalizeSupercell() got an unexpected keyword argument 'q_only' that broke CI on #428.

Proposal: keep the probes for now with the note that is in the code, and open a follow-up that pins the versions and deletes probes and fallbacks the day #126 ships. If you'd rather have the pin now and merge the three PRs together, that works too — just say which order you want them in.

Separately, requirements.txt and pyproject.toml disagree (spglib and mpi4py missing from one, and julia is listed while the code uses juliacall). Happy to fix that in its own PR since it is independent of this one.

[SorBalda]

Corrections to my previous comment, found by re-reviewing my own tests before you got to them. Pushed in 7e0753b.

1. I was wrong to downgrade the scale point to "structural". I had built and measured a case where it bites, and then broke it while moving it into the committed test: the spectator block got a small constant instead of a large one, so its column norm was the smallest of the three, the running maximum never moved and the test passed on every variant — it guarded nothing.

Corrected, with the spectator's column norm ~10³ times the coupled pair:

                        single-block   soft-mode   order
old (before this PR)    FAIL           pass        FAIL
single-loop detection   pass           FAIL        FAIL
two-pass detection      pass           pass        pass

So the single-loop version loses a real 3e-4 coupling entirely (err = 3.9e-04 against 2e-13). Both problems are measured, not one. Caveat: with bands sorted by increasing frequency the largest-norm block tends to come first, which reduces — but does not remove — the occasions where the running maximum bites on real data. Measured on a synthetic case.

2. d = 3 is not a bound, and the scope was too broad. "No crystallographic irrep exceeds d=3" holds for the 32 point groups, hence at Γ and for symmorphic little groups. At BZ-boundary q of non-symmorphic groups the small representations are projective and go higher — verified d=6 at R of Pm-3n (#223). And len(block) > 3 is a cheap trigger, not an upper bound: the zero-leakage case I mentioned has d=2, so it would not be caught anyway.

3. 1.396 is not an exact reproduction. 0.0000 and 0.7000 are; the E+A leakage depends on the arbitrary basis eigh returns. The deterministic quantity is the bound |c_E − c_A| / 2 = 1.400, which matches your "~1.4".

4. This PR changes a number under MPI, and the new one is the correct one. test_distributed_hessian started failing with AssertionError: Negative eigenvalues in Hessian. It builds with use_symmetries=False but left use_mode_symmetry=True, and it passed only because the scalar fill fabricated a degeneracy:

this PR,     use_mode_symmetry=True   evals = [-4.078e-06, -3.239e-07, 0, 0, 0, 7.17e-08]
this PR,     use_mode_symmetry=False  evals = [-4.078e-06, -3.239e-07, ...]   <- identical
before,      use_mode_symmetry=True   evals = [0, 0, 0, 1.332e-07 x3]         <- fabricated

The corrected fill reproduces the no-symmetry ground truth exactly. For a 10-configuration ensemble with symmetries off, the two lowest eigenvalues really are slightly negative, so I turned the mode symmetry off in that test to match the rest of the object and replaced the positivity assert with finiteness plus a loose bound, with the reason in a comment. mpirun -np 2 pytest tests/test_qspace/test_distributed.py is 6 passed.

5. Two more limits of the criterion, now documented in the docstring: the coupling between a singleton mode and an accidentally degenerate block is zeroed without ever being measured (Schur forbids coupling between different irreps, not different frequencies), and min(50*tol, 1e-5) is a floor on sensitivity. I have a patch for the first one that costs no extra solves; happy to add it here or keep it for a follow-up, whichever you prefer.

Unchanged and still correct: the 7.60e-01 → 9.16e-16 figures, the k(d−1) cost, and everything in the pinning discussion.