Add symmetry-aware Raman and IR spectroscopy - #36
Open
mesonepigreco wants to merge 23 commits into
Open
Conversation
Run the full q-space Lanczos on a fine uniform q-mesh not commensurate
with the ensemble supercell, keeping the stochastic ensemble on the
coarse grid and never forming 3rd/4th order force-constant tensors.
Design document: Interpolation_plan.md.
- Modules/QSpaceInterpolation.py: QSpaceLanczosInterp(QSpaceLanczos)
* per-configuration NUDFT of displacements and force residuals at the
fine q-points (plain window; exact at commensurate q)
* f_Y field pre-filtering: strips the phonon-propagator dressing of
every displacement leg on the coarse grid, so the interpolated
correlations decay with the range of the anharmonic force constants
(essential: without it the renormalization is systematically ~0.5x)
* interpolated dynamical matrix via ForceTensor.Tensor2 centering +
Apply_ASR, TRI gauge e(-q)=conj(e(q)), commensurate points reused
* acoustic-sum-rule zero-mode projection of the fields (mass-metric
translations only: optical components untouched)
* O(N_f) pair map and hash q-index lookup
- tdscha_qspace.jl: scale3=sqrt(N_c/N_f), scale4=N_c/N_f vertex
renormalization (mode-space d3 ~ N^-1/2, d4 ~ N^-1) and `prefiltered`
flag; defaults preserve the previous behavior exactly.
- tests/test_interpolation/: anharmonic diatomic-chain toy (topology-
based, identical physics for every supercell length); V0 identity at
machine precision, Tensor2 sign-convention pin, TRI/ASR/pair-map unit
tests, and physics validation of interpolated vs direct fine-supercell
renormalizations (incl. non-TRI q and a scale-factor necessity test).
Measured on the toy (coarse L=3 -> fine L=6, N=4000): ~10% noise floor,
~30% plain-window aggregate error on the anharmonic renormalization, no
systematic factor. Next: designed multitaper windows (plan section 5.2).
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…l xfail Fix the last failing test: - The strict-xfail force-Parseval "negative control" was NOT a hand-built- dyn artifact: ens.generate() zeroes forces_qspace and assigning ens.forces afterwards never refreshes it (QSpaceLanczos skips ens.init() when u_disps_qspace exists). _build_nontri.make_ensemble now calls ens.init(); the test passes as a positive regression guard. M3 (designed multitaper windows, plan sections 5.2-5.5): - Window design toolbox: per-dimension ALS fit of slot-resolved window triples to the minimal-image (zero-padded centering) kernel under a hard partition-of-unity constraint (exact commensurate limit). Numerically exact designs for L=3 (K=3) and L=4 (K=2). - Slot-resolved Julia kernel (z/w/v field slots, compute_d3/compute_d4 flags); w<->v orientation averaging preserves the block transpose symmetry and Hermiticity. D4 terms keep the plain window. - Origin averaging by data permutation (rolling the configuration), NOT by cyclic window shifts (which wrap the support with inconsistent phases). - ASR zero-mode projection restricted to (near-)uniform windows: for oscillating windows it is a rank-one modification of the effective window and biases the kernel (measured ~30% on the toy). - Measured (chain toy, Lc=3->Lf=6, N=4000, vs direct fine ensembles): noise floor 0.103, plain window 0.299, windows 0.158, windows+3 origins 0.122 -- designed windows reach the statistical noise floor. M4 (performance): - Batched BLAS-3 slot kernel (sparse x dense rotations, GEMM alpha1 intermediates and pair accumulations, chunked over configs at fixed symmetry); equivalent to the scalar kernel to 1e-20. Interp plain path routed through it. Near-linear wall-time scaling measured for N_f = 8..128. Tests: tests/test_interpolation/test_windows.py (designer exactness, partition of unity, batched==scalar, trivial-window==plain identity, origin-shift sanity, and the decisive physics comparison windows-vs-plain against a direct fine-supercell ensemble). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
q-mesh interpolation report/interpolation/main.tex (+ compiled PDF, figures, data, scripts): - full derivation of the stochastic estimator, the effective windowed interpolation kernel (with the partition-of-unity/commensurate-limit proof and the image-decomposition appendix), the equivalence with centered-tensor Fourier interpolation, the N_c->N_f vertex rescaling, the Gaussian-IBP pre-filter identity, and the acoustic-sum-rule treatment; - benchmarks: machine-precision back-compatibility, mode-resolved renormalization vs direct fine-supercell ensembles (noise floor 0.103, plain 0.299, minimal-image windows + 3 origins 0.122), near-linear wall-time scaling N_f = 8..128; - example application: one-phonon spectral functions at interpolated wavevectors from a coarse (1,1,4) ensemble vs the standard d3 dynamic bubble of cellconstructor.Spectral (exact toy Phi3, Tensor3 Center+Apply_ASR, same k-mesh/smearing): peak positions agree within one energy-grid step (~2.3 cm^-1) on 15-40 cm^-1 anharmonic shifts. Interpolation_plan.md: status updated through M3/M4 with the deviations found during implementation (windowed-ASR bias, origin-shift data permutation). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…gle constrained fit The acoustic sum rule of the interpolated D3 vertex maps onto linear image-sum constancy constraints on the interpolation kernel. One-period windows cannot satisfy them (uniform class sums <=> plain window: the no-go behind the M3 field-projection bias). Doubled-support (2L) windows parameterized with exactly uniform class sums, w(d+L) = c - h(d), make ASR on all three vertex legs structural, collapse the partition of unity to one scalar KKT constraint (machine-precision commensurate limit, verified 1e-22 at field level), and are fitted to the closed-form ASR projection of the minimal-image kernel - where ForceTensor must alternate Apply_ASR with re-symmetrization iteratively, all constraints here hold simultaneously in one constrained ALS fit. Runtime: per-q complex effective weights (second period = extra Bloch phase across the supercell); no Julia changes. Optional decay-weighted metric = analogue of Apply_ASR(power). Validation (new figs/data + 8 tests in test_asr_windows.py): - deterministic kernel leak: minimal-image O(1) at finite q1 (vanishes only when both pair legs -> Gamma together), asr design 1e-15; - pairwise potentials are structurally immune to the leak (violating weights touch only three-cell-spanning entries): added a three-body g3b*s1^2*s2 term to the toy chain to expose it; - stochastic zone-boundary probe: minimal-image acoustic vertex plateaus (7.6x plain at q2=1/48), asr tracks the plain decay; - response level: under-resolved coarse mesh (spread-2 tensor at the L=3 WS tie boundary) breaks ALL oscillating designs -> use plain; at resolvable range (Lc=5) plain fails/halves the large renorms while minimal-image and asr track the direct reference; - pairwise-toy precision kept: asr+3origins 0.150 vs minimal-image 0.122, floor 0.103, plain 0.299. Report: new Sec. 6 (kernel-level ASR conditions, no-go theorem, one-shot construction, leak-vs-resolution), asr_kernel/asr_leak figures, updated abstract, benchmark table rows; fixed mangled sys.path lines in the benchmark scripts. Full suite: 55 passed, 7 skipped. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… + tests D4 is interpolated (off-grid + scale4=N_c/N_f) on a plain-window pass, never the designed/ASR windows. Adds an end-to-end report section validating this on a D4-dominated regime of the anharmonic chain: interp-from-coarse tracks the direct fine-mesh spectral function to <=2 cm-1 while dropping D4 is ~10x worse. Physics note: a purely quartic model renormalizes exactly zero (the D4 tadpole is resummed into w_SSCHA), so D4 is exercised via a mixed model. - report §"Fourth order": text, table, figs/d4_spectral.pdf, bench script - tests/test_interpolation/test_d4_interp.py: renorm-level interp-vs-direct (+ guard that dropping D4 breaks it), operator-level exactness (scale4 linear, bit-identical across window designs), pure-quartic degeneracy guard Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
- atomic_delta q-space ASR projector: subtracts uniform acoustic component from windowed Bloch fields with Tensor3-style sinc factor f_N(q) - SnTe 2x2x2->4x4x4 benchmark: atomic_delta TO peak = 37.52 cm-1 (direct 4x4x4 = 37.56, tensor-D3 oracle = 37.65) - ASR validation: near-Gamma acoustic D3 vertex reduced by ~14,000x (3.4e-17 vs 5.0e-13 for raw atomic) at |q|=0.066 A^-1 - Updated report/interpolation/main.tex: SnTe spectral figure, ASR figure, updated benchmark table with atomic_delta production results - Added figures: snte_spectral.pdf, snte_zoom.pdf, snte_asr.pdf - Added scripts/snte_figs.py for reproducible figure generation - Updated Implementation_plan.md with full production validation log
the interpolation strategy is inherently bugged and wrong.
…constrained geometry fit Major new module Modules/QSpaceFactorKernel.py implementing the equal-factor / symmetric-power factorization of the image-assignment kernel from report/interpolation/new_plan.tex. Key additions: - Two fit modes: 'individual' (uniform-class-sum manifold, unconstrained OMP) and 'constrained' (full C_n U Z = 0 pipeline with streaming target oracle) - Matrix-free Gram evaluation via translation-orbit correlations - Greedy rank selection with coverage/rel-residual diagnostics - In-process cache by structure fingerprint - QSpaceInterpolation runtime wiring: window_design='factor', d4_center='factor', factor_rank/far/cost/fit_mode options - Unit tests (test_factor_kernel.py, test_factor_kernel_snte.py, test_d4_center.py) - Factor locality audit script and results - Updated Implementation_plan.md and Interpolation_plan.md with detailed kernel-space yardsticks, locality audit, and next steps WARNING: This is a work-in-progress implementation. The constrained mode quality on SnTe D3 is currently plain-like (rel ~0.998); D4 reaches validation rel 0.914 but only 0.295 coverage. The individual mode beats the validated atomic kernel in kernel space but has not been tested on the full N=2000 physical benchmark. Needs per-class nonlocal leakage gates, cancellation convergence, q-smoothness checks, bootstrap variance analysis, and production spectrum benchmarks before any physics use. ABSOLUTELY NOT PRODUCTION READY.
(cherry picked from commit 620c527)
Propagate Raman invariant fix to atom-Fourier interpolation
The IR dielectric function was built as eps_inf + (4*pi/Omega)*chi_ionic with Omega = unit-cell volume, but the Lanczos perturbation carries a sqrt(n_cell) factor (prepare_ir scales the unit-cell charge vector by sqrt(n_cell)). The Green function therefore already includes n_cell, so Omega must be the supercell volume V = n_cell * V_unit_cell. Together with the factor of two in ir_susceptibility (Rydberg e^2 = 2), the total prefactor is 8*pi/V_supercell, matching the CellConstructor non-analytic LO-TO term. With the old unit-cell volume the LO mode was off by n_cell (verified on C2/c-24 H2 at 200 GPa: LO_x = 4111.1 cm-1 only reproduces the library value with V_supercell). - Store supercell_volume_angstrom3 in the execution manifest (backward-compatible fallback to n_cell * V_unit for old manifests) - Document the sqrt(n_cell) / Rydberg-8pi convention in the ir_susceptibility and dielectric_function docstrings - Update docs/spectroscopy.md - Add a regression test for the supercell-volume prefactor
Spectroscopy handed a fully loaded sscha.Ensemble.Ensemble to the backend
constructor on every MPI rank. With a 160k-configuration ensemble that is a
hard memory blocker: the configurations, and the X_q/Y_q arrays built from
them, are replicated n_procs times.
The driver now takes the *location* of the ensemble (EnsembleSource, or
Spectroscopy.from_ensemble_path) and routes the q-space backends through the
existing distributed loaders, so the master reads the configurations once and
scatters them and every rank keeps only N/n_procs of them.
Why the interpolated backend could not do this before, and what changed
----------------------------------------------------------------------
Building QSpaceAtomFourierLanczos interpolates the dynamical matrix, which
goes through CellConstructor's ForceTensor.Tensor2: Center() ends with five
Settings.broadcast and Apply_ASR() with one. Those assume every rank calls
them together. In a master-only build only the master enters them.
A stack dump shows this is worse than a deadlock. The master blocks in the
ASR broadcast, but the worker is already *past* its metadata bcast: MPI
matched it against Center()'s first broadcast and handed it the force-constant
tensor as if it were the loader's metadata. Nothing raises. A slightly
different arrangement of collectives does not hang at all and simply returns a
plausible, wrong spectrum.
The fix moves the collective work instead of suppressing it. The
interpolation depends only on the dynamical matrix, never on the
configurations, so every rank runs it before the master/worker split; only the
ensemble read stays master-only. Collectives match because every rank really
does execute them. CellConstructor is unchanged, and the private
monkeypatching workaround (tdscha_master_load.py in the hydrogen scripts) is
no longer needed by anything.
The alternative -- a CellConstructor suppress_collectives() context manager --
was rejected: a global mutable switch whose safety contract can only be
checked by reading every function reachable from inside the block, including
future ones, and it would tie tdscha to a new CellConstructor release.
Changes
-------
QSpaceInterpolation
FineHarmonicInterpolation and build_fine_harmonic(): the harmonic
interpolation as one immutable value, carrying a digest of every input it
read. validate_for() refuses an interpolation from a different matrix,
mesh, or ASR/effective-charge/LO-TO setting, so it can never contract the
ensemble in a mode basis that is not its own. The digest casts to one
canonical numeric type on purpose: passing a dynamical matrix through an
Ensemble demotes its Gamma block from complex128 to float64 without
changing a value, and a byte-level digest would reject every distributed
interpolated run.
QSpaceLanczos
prepare_distributed_construction(): a classmethod returning extra
constructor arguments, run by all ranks before the master/worker split.
The base returns {}; QSpaceAtomFourierLanczos overrides it.
load_distributed_atom_fourier_tdscha no longer passes
build_on_all_ranks=True, which survives only as a diagnostic oracle.
The metadata broadcast now carries a sentinel the workers check: a
mismatched collective aborts with the diagnosis instead of hanging or
returning a wrong spectrum.
Spectroscopy
EnsembleSource(data_dir, population, dyn, T, nqirr=, n_configs=,
final_dyn=, final_nqirr=, final_T=) and from_ensemble_path(). A loaded
ensemble is still accepted -- the replicated path, right for small systems
and required by backend="real", whose real-space Lanczos parallelizes over
a replicated ensemble by design; a q-space backend given one under mpirun
now warns. reference_dyn/reference_temperature replace the direct
self.ensemble.current_dyn reads.
The engine is built once per run() rather than once per perturbation:
preparing a perturbation resets the whole Lanczos state, and reading a
production ensemble is minutes of I/O. It is built lazily, so a fully
restored calculation reloads nothing.
Manifest schema 3 records the ensemble identity (population, n_configs,
temperatures, reweighted, directory name) next to the fingerprint of the
reference dynamical matrix; previously two different populations from the
same matrix fingerprinted identically. Absolute paths go in
ensemble_provenance, which is not validated, so moving a finished
calculation does not invalidate its checkpoint. ensemble_fingerprint()
became reference_fingerprint(dyn, T) and no longer touches the
configurations, which after a distributed load exist only on the master.
Tests
-----
New: the master-only interpolated path against the build-everywhere oracle
(this is what would catch a mode-basis gauge disagreement); a subclass that
deliberately leaves the collective in its constructor must abort, not hang and
not succeed; Spectroscopy from a source on one rank vs two, both q-space
backends, with every rank's N < N_global and the shares summing to N_global;
one engine for two runs and none once restored; restarted analysis parity;
FineHarmonicInterpolation validation including the Ensemble round trip;
EnsembleSource validation and backend routing.
One rank vs two agree to 4e-17 (q-space) and 1.3e-15 (atom-Fourier) relative.
Full suite: 172 passed, 8 skipped.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013ccBq6B2HFRZzmt71qc71b
Collaborator
Author
Distributed ensemble loading (pushed in efe9424)
Why the interpolated backend could not do this beforeBuilding A stack dump shows this is worse than a deadlock: the master blocks in the ASR broadcast, but the worker is already past its metadata The fix moves the collective work rather than suppressing it: the interpolation depends only on the dynamical matrix, so every rank runs it before the master/worker split, and only the ensemble read stays master-only. CellConstructor is unchanged. The alternative considered — a A sentinel in the metadata broadcast now makes the silent-corruption mode impossible: if a collective is ever mismatched again the job aborts with the diagnosis. There is a test that reintroduces the defect on purpose and asserts it aborts. CI diagnosis (run 31111882964, before this push)3 failed, 145 passed, 17 skipped. All three failures have the same root cause and none is caused by this PR's code: Every one is a test that shells out to In-process Julia works — all the other Julia tests pass in the same job. Julia only fails to initialise in Python processes launched by The likely cause is the CI's Julia backend choice. This is pre-existing and repo-wide. CI has been red on Two candidate fixes, neither applied here (they change CI/test semantics and are separable from this PR):
Expect this push to add failures, not remove them. It adds four more Not validated at production scale: everything above is the 2×2×2 test ensemble. The 250 GPa 8×6×8 hydrogen run is the real check of the memory envelope. 🤖 Generated with Claude Code |
The testsuite job failed with three MPI tests dying the same way in the child process: julia.core.JuliaError: Exception 'InitError' occurred while calling julia code: const PyCall = Base.require(Base.PkgId(... "PyCall")) In-process Julia worked -- every other Julia test passed in the same job. Julia only failed to initialise in interpreters launched through mpirun. The workflow installed PyJulia (pip install julia, then julia.install()) and never installed juliacall, so both python-sscha and tdscha fell back to the PyCall backend. That backend needs PyCall built against the running interpreter and shares one ~/.julia/compiled cache: test_mpi_parallel runs first with -np 2, and every later mpirun child fails afterwards, including the single-rank ones. The runner's Julia moved to 1.12 in the meantime, which PyCall/PyJulia handle poorly; this used to pass in June. Switch to juliacall, which is what pyproject declares as the julia extra and what JuliaExt prefers when both are present. It has no libpython coupling, so it works in any interpreter, and it does not pay PyJulia's import cost. SSCHA_JULIA_BACKEND=juliacall is set at the workflow level so a transitively installed pyjulia can never silently take over. A warm-up step boots Julia once, serially, before pytest. The mpirun tests launch several interpreters at the same time; if those were the first to touch ~/.julia they would race on the same depot, which is the same failure mode again with a different backend. Also: - Drop aiida-core, aiida-quantumespresso, pgtest and the rabbitmq service. Nothing in tests/, Modules/ or scripts/ references any of them; they were template leftovers costing several minutes of install and a pseudopotential download. matplotlib is now installed explicitly because it was arriving transitively through aiida and is imported at module level by collected tests. - Install tdscha.testing (meson.build). scripts/validate_docs.py imports tdscha.testing.test_data from the *installed* package, and it was never installed, so the "Validate documentation" step fails as soon as pytest stops failing before it. With this it passes. - checkout@v2 -> v4 and setup-python@v2 -> v5; both were being force-run on Node 24 with a deprecation annotation. - requirements.txt: julia -> juliacall, matching pyproject. Verified locally with SSCHA_JULIA_BACKEND=juliacall, OMP_NUM_THREADS=1 and JULIA_NUM_THREADS=1: 172 passed, 8 skipped, and validate_docs.py --test-imports --doctest exits 0. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013ccBq6B2HFRZzmt71qc71b
The juliacall switch fixed the Julia failures (145 -> 165 passing, and test_mpi_parallel now runs), but the distributed tests still failed: assert 10 < 10 # each rank should hold a strict subset On a two-rank run every rank owned the whole ensemble. The log explains why: the string "rank 1" does not appear anywhere in it. "mpirun -np 2" was starting two INDEPENDENT single-rank jobs, each with its own singleton COMM_WORLD. mpi4py came from a pip wheel that bundles its own MPI runtime, installed next to the apt mpich that owns mpirun. The two do not interoperate, so MPI_Init in each child produced a communicator of size 1. Build mpi4py from source instead, and assert in the install step that a two-rank job really reports size 2. This was not new. It means every MPI test in this repository has been passing vacuously for as long as the wheel has been used: test_julia/test_mpi_parallel compares a serial run against an "MPI" run that was also serial, so it agreed trivially. The distributed tests added in the previous commit are the first ones that assert the configurations were actually partitioned, which is why they are the first to notice. A silently vacuous MPI suite is worse than a failing one, so tests/conftest.py adds a multi_rank_mpirun fixture: it probes once per session, skips when there is no MPI at all, and fails with the diagnosis and the fix when there is an mpirun that does not produce a multi-rank job. The distributed tests take the launcher from that fixture. This changes no pass/fail outcome -- it turns "assert 10 < 10" into an explanation. Verified locally: the eight distributed tests pass against a real two-rank mpirun, and against a stub launcher that reproduces the CI condition the fixture reports exactly what is wrong and how to fix it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013ccBq6B2HFRZzmt71qc71b
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds a restartable, backend-neutral
Spectroscopydriver for one-phonon Raman and IR responses (real, q-space, and atom-Fourier backends), with symmetry reduction that deduplicates symmetry-equivalent perturbations into a single Lanczos run and reconstructs the full observable.Key components:
Modules/Spectroscopy.py— public request/perturbation API, symmetry orbits and cosets, manifest-driven checkpointing, analysis helpers.Modules/_SpectroscopyWorkflow.py— execution, persistence, and backend orchestration.eps_inf + (4*pi/Omega) * chi_ionicwithOmega= supercell volume, so the total prefactor is8*pi/V_supercell(Rydberge^2 = 2, perturbation carriessqrt(n_cell)), matching the CellConstructor non-analytic LO-TO term.Verification
tests/test_spectroscopy/passes (33 tests). The IR prefactor is cross-checked on C2/c-24 H2 at 200 GPa: LO_x = 4111.1 cm-1 reproduces the CellConstructor library value only with the supercell volume.