Skip to content

Give SCIState its amplitudes in PySCF's sign convention (backport #62) - #63

Merged
Jim Garrison (garrison) merged 7 commits into
stable/1.7from
mergify/bp/stable/1.7/pr-62
Oct 16, 2026
Merged

Jim Garrison (garrison) merged 7 commits into
stable/1.7from
mergify/bp/stable/1.7/pr-62

Conversation

@mergify

@mergify mergify Bot commented Oct 16, 2026

Copy link
Copy Markdown

Summary

SCIResult.sci_state from sbd.sbd_solver.solve_sci / solve_sci_batch held SBD's amplitudes in SBD's own sign convention, but qiskit-addon-sqd's SCIState hands them to PySCF, which assumes a different one. Everything SCIState computes from the amplitudes was therefore wrong, without any error:

full space sci_state.spin_square() energy rebuilt from sci_state.rdm()
H3 doublet 0.7599 (should be 0.75) 85 mHa off
H4 singlet 0.0356 (should be 0) 89 mHa off

Not affected: energies, SCIResult.rdm1/rdm2 (SBD's own RDMs), orbital occupancies, carryover selection (amplitude magnitudes only), and everything inside qiskit-addon-sqd's loop. None of our examples or the qiskit-addon-sqd guides call the affected methods. It only shows when a user calls sci_state.rdm() or sci_state.spin_square(), or uses the amplitudes with PySCF.

Cause

The two conventions order a determinant's creation operators differently. SBD goes orbital by orbital (alpha 0, beta 0, alpha 1, ...); PySCF puts every alpha operator before every beta one. The magnitudes agree, but each determinant differs by the reordering phase (-1)^n, where n counts, over the occupied beta orbitals p, the occupied alpha orbitals above p.

Fix

  • _solve_sci_core applies that phase when it builds SCIState, in row blocks so the temporary stays bounded at large subspaces.
  • The tpb_diag docstring now states that the raw wavefunction files SBD writes keep its own convention, and gives the conversion.

Tests

test/test_sci_state.py solves each subspace with SBD through solve_sci and with qiskit-addon-sqd's PySCF solver, then compares the states themselves, not only the energies:

  • elementwise, up to one global sign;
  • spin_square();
  • the energy from rdm().

Cases are full spaces (an H3 doublet, H4 and H6 singlets) and a seeded random truncation of H6. Each runs single-process and under MPI. Full spaces are included because a small truncation can happen to need no sign change, so a lone small case can pass by accident.

Without the fix all four cases fail (state overlap 0.95–0.97 while energies agree). With it, they pass:

  • locally under tox -e py,notebook,mpi,docs;
  • on an H100 node for cpu, gpu (Thrust) and gpu-omp, single-process and on 2 ranks, so all three backends share the convention.

Compared test by test with main on that node, no existing test changed outcome.

🤖 Generated with Claude Code


This is an automatic backport of pull request #62 done by [Mergify](https://mergify.com).

…der (#55)

* Document how to retrieve TPB amplitudes, and the bitstring packing order

Two things a user has to discover by reading upstream C++ or experimenting.

TPB has two wavefunction output paths that write different formats, and
nothing said so. `savename` -- the one named in `tpb_diag`'s signature --
goes to `SaveWavefunction`, which writes one file per b_comm position
holding only that rank's block, behind a three-`size_t` header and the
block's determinant words. `sbd_data.dump_matrix_form_wf` goes to
`SaveMatrixFormWF`, which gathers the full adet x bdet subspace onto one
rank and writes a bare row-major float64 array with no header at all, or,
for a `.dat`/`.txt` path, a text table labelling each amplitude with its
own alpha and beta bitstring.

The second is what almost every caller wants, and is what this package's
own SQD integration already uses (`sbd_solver.py` sets it and reads the
result with a bare `np.fromfile`). It appeared only as a one-line pybind
docstring and a CLI flag in `run_sbd_diag.py`, absent from the README's
TPB section and from the `tpb_diag` docstring, both of which mention
`savename` alone -- so following the signature led to parsing the wrong
layout. Document both, side by side, and note that the text form is
written with default stream formatting and so carries ~6 significant
digits rather than full float64.

`from_string` and `makestring` packed words with no stated convention.
Bitstrings are packed from the right: the trailing `bit_length` characters
become word 0, so the leading characters of a multi-word string land in
the last word, not the first. Easy to get backwards, and a wrong guess
yields a valid-looking but wrong subspace rather than an error.

GDB needed no changes here: #46 documented its amplitude retrieval and
file layout in examples/gdb/README.md.

Assisted-by: Claude Opus 5

* Fit the amplitude and packing docs to the API surface #53 defined

#53 curated the autodoc list and added field tables, which changes where
this branch's documentation belongs.

The packing convention moves from `from_string` to `from_strings`. #53
documents the bulk form and deliberately leaves the singular one out, so
the convention was written on a function the rendered docs never show --
and `from_strings`, the form users are steered to, had no packing
documentation of its own. `from_string` and `makestring` now point at it
instead of restating it. The examples are shown as the `uint64` ndarray
`from_strings` actually returns rather than as nested lists.

`dump_matrix_form_wf` now appears in #53's `TPB_SBD` field table as well
as the README's, so the README row defers to the format section instead of
describing the field a second time.

Also corrects a claim this branch made about `sbd_solver.py`: since #49 the
wavefunction write is skipped when SBD selects the carryover itself
(`carryover_type` 1 against an addon that accepts it), because the next
subspace then comes back in `carryover_adet`/`carryover_bdet` and the
amplitudes never reach Python. Saying it "uses" the dump without that
qualification implied the write always happens.

`tox -e docs` passes with no warnings, and the packing convention renders
on the `from_strings` entry.

Assisted-by: Claude Opus 5
The Darwin block probes $CONDA_PREFIX/include/omp.h and falls back to Homebrew.
That misses the case where the activated environment is not the one being built
into: conda-build and rattler-build activate their own base installation and
install into $PREFIX, so a recipe that puts llvm-openmp in the host environment
still fails with

    Error: no OpenMP runtime found on this macOS host.
           Looked in $CONDA_PREFIX/include and /opt/homebrew/opt/libomp/include.

on a builder with no Homebrew libomp -- while the header sits unexamined in
$PREFIX. This is what the conda-forge osx_64 build hits.

Consult sys.prefix, then sys.base_prefix, after CONDA_PREFIX and before
Homebrew. sys.prefix is the prefix of the interpreter running setup.py, which
under both build backends is the environment the extension is built for; it is
also the convention _mpi_prefix_from_env_prefix() already uses, which is why MPI
detection succeeded in the same build where OpenMP failed.

CONDA_PREFIX stays first, so an ordinary activated environment behaves exactly
as before. The new candidates can only be reached where the code previously
fell through to Homebrew, so no build that passes today changes behavior.

This serves the intent of the original change rather than reversing it: build
against the libomp that will be LOADED at import time, resolved through the
interpreter's rpath. Linking the base installation's copy while installing into
$PREFIX is how a process ends up with two LLVM OpenMP runtimes, the abort
described in the README and tracked as #27.

The error message listed $CONDA_PREFIX literally, which no longer describes the
set of prefixes searched, so report the paths actually tried, deduplicated.

Assisted-by: Claude Opus 5
* docs: correct outdated and misleading API docstrings

- tpb_diag, tpb_diag_from_files, gdb_diag and get_backend list 'gpu-omp' and
  the aliases; gdb_diag says gpu-omp runs GDB on the host.
- get_backend says what None resolves to ('auto' before any init, 'cpu' once
  a diag call has initialized) and that cpu + gpu-omp must not share a process.
- init lists gdb_diag among the per-call device overrides.
- tpb_diag_from_files says its one file supplies both alpha and beta.
- TPB_SBD, GDB_SBD, LoadFCIDump and LoadAlphaDets get real docstrings.
- assemble_rdms cites upstream SBD's docs/user-guide.md, not a private repo.
- DeviceConfig's example passes the config to solve_sci_batch instead of the
  undocumented apply(); auto() drops development history.
- References to functions not on the API page become literals, so they no
  longer render as dead links; from_strings examples use bit_length 62.
- The bit_length row covers odd values in all GDB use, not only the heatbath
  expansion: upstream's own diag gives a wrong plain-diagonalization energy at
  31 on Fe4S4 (r-ccs-cms/sbd#104).

Docstrings and one table row only; no code changes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

* docs: LoadAlphaDets sorts and removes duplicates

Upstream's LoadAlphaDets ends with an unconditional sort_bitarray
(makedeterminants.h:190), so the result is sorted and distinct, not in file
order, and needs no further sort. Also note that only .txt and .bin are read:
any other extension returns an empty list without an error.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
* sbd_solver: give SCIState its amplitudes in PySCF's sign convention

SBD orders a determinant's creation operators orbital by orbital (alpha 0,
beta 0, alpha 1, ...); PySCF, and so qiskit-addon-sqd's SCIState, puts every
alpha operator before every beta one. The magnitudes agree but determinants
differ by (-1)^(sum over occupied beta p of the occupied alpha orbitals above
p), so what SCIState computes through PySCF from the amplitudes was wrong:
spin_square() 0.036 for an exact singlet (0.7599 for a doublet), and an
energy rebuilt from rdm() 85-89 mHa off. Energies, SBD's own RDMs,
occupancies and magnitude-based carryover were unaffected, which is why no
existing test saw it.

_solve_sci_core now applies that phase, in row blocks to bound memory.
test_sci_state.py compares the whole state with qiskit-addon-sqd's PySCF
solver on full and truncated subspaces (a doublet and two singlets): elementwise
up to a global sign, spin_square(), and the energy from rdm(). Without the
fix the overlaps are 0.95-0.97 and all four cases fail. The tpb_diag
docstring documents the convention of the files SBD writes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

* Apply batched suggestions from code review

Co-authored-by: Jim Garrison <jim@garrison.cc>

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Co-authored-by: Jim Garrison <garrison@ibm.com>
Co-authored-by: Jim Garrison <jim@garrison.cc>
* sbd_solver: give SCIState its amplitudes in PySCF's sign convention

SBD orders a determinant's creation operators orbital by orbital (alpha 0,
beta 0, alpha 1, ...); PySCF, and so qiskit-addon-sqd's SCIState, puts every
alpha operator before every beta one. The magnitudes agree but determinants
differ by (-1)^(sum over occupied beta p of the occupied alpha orbitals above
p), so what SCIState computes through PySCF from the amplitudes was wrong:
spin_square() 0.036 for an exact singlet (0.7599 for a doublet), and an
energy rebuilt from rdm() 85-89 mHa off. Energies, SBD's own RDMs,
occupancies and magnitude-based carryover were unaffected, which is why no
existing test saw it.

_solve_sci_core now applies that phase, in row blocks to bound memory.
test_sci_state.py compares the whole state with qiskit-addon-sqd's PySCF
solver on full and truncated subspaces (a doublet and two singlets): elementwise
up to a global sign, spin_square(), and the energy from rdm(). Without the
fix the overlaps are 0.95-0.97 and all four cases fail. The tpb_diag
docstring documents the convention of the files SBD writes.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

* Apply batched suggestions from code review

Co-authored-by: Jim Garrison <jim@garrison.cc>

---------

Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Co-authored-by: Jim Garrison <garrison@ibm.com>
Co-authored-by: Jim Garrison <jim@garrison.cc>
(cherry picked from commit 73c71a4)
Since qiskit 2.6.0 reached PyPI, every Linux CI job fails before running a
test: importing our extension raises "cannot allocate memory in static TLS
block" from libgomp, and the whole suite errors out. macOS is unaffected.

glibc sizes the per-thread TLS block at startup and keeps only a small
surplus for libraries dlopen'ed later that need initial-exec TLS. Every
Python extension module arrives that way, so they all compete for that one
surplus, and the library that fails is whichever asks last rather than
whichever took the space. qiskit 2.6.0's _accelerate grew from 736 to 1544
bytes of it, which leaves too little for the libgomp we pull in afterwards.
The growth is matrixmultiply 0.3.11 enabling its AVX-512 kernels, which
quadruple one thread-local scratch buffer; it takes effect at rustc >= 1.89,
and 2.6.0 also raised its MSRV to exactly that.

Raising the surplus restores the headroom. It goes in tox.ini rather than the
workflows so local runs get it too; passenv would only forward the variable
when the caller already set it, which is no help in CI. tox replaces setenv
rather than merging it in a child environment, so the mpi env repeats the
setting.

Measured on glibc 2.42: a single later dlopen can claim 1664 bytes of
initial-exec TLS, so 2.6.0 clears it by about seven bytes once alignment is
counted. That thin margin is why the same build fails in one environment and
passes in another. Preloading libgomp also works, but only for libgomp,
since it removes it from the queue without enlarging the surplus.

Reported upstream as Qiskit/qiskit#17213; revert this once that is fixed.

Assisted-by: Claude Opus 5
@garrison
Jim Garrison (garrison) merged commit bae8c33 into stable/1.7 Oct 16, 2026
17 checks passed
@garrison
Jim Garrison (garrison) deleted the mergify/bp/stable/1.7/pr-62 branch October 16, 2026 19:53
@garrison Jim Garrison (garrison) added this to the 1.7.1 milestone Oct 16, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants