HPhi solver: off-diagonal Gf fix, speedups, canonical & spin-orbit paths#174
Open
k-yoshimi wants to merge 18 commits into
Open
HPhi solver: off-diagonal Gf fix, speedups, canonical & spin-orbit paths#174k-yoshimi wants to merge 18 commits into
k-yoshimi wants to merge 18 commits into
Conversation
…imize einsum HPhi (hphi_spectrum.py): - Fix the multi-orbital off-diagonal G reconstruction. The previous code assigned G_ij twice to the same element (dead store), never set the transposed G_ji, subtracted the diagonal twice inside the idx loop, and could overwrite the diagonal block. Now G_ij and G_ji are both filled from the same excitation data with the diagonal subtracted once from each combination (B=c_i+c_j, A=c_i+i c_j). - Speed up: compute only the upper-triangular (a_i <= a_j) excitation tasks and reconstruct G_ji from the same data; skip the unused diagonal idx=1 tasks. This roughly halves the number of HPhi subprocess invocations, which scale as O(n_orb^2). - Guard against an empty task set (n_site == 0) with a clear error. - Remove the dead, now-divergent CalcSpectrum class and its dangling commented-out references. Interaction basis rotation (base.py, interaction.py): - Add optimize=True to the four-rotation U/X/chi einsum contractions, reducing them from O(N^8) to O(N^5) for large spin-orbital dimension N. Tests (tests/non-mpi/hphi/test_reconstruction.py): - New HPhi-binary-free unit tests for the off-diagonal/diagonal reconstruction and the reduced parallel driver (task pruning + assembly), plus the empty-task guard. Note: the off-diagonal reconstruction convention is verified algebraically and by these unit tests; an end-to-end multi-orbital HPhi run is still pending external validation data. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
HPhi is a registered solver but had no documentation. Add a reference page (features, install, the exec_path/n_bath/exct/fit_gtol parameters, and the power-of-four MPI requirement) and link it from the impurity-solver toctree. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Two ED solvers must give the same self-energy for the same impurity model. This cross-check catches bugs that single-solver unit tests miss -- notably an insufficient exct: the 2-orbital Kanamori atom has a 4-fold degenerate 1-electron ground state, so with the default exct=1 HPhi samples only one member of the multiplet and spuriously breaks orbital symmetry (asymmetric diagonal, large fake off-diagonal). With exct spanning the full space the finite-T thermal trace is complete and HPhi reproduces the scipy/sparse reference to ~1e-6. The test runs only when an HPhi executable is available (DCORE_HPHI_EXEC env var or HPhi on PATH); otherwise it skips, so it is safe in CI. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
A degenerate ground state (common in multi-orbital models) needs exct large enough to cover the whole low-energy multiplet at finite T; the default exct=1 otherwise samples a single member and silently breaks orbital symmetry. After the eigenvalue step the solver now inspects the eigenenergies and warns when the highest computed state still carries Boltzmann weight above exct_weight_threshold (default 1e-3). Documented in the HPhi manual page with a dedicated section, and covered by pure-Python unit tests. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The off-diagonal Green's function is reconstructed from two excitations, B = c_i + c_j and A = c_i + i c_j. The creation channel (c^dag) of A must use the conjugate operator A^dag = c_i^dag - i c_j^dag, but the imaginary coefficient was written as +i for both the annihilation and creation channels. This left the anti-symmetric (imaginary) part of the off-diagonal Gf with a wrong high-frequency tail: the off-diagonal self-energy stayed correct at low frequency but diverged linearly with the Matsubara frequency, and the 2x2 inversion spread the error onto the diagonal too. A purely real off-diagonal (B channel) was unaffected. Verified against the scipy/sparse ED solver: with the fix HPhi reproduces the full omega_n dependence of a complex off-diagonal crystal-field model to ~1e-5 (was a divergence of ~90). Added an off-diagonal cross-check test that compares over the whole Matsubara axis, which is what exposes this class of tail bug. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
The one-body Green's function is computed from many independent HPhi runs. They were always run serially inside a node-local ProcessPool of size np. Add an n_procs_per_hphi (n_inner) option so each HPhi run can use n_inner MPI ranks while n_outer = np // n_inner runs execute concurrently, splitting the total np between the two levels. Each task already runs in its own work directory, so concurrent mpirun launches do not collide. n_procs_per_hphi defaults to 1, which reproduces the previous behaviour exactly (serial HPhi, np concurrent runs) and is exercised end-to-end by the cross-check tests. n_inner is clamped to a power of four (HPhi requirement). The split arithmetic and launcher-string construction are covered by unit tests; running n_inner > 1 end-to-end needs an MPI build of HPhi. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- Reword the exct warning to reflect what is actually provable: w_last is an UPPER BOUND on every omitted state's Boltzmann weight (all omitted states lie at energy >= e_last), so exceeding the threshold means truncation cannot be ruled out, not that states are definitely missing. Avoids over-claiming on the near-degenerate-then-gap false-positive case. - Add unit tests that guard the c_i + i c_j conjugation fix by inspecting the generated single_ex.def for both the annihilation (ex_state=0, +i) and creation (ex_state=1, -i) channels, and that the real B-combination is unaffected. - Extend the Gf parallel-layout tests with a non-divisible split and with mpi_prefix propagation through a 9-element p_common into the HPhi command (serial vs mpirun), plus a test documenting the conservative warning bound. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
- Add an opt-in integration test that verifies the two-level Gf parallelism is split-invariant: running each HPhi of the Gf step under mpirun -np 4 (n_procs_per_hphi=4) reproduces the serial layout's self-energy on a complex off-diagonal model. Gated on DCORE_HPHI_MPI_EXEC + DCORE_HPHI_EXEC + mpirun, so it is skipped in CI; ~35s when enabled. - Document running HPhi in parallel: the two MPI phases (eigenvalue step uses all np; Gf step splits into n_outer x n_procs_per_hphi), the MPI-only-binary launcher caveat, and a SLURM/ohtaka (ISSP System B) job-script example. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
HPhi writes the eigenvectors MPI-distributed (zvo_eigenvec_{i}_rank_{r}.dat, one
file per rank) in the eigenvalue step, and the Green's-function step reads them
back. An MPI-distributed eigenvector can only be read by the same number of
ranks that wrote it; otherwise HPhi stops at "An Eigenvector is inputted in
CalcSpectrum" and no DynamicalGreen is produced.
Previously the eigenvalue step used all np (rounded to a power of four) while the
Gf step used n_inner = n_procs_per_hphi, so any run with n_inner != np_power4
failed -- including the default n_procs_per_hphi=1 with np>1. gf_parallel_layout
now returns a single launcher used for BOTH phases (n_inner ranks each), with
n_outer = np // n_inner concurrent Gf runs. np itself need not be a power of
four; only n_inner must be.
This also corrects the earlier "concurrent bare runs collide" diagnosis: the
failure was the rank-count mismatch, not concurrency -- n_inner=1 with n_outer=4
now runs and matches the serial reference to ~3e-5.
For an unparseable launcher (rank count not the last token, e.g.
"mpirun -np 16 --bind-to core") both phases use the same given launcher and the
Gf step is serialized (n_outer=1) to keep the rank invariant. Docs updated.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
From the System B user guide: Slurm with srun, and concurrent srun steps inside one allocation need the "bulk job" pattern (srun --exclusive --mem-per-cpu). Update the HPhi parallel section to use [mpi] command = srun --exclusive --mem-per-cpu=1840 -n # and note the i8cpu interactive/debug queue (<=8 nodes, 30 min) and salloc-based interactive testing. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
…tion Add an opt-in fast path to the HPhi spectrum driver that uses HPhi's new SpectrumLoopExct feature: instead of relaunching HPhi once per (eigenstate x excitation operator), HPhi loops over the thermally relevant eigenstates internally, so the Gf step launches HPhi once per operator. This removes the per-eigenstate process-startup overhead, which dominates the cost for small/medium Hilbert spaces. Enabled by the environment variable DCORE_HPHI_INTERNAL_LOOP=1; default off keeps the validated per-(idx) relaunch path as the reference fallback, so this is a no-op unless explicitly requested (and an HPhi build providing SpectrumLoopExct is used). In loop mode CalcSpectrumCore: - writes a single namelist with SpectrumVec as the common eigenvector base name (HPhi appends _<idx>_rank_<r>.dat per eigenstate) and SpectrumLoopExct set to exct_cut in modpara; - stages <header>_energy.dat into the per-element calc dir, since HPhi reads each eigen-energy E_idx from there (childfopenMPI prepends "output/"); - launches HPhi once and reads back the per-eigenstate <header>_DynamicalGreen_<idx>.dat files it writes directly. Validated against the independent scipy/sparse ED solver: loop mode reproduces the diagonal and complex off-diagonal Green's functions across the full Matsubara axis, is invariant to the MPI rank split, and is ~3.3x faster (95s -> 29s on the cross-check) than the per-launch path. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
HPhi-level golden test for the SpectrumNumOp op-inner eigenvector reuse: drives HPhi directly (standard mode for the eigenvalue step, expert mode for the spectrum) on a 2-site Hubbard dimer and checks that evaluating two operator sets in one run (SpectrumNumOp=2) is byte-identical to running each operator separately (SpectrumNumOp=1), for genuinely distinct same-sector operators. Two negative tests confirm a clean abort with no partial output when set 1 is in a different Hilbert sector or its def file is truncated. Gated on DCORE_HPHI_EXEC (skipped in CI without an HPhi executable). Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
…ctor reuse)
Builds on the internal eigenstate loop: in DCORE_HPHI_INTERNAL_LOOP mode the
one-body Green's function operators are now grouped by Hilbert sector and each
group is computed in a single HPhi launch using SpectrumNumOp, so every
eigenvector is read once and reused across all operators of that sector
(op-inner). The eigenvector reads drop from exct*(#operators) to
exct*(#sectors), which is the dominant cost for large Hilbert spaces where the
eigenvector I/O — not the BiCG solve — limits the per-operator path.
Grouping: a single-excitation operator changes the electron number by +-1 in
spin sigma, so diagonal c_{i,sigma} and same-spin off-diagonal
c_{i,sigma}+i c_{j,sigma} share the sector keyed by (ex_state, sigma) and are
batched. Cross-spin off-diagonal operators (sigma_i != sigma_j) mix Sz and are
sector-inconsistent, so they run as singletons (one launch each). HPhi verifies
the shared sector at run time (and, for grand-canonical models, simply reuses the
common space).
Implementation (internal-loop path only; default path unchanged):
- _group_tasks_by_sector / calc_one_body_green_batch: group tasks, run one
HPhi launch per batch, route results back to task order via a hashable
task key (task[:6], dropping the unhashable p_common).
- CalcSpectrumCore.get_one_body_green_batch: write operator set 0 to
single_ex.def (the namelist SingleExcitation) and sets 1.. to
single_ex_<op>.def, set SpectrumNumOp, run once, and Boltzmann-sum the
per-(idx, op) DynamicalGreen output per operator.
- Error handling: _run_HPhi raises on a non-zero HPhi exit; batch directories
are cleaned in a finally; assert exact-once task coverage and op/result
alignment so any routing regression fails immediately and locally.
- New pure-Python routing test (no HPhi needed) covering exact-once coverage,
sector consistency, cross-spin singletons, and task-order reassembly.
Validated end-to-end against the scipy/sparse ED solver (grand-canonical model):
the batched path reproduces the diagonal and complex off-diagonal self-energy
across the full Matsubara axis, is MPI split-invariant, and the default
(non-loop) path is unchanged.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
Add an opt-in canonical path to the HPhi solver that replaces the single grand-canonical run with per-(Ne, 2Sz)-sector canonical runs. The Hamiltonian is block-diagonal in (Ne, 2Sz) when 2Sz is conserved, so the grand-canonical finite-temperature trace equals the weighted sum of the per-sector contributions. Each sector spans a much smaller Hilbert space than the full 4**n_site grand-canonical space, which is where the dominant Green's-function (BiCG) cost shrinks for large systems -- the same idea that makes the scipy/sparse ED solver efficient. Enabled by DCORE_HPHI_CANONICAL_SECTORS=1; default (unset) keeps the unchanged grand-canonical path. Guarded on `not use_spin_orbit`: spin-orbit one-body terms break 2Sz (Ne is always conserved by the hopping + density interactions, and the supported interactions -- Kanamori/Slater/Respack -- are spin-diagonal). - enumerate_particle_sectors(n_site): the (Ne, 2Sz) sectors that partition the grand-canonical space (with unit tests). - solve(): per sector, write Ncond/2Sz + CalcModel=0 (canonical Hubbard), run the eigenvalue step, then the Gf step, and recombine exactly with G = sum_s w_s G_s / sum_s w_s, w_s = Z_s exp(-beta (Emin_s - Emin_global)), reproducing the single global Boltzmann sum. The per-sector thermal-truncation warning is preserved. - calc_one_body_green_core_parallel gains an optional sector_occupancy filter: operators with an exactly-zero contribution in the sector (cross-spin off-diagonal -- zero by spin conservation; annihilation/creation on an empty/full spin) are skipped, which also avoids HPhi excitations into non-existent sectors. Unit tests cover the filter. - set_energies() uses min/max so the sector-local Boltzmann normalization matches the cross-sector recombination regardless of energy ordering. Requires the companion HPhi fix that re-orthonormalizes degenerate eigenvectors: the small canonical sectors are exactly degenerate, and non-orthonormal degenerate eigenvectors would otherwise break the per-sector trace (spurious orbital-symmetry breaking). Validated end-to-end against the independent scipy/sparse ED solver: the canonical path reproduces the diagonal and complex off-diagonal self-energy across the full Matsubara axis; the grand-canonical path is unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
Phase 1 of the canonical path computed every (Ne, 2Sz) sector (exact but with no speedup -- many sectors carry negligible Boltzmann weight). Phase 2 adds a cheap ground-state-only pre-pass per sector, then runs the full eigenvalue + dominant Green's-function step only for the thermally relevant sectors. Only the thermally occupied sectors need a run: the one-body spectrum reaches the (Ne+-1) excited sectors as intermediate states resolved internally by HPhi's BiCG, not as separately summed eigenstates. A sector is kept when its (bounded) thermal contribution exceeds `canonical_sector_weight_threshold` (default 1e-8; 0 keeps every sector = the exact Phase-1 behavior). The bound is `min(exct, dim) * exp(-beta (Egs_sec - Egs_global))`: since the sector partition factor Z_sec is at most the number of states summed, this never drops a sector whose many low-lying / nearly degenerate states make its full contribution non-negligible (a ground-state -weight-only filter could). The threshold is validated to lie in [0, 1). The per-sector eigenvalue runs reuse one work dir sequentially (output/ read immediately after each run; only the Gf step within a sector is parallel). Validated vs the scipy/sparse ED solver: the cross-check keeps 6 of 9 sectors and still reproduces the diagonal and complex off-diagonal self-energy; for larger systems at fixed filling the kept fraction shrinks, which is the speedup. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
Golden test for the HPhi SpectrumNumBra feature (Stage-3 bra/ket reuse): one ket BiCG solve projected onto L bras must equal L separate single-bra runs to <1e-10, and a bra set in a different Hilbert sector than the ket must be rejected with no partial output. Drives HPhi directly on a 2-site Hubbard dimer; gated on DCORE_HPHI_EXEC so it skips in CI. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
…b^2 -> n_orb
Stage-3 opt-in path that replaces the c_i + i c_j combination trick for the
off-diagonal one-body Green's function with a direct bra/ket reconstruction:
each ket c_{j,sigma} is solved by BiCG ONCE and projected onto every bra
c_{i,sigma} via HPhi's SpectrumNumBra, so the whole same-spin G block comes
from n_site solves instead of ~n_site^2 combination solves.
- hphi_spectrum.py: env DCORE_HPHI_BRAKET=1 routes to _braket_one_body_green.
_braket_sector_jobs enumerates one launch per (ex_state, spin) sector
(kets = bras = all same-spin sites); get_one_body_green_braket_batch writes
the ket/bra single_ex(.bra) files, sets SpectrumNumOp/SpectrumNumBra, and
reads DynamicalGreen_<idx>_<op>_<bra>.dat. The creation channel (ex_state 1)
builds bra/ket from c^dag, which reverses the two indices, so its
off-diagonal is emitted transposed; the result is stored at
one_body_green[ket][bra] to match the validated combination-trick
convention. _finite_T_spectrum_braket handles the unsuffixed single-site
(num_op == num_bra == 1) file name. Cross-spin elements are not computed
(zero by spin conservation; the final Gimp keeps only same-spin blocks).
- hphi.py: DCORE_HPHI_BRAKET=1 implies the canonical (Ne, 2Sz) sector path it
is built on, and is rejected with spin-orbit (2Sz not conserved). The
grand-canonical bra/ket route and n_site < 1 are rejected explicitly.
- Default (combination-trick) path is unchanged. Validated end-to-end vs the
scipy/sparse ED solver (diagonal symmetric + complex off-diagonal, to the
test tolerance over the whole Matsubara axis); element-wise it matches the
combination-trick same-spin blocks to ~1e-10.
Requires an HPhi build with SpectrumNumBra and the HubbardGC single-excitation
off-diagonal shift. New CI-safe routing tests in test_braket_routing.py.
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
…ed sectors)
Generalize the Stage-3 direct bra/ket Green's function to spin-orbit-coupled
impurities (2Sz NOT conserved). When 2Sz is broken but Ne is still conserved,
sector by Ne ONLY (HubbardNConserved: write Ncond, omit 2Sz) instead of by
(Ne, 2Sz): both spins then share each Ne+-1 excited space, so one ket BiCG solve
yields the same-spin AND cross-spin G_{i,sigma_i; j,sigma_j} from one launch.
- hphi_spectrum.py: enumerate_ne_sectors; _braket_ne_sector_jobs (both-spin
kets=bras per Ne sector); _braket_store_op (pure helper centralizing the
creation-channel transpose, same convention for same- and cross-spin);
_braket_assemble (pure tensor accumulation, cross-spin tested);
_braket_one_body_green / calc_one_body_green_core_parallel gain ne_only.
- hphi.py: use_ne_canonical branch (spin_orbit + braket -> Ncond-only modpara ->
HubbardNConserved, ne_only Gf, exact Ne-sector recombination). The boundary Ne
sectors (Ne=0 vacuum; Ne=2*n_site full) can't run via HubbardNConserved and are
dropped with a loud warning (exact when those near-empty/near-full states are
thermally negligible). The spin-orbit self-energy now computes: the generic
Dyson Sigma=G0^-1-Gimp^-1 block code already handles the single combined 'ud'
block, so the previous NotImplementedError is removed.
DCORE_HPHI_FORCE_NE_SECTORS=1 forces the Ne-only route on a 2Sz-conserving model
(testing / alternative decomposition).
Validated end to end against the scipy/sparse ED solver: a spin-orbit model with
a spin-flip crystal field reproduces the full 2*norb x 2*norb self-energy --
including the genuinely non-zero cross-spin blocks -- to ~1e-4
(test_hphi_braket_spin_orbit_matches_scipy_sparse). The force-NE route reproduces
the (Ne,2Sz) result on a spin-diagonal model
(test_hphi_braket_ne_only_matches_scipy_sparse). CI-safe routing/enumeration unit
tests in test_braket_routing.py + test_ne_sectors.py. Default combination-trick
and the (Ne,2Sz) braket path are unchanged.
Requires an HPhi build with the HubbardNConserved single-excitation off-diagonal
spectrum + sz() boundary fix (companion HPhi change).
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_018pa872xYjHXtjFRyUn5ehj
k-yoshimi
changed the title
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
Improves the HPhi impurity solver's finite-temperature one-body Green's function: a correctness fix for the off-diagonal components, substantial speedups (sector batching, two-level parallelism, internal eigenstate loop, direct bra/ket projection), a particle-number-resolved canonical path, and a spin-orbit-capable route — all cross-validated against the independent
scipy/sparseED solver.Highlights
Correctness
c_i + i·c_jexcitation was mis-conjugated) — the diagonal-only result was already correct; off-diagonal blocks now matchscipy/sparseED.exctis too small to capture the thermal trace at the requestedT.Performance (all opt-in via
DCORE_HPHI_*env vars; default path unchanged)SpectrumNumOpinternal loop).n_inner × n_outer) MPI parallelism for the Gf step, with the eigenvalue and Gf steps pinned to the same rank count (split-invariance tested).n_orbbras, cutting the shifted-BiCG count fromn_orb²ton_orb(SpectrumNumBra).New capabilities
(Ne, 2Sz)sector selection.2·norb × 2·norbself-energy including non-zero cross-spin blocks (validated to ~1e-4 vsscipy/sparsewith a spin-flip crystal field).Tests & docs
scipy/sparseED cross-validation test; golden checks for theSpectrumNumOp/SpectrumNumBrafinite-T loops; sector-enumeration and reconstruction unit tests (CI-safe, no external solver needed for the unit ones).Dependency
The bra/ket and spin-orbit paths exercise HPhi spectrum features upstreamed separately — HPhi PR #254 (
SpectrumNumBra) and the stacked HubbardNConserved off-diagonal PR. They are gated behindDCORE_HPHI_*env vars, so this PR is safe to merge ahead of an HPhi release; those code paths simply require a sufficiently new HPhi build at runtime. The always-on off-diagonal correctness fix has no such dependency.🤖 Generated with Claude Code