From 90afb07d7a383cda15ccc34dea557dac69e4c57e Mon Sep 17 00:00:00 2001 From: Jim Garrison Date: Fri, 9 Oct 2026 16:44:25 -0400 Subject: [PATCH 1/6] Document how to retrieve TPB amplitudes, and the bitstring packing order (#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 --- README.md | 34 +++++++++ python/__init__.py | 74 +++++++++++++++++-- python/bindings.cpp | 7 +- ...tpb-amplitude-output-2d9f4a1c76e8b305.yaml | 20 +++++ 4 files changed, 128 insertions(+), 7 deletions(-) create mode 100644 releasenotes/notes/document-tpb-amplitude-output-2d9f4a1c76e8b305.yaml diff --git a/README.md b/README.md index c8f46c3..fbd9d28 100644 --- a/README.md +++ b/README.md @@ -224,6 +224,7 @@ config = sbd.TPB_SBD() | `adet_comm_size` | 1 | Alpha determinant communicator size | | `bdet_comm_size` | 1 | Beta determinant communicator size | | `task_comm_size` | 1 | Task communicator size | +| `dump_matrix_form_wf` | `""` | Path to write the amplitudes to — see [Getting the TPB amplitudes](#getting-the-tpb-amplitudes) for the file format | Total MPI ranks = `task_comm_size × adet_comm_size × bdet_comm_size`. @@ -267,6 +268,39 @@ results = sbd.tpb_diag(fcidump, adet, bdet, sbd_data, **Returns:** `dict` with keys `energy`, `density`, `carryover_adet`, `carryover_bdet`, `one_p_rdm`, `two_p_rdm`. +#### Getting the TPB amplitudes + +Neither return dict holds the wavefunction amplitudes. Two separate mechanisms write them, in **different formats**, and the one named in the function signature is not the one most callers want: + +| | `sbd_data.dump_matrix_form_wf` | `savename` argument | +|---|---|---| +| Writes | the full `len(adet) × len(bdet)` subspace, gathered onto one rank | one file per b_comm position, each holding **only that rank's block** | +| Path | exactly what you set | `f"{savename}{rank_b:06d}.bin"` | +| Header | **none** | three `size_t` (`adet_range`, `bdet_range`, `det_length`), then that block's alpha and beta determinant words | +| Payload | `float64`, row-major over α×β | `adet_range × bdet_range` `float64` | +| Intended for | analysis | checkpoint/restart, paired with `loadname` | + +So the whole-subspace read is just: + +```python +config.dump_matrix_form_wf = "wf.bin" +results = sbd.tpb_diag(fcidump, adet, bdet, config) + +amplitudes = np.fromfile("wf.bin").reshape(len(adet), len(bdet)) +``` + +Naming that file `.dat` or `.txt` instead writes a text table labelling every amplitude with its own alpha and beta bitstring — the easiest form to inspect by hand, and it sidesteps having to reconstruct the ordering yourself: + +``` + # : : +``` + +Be aware that this text form is written with the stream's default formatting, so amplitudes land with only ~6 significant digits; use the binary form when you need full `float64` precision. + +Rows and columns follow `adet`/`bdet` as the diagonalizer saw them. `sort_bitarray` removes duplicates as well as sorting, so pass lists already in that form if you intend to label the axes yourself — `from_strings` documents the word packing. + +This is the mechanism [`sbd_solver.py`](python/sbd_solver.py) uses to populate `SCIState.amplitudes` — though it skips the write entirely when SBD is picking the carryover itself (`carryover_type` 1, against a qiskit-addon-sqd that accepts it), since then the next subspace comes back in `carryover_adet`/`carryover_bdet` and the amplitudes never need to reach Python. + ```python # GDB: over an explicit list of full determinants rather than a product space results = sbd.gdb_diag(fcidump, det, sbd_data, diff --git a/python/__init__.py b/python/__init__.py index f8ccab8..fce06f9 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -526,13 +526,23 @@ def LoadAlphaDets(filename, bit_length, total_bit_length, device=None): def makestring(config, bit_length, total_bit_length, device=None): - """Convert determinant to string representation.""" + """Render a packed determinant as a ``total_bit_length``-character string. + + The inverse of :func:`from_string`, and subject to the same right-to-left + word packing, which :func:`from_strings` describes. + """ _ensure_initialized() return get_backend(device).makestring(config, bit_length, total_bit_length) def from_string(s, bit_length, total_bit_length, device=None): - """Convert binary string to determinant format.""" + """Convert a binary string to SBD's packed determinant format. + + Returns a list of ``size_t`` words, e.g. ``from_string("0011", 20, 4)`` is + ``[3]``. :func:`makestring` is the inverse. See :func:`from_strings`, the + bulk form, for the word-packing convention and for why it is worth + preferring that form outside of one-off calls. + """ _ensure_initialized() return get_backend(device).from_string(s, bit_length, total_bit_length) @@ -544,6 +554,25 @@ def from_strings(strings, bit_length, total_bit_length, device=None): determinant across the Python boundary. Prefer it for anything larger than a handful: the per-call form costs tens of thousands of round trips on the vendored inputs alone. + + Each row holds ``total_bit_length`` bits packed ``bit_length`` to a word, so + a bitstring longer than ``bit_length`` occupies more than one word. The + packing runs from the **right end** of the string: its trailing + ``bit_length`` characters become word 0, the next ``bit_length`` to the left + become word 1, and so on, so the leading characters of a multi-word string + land in the *last* word rather than the first:: + + from_strings(["0011"], 20, 4) # -> array([[3]], dtype=uint64) + from_strings(["1" + "0" * 69], 63, 70) # -> array([[ 0, 64]], ...) + from_strings(["0" * 69 + "1"], 63, 70) # -> array([[1, 0]], ...) + + Within a word the rightmost character is the least significant bit, which is + why ``"0011"`` packs to ``3`` rather than ``12``. Getting this backwards + produces a valid-looking but wrong subspace rather than an error, so check + round trips with :func:`makestring` when building determinants by hand. + + ``bit_length`` must match the ``bit_length`` of the config passed to the + diagonalizers, since it sets how many words each determinant occupies. """ _ensure_initialized() return get_backend(device).from_strings(list(strings), bit_length, total_bit_length) @@ -579,13 +608,17 @@ def tpb_diag_from_files(fcidumpfile, adetfile, sbd_data, fcidumpfile: Path to FCIDUMP file. adetfile: Path to alpha determinants file. sbd_data: TPB_SBD configuration object. - loadname: Path to load initial wavefunction (optional). - savename: Path to save final wavefunction (optional). + loadname: Path prefix to load an initial wavefunction from (optional). + savename: Path prefix to save the final wavefunction to (optional). device: Override device ('cpu', 'gpu', or None for default). Returns: dict with keys: energy, density, carryover_adet, carryover_bdet, one_p_rdm, two_p_rdm. + + Note: + As with :func:`tpb_diag`, no amplitudes are returned in the dict; see + that function for the two ways of writing them to disk. """ _ensure_initialized() backend = get_backend(device) @@ -604,13 +637,42 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, adet: Alpha determinants. bdet: Beta determinants. sbd_data: TPB_SBD configuration object. - loadname: Path to load initial wavefunction (optional). - savename: Path to save final wavefunction (optional). + loadname: Path prefix to load an initial wavefunction from (optional). + savename: Path prefix to save the final wavefunction to (optional). See + the note below: for reading amplitudes back, + ``sbd_data.dump_matrix_form_wf`` is usually the one you want. device: Override device ('cpu', 'gpu', or None for default). Returns: dict with keys: energy, density, carryover_adet, carryover_bdet, one_p_rdm, two_p_rdm. + + Note: + The returned dict holds no wavefunction amplitudes. There are two ways + to get them, and they write *different* formats -- do not mix them up: + + ``sbd_data.dump_matrix_form_wf = "wf.bin"`` writes the full + ``len(adet) x len(bdet)`` subspace, gathered onto one rank, as a bare + row-major array of ``float64`` with **no header** -- so + ``np.fromfile(path).reshape(len(adet), len(bdet))`` reads it. Naming it + ``.dat`` or ``.txt`` instead writes a text table in which every + amplitude is labelled with its own alpha and beta bitstring, which is + the easiest form to inspect by hand. This is the mechanism this + package's own SQD integration uses. + + ``savename`` instead writes one file per b_comm position, + ``f"{savename}{rank_b:06d}.bin"``, each holding only **that rank's + block** of the subspace: three ``size_t`` headers + (``adet_range``, ``bdet_range``, ``det_length``), then that block's + alpha and then beta determinant words, then + ``adet_range x bdet_range`` ``float64`` amplitudes. It pairs with + ``loadname`` for restarting a run, so it is the better choice for + checkpointing and the worse one for analysis. + + Amplitudes are ordered by ``adet``/``bdet`` as the diagonalizer saw + them. ``sort_bitarray`` removes duplicates as well as sorting, so pass + lists already in that form if you need to label the rows and columns + yourself. """ _ensure_initialized() backend = get_backend(device) diff --git a/python/bindings.cpp b/python/bindings.cpp index 19fbb4b..266d3a6 100644 --- a/python/bindings.cpp +++ b/python/bindings.cpp @@ -319,7 +319,12 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { .def_readwrite("bit_length", &sbd::tpb::SBD::bit_length, "Bit length for determinant representation") .def_readwrite("dump_matrix_form_wf", &sbd::tpb::SBD::dump_matrix_form_wf, - "Filename to dump wavefunction in matrix form") + "Path to write the full adet x bdet amplitude matrix to, " + "gathered onto one rank. A .bin path writes a bare " + "row-major float64 array with no header; .dat/.txt write " + "a labelled text table. Distinct from tpb_diag's " + "savename, which writes per-rank blocks with a header " + "for restart -- see the tpb_diag docstring.") #ifdef SBD_THRUST .def_readwrite("use_precalculated_dets", &sbd::tpb::SBD::use_precalculated_dets, "Use precalculated determinants (THRUST)") diff --git a/releasenotes/notes/document-tpb-amplitude-output-2d9f4a1c76e8b305.yaml b/releasenotes/notes/document-tpb-amplitude-output-2d9f4a1c76e8b305.yaml new file mode 100644 index 0000000..c45b61e --- /dev/null +++ b/releasenotes/notes/document-tpb-amplitude-output-2d9f4a1c76e8b305.yaml @@ -0,0 +1,20 @@ +--- +other: + - | + TPB's two wavefunction output paths are now documented, including the fact + that they write different formats. ``sbd_data.dump_matrix_form_wf`` writes + the full ``len(adet) x len(bdet)`` subspace, gathered onto one rank, as a + bare row-major ``float64`` array with no header, so + ``np.fromfile(path).reshape(len(adet), len(bdet))`` reads it; naming the + file ``.dat`` or ``.txt`` writes a text table labelling each amplitude with + its own alpha and beta bitstring, at the stream's default precision rather + than full ``float64``. The ``savename`` argument of ``tpb_diag`` instead + writes one file per b_comm position holding only that rank's block, behind a + three-``size_t`` header and that block's determinant words, and is meant to + pair with ``loadname`` for restart. For reading amplitudes back, + ``dump_matrix_form_wf`` is normally the one you want. + + The packing convention of ``from_string`` and ``makestring`` is documented + as well. Bitstrings are packed from the right: the trailing ``bit_length`` + characters become word 0, so the leading characters of a string longer than + one word land in the last word rather than the first. From b49c22b54813be1e0344ba52f8ffd03e3825433e Mon Sep 17 00:00:00 2001 From: Jim Garrison Date: Mon, 12 Oct 2026 09:33:59 -0400 Subject: [PATCH 2/6] macOS: find libomp in the environment being built into (#58) 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 --- ...-openmp-build-prefix-5c97decf03d0c7cc.yaml | 12 +++++++ setup.py | 31 ++++++++++++++----- 2 files changed, 35 insertions(+), 8 deletions(-) create mode 100644 releasenotes/notes/darwin-openmp-build-prefix-5c97decf03d0c7cc.yaml diff --git a/releasenotes/notes/darwin-openmp-build-prefix-5c97decf03d0c7cc.yaml b/releasenotes/notes/darwin-openmp-build-prefix-5c97decf03d0c7cc.yaml new file mode 100644 index 0000000..6eb0576 --- /dev/null +++ b/releasenotes/notes/darwin-openmp-build-prefix-5c97decf03d0c7cc.yaml @@ -0,0 +1,12 @@ +--- +fixes: + - | + Building on macOS inside a conda-build or ``rattler-build`` recipe no longer + fails with ``no OpenMP runtime found on this macOS host`` when the recipe + provides ``llvm-openmp``. The build located ``libomp`` by probing + ``$CONDA_PREFIX``, which under those tools names the base installation + rather than the environment being built into, so the ``omp.h`` installed in + the build prefix went unexamined and the build fell through to Homebrew. + ``sys.prefix`` and ``sys.base_prefix`` are now consulted as well, after + ``CONDA_PREFIX`` and before Homebrew, so an ordinary activated environment + behaves as before. diff --git a/setup.py b/setup.py index 394d1a9..3993594 100644 --- a/setup.py +++ b/setup.py @@ -774,16 +774,26 @@ def detect_gpu_toolchain(): print("\nConfiguring CPU backend (_core_cpu)") if platform.system() == 'Darwin': # macOS has no system OpenMP, so libomp comes from a package manager. - # Prefer the conda env when it has one: those are the libraries actually + # Prefer an environment that has one: those are the libraries actually # LOADED at import time (resolved via the python executable's # @loader_path/../lib), so building against Homebrew's copies instead # means compiling against different libraries than the process runs on. - conda_prefix = os.environ.get('CONDA_PREFIX') - if conda_prefix and os.path.exists( - os.path.join(conda_prefix, 'include', 'omp.h')): - omp_inc = os.path.join(conda_prefix, 'include') - omp_lib = openblas_lib = os.path.join(conda_prefix, 'lib') - print(f"Darwin: libomp and BLAS from conda env {conda_prefix}") + # + # CONDA_PREFIX names the *activated* environment, which is not always + # the one being built into: conda-build and rattler-build activate their + # own base installation and install into $PREFIX, whose python is the + # interpreter running this script. So fall back to sys.prefix, which + # identifies the environment the extension is actually built for -- the + # same convention _mpi_prefix_from_env_prefix() already uses. + omp_prefix = next( + (p for p in (os.environ.get('CONDA_PREFIX'), sys.prefix, + getattr(sys, 'base_prefix', None)) + if p and os.path.exists(os.path.join(p, 'include', 'omp.h'))), + None) + if omp_prefix: + omp_inc = os.path.join(omp_prefix, 'include') + omp_lib = openblas_lib = os.path.join(omp_prefix, 'lib') + print(f"Darwin: libomp and BLAS from env prefix {omp_prefix}") else: # Ask brew for its prefix rather than assuming: it is /opt/homebrew # on Apple silicon and /usr/local on Intel, so either one hardcoded @@ -797,8 +807,13 @@ def detect_gpu_toolchain(): if not os.path.exists(os.path.join(omp_inc, 'omp.h')): # Fail here with the fix, rather than 100 lines later with # "'omp.h' file not found" from the middle of a compile. + _searched = ', '.join(dict.fromkeys( + [os.path.join(p, 'include') + for p in (os.environ.get('CONDA_PREFIX'), sys.prefix, + getattr(sys, 'base_prefix', None)) if p] + + [omp_inc])) print("Error: no OpenMP runtime found on this macOS host.\n" - f" Looked in $CONDA_PREFIX/include and {omp_inc}.\n" + f" Looked in: {_searched}.\n" " Apple clang ships without OpenMP, so install one:\n" " conda install -c conda-forge llvm-openmp (preferred)\n" " brew install libomp") From b3bf3fa1837c9154227928d5ee9c3df03c3a7c9c Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Wed, 14 Oct 2026 16:24:30 -0400 Subject: [PATCH 3/6] docs: correct outdated and misleading API docstrings (#60) * 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 * 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 --------- Co-authored-by: Claude Opus 5.5 --- docs/apidocs/sbd.rst | 9 ++-- python/__init__.py | 92 +++++++++++++++++++++++++++++++++-------- python/device_config.py | 36 ++++++---------- python/sbd_solver.py | 21 ++++------ 4 files changed, 101 insertions(+), 57 deletions(-) diff --git a/docs/apidocs/sbd.rst b/docs/apidocs/sbd.rst index 1a5e968..7484075 100644 --- a/docs/apidocs/sbd.rst +++ b/docs/apidocs/sbd.rst @@ -67,10 +67,11 @@ Fields shared by ``TPB_SBD`` and ``GDB_SBD``: reads differs, and SBD silently ignores the other. * - ``bit_length`` - Bits per word when packing determinants. At most 63; 64 is undefined - behavior in SBD's multi-rank redistribution. For GDB's heatbath expansion - (``carryover_type`` 2 or 3) it must also be even whenever a determinant spans - more than one word, or the expansion crashes or returns a wrong energy. The - ``examples/gdb`` drivers require an even value of at most 62. + behavior in SBD's multi-rank redistribution. For GDB it must also be even + whenever a determinant spans more than one word: an odd value breaks SBD's + alpha/beta conversion, which can crash or return a wrong energy, in plain + diagonalization as well as in the heatbath expansion. The ``examples/gdb`` + drivers require an even value of at most 62. Fields of ``TPB_SBD`` only: diff --git a/python/__init__.py b/python/__init__.py index fce06f9..d0a9b78 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -318,8 +318,9 @@ def init(device='cpu', comm_backend='mpi'): first use with ``device='cpu'`` and ``comm_backend='mpi'``. Call it explicitly only when you need GPU or want to control startup timing. - The device can be overridden per-call via the ``device`` parameter on - ``tpb_diag()``, ``tpb_diag_from_files()``, and ``get_backend()``. + The device can be overridden per call via the ``device`` parameter of + ``tpb_diag()``, ``tpb_diag_from_files()``, ``gdb_diag()`` and + ``get_backend()``. Args: device: Default compute device — 'cpu', 'gpu', 'gpu-omp', or 'auto'. @@ -413,11 +414,20 @@ def get_backend(device=None): Get the backend module for the given device. Auto-initializes SBD if needed. Passing ``device`` overrides the - default — this is how you switch between CPU and GPU within the - same process. + default, which is how you switch backends within one process. The one + combination to avoid is the CPU backend together with the OpenMP-offload + one ('gpu-omp'): they share an OpenMP runtime, and loading both leaves + offload regions running on the host. ``sbd.has_backend_conflict()`` + reports it. Args: - device: 'cpu', 'gpu', 'auto', or None (use default). + device: 'cpu', 'gpu' (Thrust), 'gpu-omp' (OpenMP offload), 'auto', an + alias accepted by :func:`init`, or None. None means the device set + by :func:`init`. If nothing has initialized yet it means 'auto', + which prefers Thrust, then OpenMP offload, then CPU, whereas the + first call of :func:`tpb_diag` or :func:`gdb_diag` initializes with + 'cpu'. Call :func:`init` with the device you want to avoid + depending on call order. Returns: The pybind11 backend module (_core_cpu, _core_gpu_thrust, or @@ -496,13 +506,19 @@ def barrier(): # --------------------------------------------------------------------------- def TPB_SBD(device=None): - """Create TPB_SBD configuration object.""" + """Create a TPB_SBD configuration object, for :func:`tpb_diag`. + + Its fields are listed in the tables below. + """ _ensure_initialized() return get_backend(device).TPB_SBD() def GDB_SBD(device=None): - """Create GDB_SBD configuration object.""" + """Create a GDB_SBD configuration object, for :func:`gdb_diag`. + + Its fields are listed in the tables below. + """ _ensure_initialized() return get_backend(device).GDB_SBD() @@ -514,13 +530,43 @@ def FCIDump(device=None): def LoadFCIDump(filename, device=None): - """Load FCIDUMP file.""" + """Load an FCIDUMP file. + + Args: + filename: Path to the FCIDUMP file. + device: Backend whose ``FCIDump`` type to build; None uses the default. + + Returns: + The ``FCIDump`` object that :func:`tpb_diag` and :func:`gdb_diag` take. + """ _ensure_initialized() return get_backend(device).LoadFCIDump(filename) def LoadAlphaDets(filename, bit_length, total_bit_length, device=None): - """Load alpha determinants from file.""" + """Load determinants from a file of bitstrings, one per line. + + Each line is packed as :func:`from_strings` does. Upstream sorts the + determinants and removes duplicates, so the result is ready to use as a + TPB determinant list and may hold fewer determinants than the file has + lines. + + The format is chosen by extension: ``.txt`` for one bitstring per line, + ``.bin`` for upstream's binary format. Any other extension returns an + empty list without an error. + + Args: + filename: Path to the file. + bit_length: Bits per packed word. Must match the config's + ``bit_length``. + total_bit_length: Length of each bitstring, ``norb`` for the alpha or + beta strings of TPB. + device: Backend to use; None uses the default. + + Returns: + A list of packed determinants, each a list of words, sorted and + distinct. + """ _ensure_initialized() return get_backend(device).LoadAlphaDets(filename, bit_length, total_bit_length) @@ -550,7 +596,7 @@ def from_string(s, bit_length, total_bit_length, device=None): def from_strings(strings, bit_length, total_bit_length, device=None): """Pack many bitstrings into one ``(n, words)`` array. - The bulk form of :func:`from_string`, looping in C++ instead of once per + The bulk form of ``sbd.from_string``, looping in C++ instead of once per determinant across the Python boundary. Prefer it for anything larger than a handful: the per-call form costs tens of thousands of round trips on the vendored inputs alone. @@ -563,13 +609,13 @@ def from_strings(strings, bit_length, total_bit_length, device=None): land in the *last* word rather than the first:: from_strings(["0011"], 20, 4) # -> array([[3]], dtype=uint64) - from_strings(["1" + "0" * 69], 63, 70) # -> array([[ 0, 64]], ...) - from_strings(["0" * 69 + "1"], 63, 70) # -> array([[1, 0]], ...) + from_strings(["1" + "0" * 69], 62, 70) # -> array([[ 0, 128]], ...) + from_strings(["0" * 69 + "1"], 62, 70) # -> array([[1, 0]], ...) Within a word the rightmost character is the least significant bit, which is why ``"0011"`` packs to ``3`` rather than ``12``. Getting this backwards produces a valid-looking but wrong subspace rather than an error, so check - round trips with :func:`makestring` when building determinants by hand. + round trips with ``sbd.makestring`` when building determinants by hand. ``bit_length`` must match the ``bit_length`` of the config passed to the diagonalizers, since it sets how many words each determinant occupies. @@ -604,13 +650,19 @@ def tpb_diag_from_files(fcidumpfile, adetfile, sbd_data, """ Perform TPB diagonalization from files. + The one determinant file supplies both the alpha and the beta list, so the + subspace is the product of that list with itself. For different alpha and + beta lists, as an open-shell case needs, load them and call :func:`tpb_diag`. + Args: fcidumpfile: Path to FCIDUMP file. - adetfile: Path to alpha determinants file. + adetfile: Path to the determinant file, used for both alpha and beta. sbd_data: TPB_SBD configuration object. loadname: Path prefix to load an initial wavefunction from (optional). savename: Path prefix to save the final wavefunction to (optional). - device: Override device ('cpu', 'gpu', or None for default). + device: Backend to run on: 'cpu', 'gpu' (Thrust), 'gpu-omp' (OpenMP + offload), or an alias accepted by :func:`init`. None uses the + default; see :func:`get_backend`. Returns: dict with keys: energy, density, carryover_adet, carryover_bdet, @@ -641,7 +693,9 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, savename: Path prefix to save the final wavefunction to (optional). See the note below: for reading amplitudes back, ``sbd_data.dump_matrix_form_wf`` is usually the one you want. - device: Override device ('cpu', 'gpu', or None for default). + device: Backend to run on: 'cpu', 'gpu' (Thrust), 'gpu-omp' (OpenMP + offload), or an alias accepted by :func:`init`. None uses the + default; see :func:`get_backend`. Returns: dict with keys: energy, density, carryover_adet, carryover_bdet, @@ -733,7 +787,11 @@ def gdb_diag(fcidump, det, sbd_data, writes one file per b_comm position, ``f"{savename}{rank_b:06d}.bin"``, each holding only that shard; rank 0's file is not the whole wavefunction. GDB has no combined matrix-form dump. - device: Override device ('cpu', 'gpu', or None for default). + device: Backend to run on: 'cpu', 'gpu' (Thrust), 'gpu-omp' (OpenMP + offload), or an alias accepted by :func:`init`. None uses the + default; see :func:`get_backend`. 'gpu-omp' is accepted, but SBD + has no OpenMP-offload GDB code, so the run executes on the host; + use 'gpu' for GDB on GPUs. determinant_distribution: How to place determinants across b_comm — one of ``input`` (keep the shards as given), ``equal-bra-a`` (default; equal distinct-alpha count per rank, which is what balances the diff --git a/python/device_config.py b/python/device_config.py index f68aa0b..477c53a 100644 --- a/python/device_config.py +++ b/python/device_config.py @@ -27,18 +27,14 @@ class DeviceConfig: Usage:: - # Auto-detect (uses GPU if available) - config = DeviceConfig.auto() + from sbd.device_config import DeviceConfig + from sbd.sbd_solver import solve_sci_batch - # Force CPU - config = DeviceConfig.cpu() + device = DeviceConfig.auto() # best backend that is built + device = DeviceConfig.cpu() # or force one + device = DeviceConfig.gpu(max_memory_gb=16) - # Force GPU with specific settings - config = DeviceConfig.gpu(max_memory_gb=16) - - # Apply to SBD configuration - sbd_config = sbd.TPB_SBD() - config.apply(sbd_config) + results = solve_sci_batch(..., device_config=device) """ def __init__(self, device: str = 'cpu', @@ -72,17 +68,11 @@ def auto(cls, max_memory_gb: int = -1) -> 'DeviceConfig': """ Auto-detect the best available backend and use it. - Resolves against the backends that were actually COMPILED, not just the - hardware that is present. Detecting a GPU and returning 'gpu' - unconditionally was wrong in two ways: on an AMD host it selected the - Thrust/CUDA backend, which cannot exist there (upstream wires Thrust to - nvc++ -cuda), so a machine that reported "GPU detected (HIP)" then failed - to load a CUDA module; and on a CPU-only build with a GPU present it - picked a backend that was never built. - - Preference order matches sbd._resolve_device(): Thrust ('gpu') first - where it exists, since it is the long-validated NVIDIA default and keeps - more phases on the device, then OpenMP offload ('gpu-omp'), then CPU. + Chooses among the backends that were actually built, not merely the + hardware present, so an AMD host or a CPU-only build is never handed a + backend it lacks. The preference order is Thrust ('gpu') where it exists, + then OpenMP offload ('gpu-omp'), then CPU, the same order + ``sbd.get_backend('auto')`` uses. Args: max_memory_gb: Maximum GPU memory in GB (-1 = auto) @@ -158,8 +148,8 @@ def gpu_omp(cls, max_memory_gb: int = -1) -> 'DeviceConfig': one combination to avoid in a single process is this backend together with the CPU one: they share an OpenMP runtime (``libnvomp`` on NVIDIA, ``libomp`` on AMD), and loading ``_core_cpu`` first leaves it initialised - host-only, after which offload regions run on the host. See - :func:`sbd.has_backend_conflict`. + host-only, after which offload regions run on the host. + ``sbd.has_backend_conflict()`` reports that combination. """ return cls(device='gpu-omp', max_memory_gb=max_memory_gb) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index e8b4508..b346158 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -470,25 +470,20 @@ def assemble_rdms(results: dict, norb: int) -> tuple[np.ndarray | None, np.ndarr which only computes the diagonal density, not the correlation functions), matching SCIResult's own ``rdm1``/``rdm2`` default. - The reshape/transpose below is not a guess: it was verified against - PySCF's ``make_rdm1``/``make_rdm2`` on a fixed subspace, on all three - SBD backends (cpu, gpu-thrust, gpu-omp-offload) -- both element-wise on - the full tensors and via the energy identity + The layout conversion was verified against PySCF's ``make_rdm1`` and + ``make_rdm2`` on all three backends (cpu, gpu, gpu-omp), both element-wise + and through the energy identity ``E = einsum("pr,pr->",rdm1,hcore) + 0.5*einsum("prqs,prqs->",rdm2,eri)`` - that ``SCIResult.rdm1``/``rdm2`` are contracted with everywhere else in - qiskit-addon-sqd (e.g. ``fermion.py``'s own ``solve_fermion``). + with which ``qiskit-addon-sqd`` contracts ``SCIResult.rdm1`` and ``rdm2``. - SBD's documented layout (sbd-ext docs/user-guide.md, matching the C++ - reference in apps/chemistry_tpb_selected_basis_diagonalization/main.cc):: + SBD's layout, as documented in upstream SBD's ``docs/user-guide.md``:: one_p_rdm[s][i + L*j] = two_p_rdm[s+2t][i + L*j + L^2*k + L^3*l] = - A Fortran-order reshape implements those flat-index formulas directly - (arr_F[i, j] / arr_F[i, j, k, l]); rdm1 needs no further transpose - (it is symmetric here regardless), and rdm2's spin-summed block sum - needs axes (0, 2, 1, 3) to land in the "prqs" slot order SCIResult's - contract expects. + A Fortran-order reshape implements those flat-index formulas directly. + rdm1 needs no further transpose, and the spin-summed rdm2 is transposed + with axes ``(0, 2, 1, 3)`` into the "prqs" order ``SCIResult`` expects. """ one_p_rdm = results.get("one_p_rdm") two_p_rdm = results.get("two_p_rdm") From 73c71a4a09abd77a6ff35da7adab7fbc308594dc Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Fri, 16 Oct 2026 09:52:00 -0400 Subject: [PATCH 4/6] Give SCIState its amplitudes in PySCF's sign convention (#62) * 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 * Apply batched suggestions from code review Co-authored-by: Jim Garrison --------- Co-authored-by: Claude Opus 5.5 Co-authored-by: Jim Garrison Co-authored-by: Jim Garrison --- python/__init__.py | 10 ++ python/sbd_solver.py | 39 +++++- ...tate-sign-convention-58b57b05535e1ca2.yaml | 15 +++ test/test_sci_state.py | 126 ++++++++++++++++++ 4 files changed, 189 insertions(+), 1 deletion(-) create mode 100644 releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml create mode 100644 test/test_sci_state.py diff --git a/python/__init__.py b/python/__init__.py index d0a9b78..e0f1c7e 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -727,6 +727,16 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, them. ``sort_bitarray`` removes duplicates as well as sorting, so pass lists already in that form if you need to label the rows and columns yourself. + + Both files use SBD's sign convention, which orders a determinant's + creation operators orbital by orbital (alpha 0, beta 0, alpha 1, ...). + PySCF, and so qiskit-addon-sqd's ``SCIState``, puts all alpha operators + before all beta ones, so some determinants differ in sign: multiply each + amplitude by ``(-1) ** n``, where ``n`` counts, over the occupied beta + orbitals ``p``, the occupied alpha orbitals above ``p``, before handing + it to PySCF. Magnitudes, and so selection by weight, are the same either + way. :func:`sbd.sbd_solver.solve_sci` converts the amplitudes in its + ``SCIState`` already. """ _ensure_initialized() backend = get_backend(device) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index b346158..59f975f 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -297,8 +297,10 @@ def _solve_sci_core( "a different size means the dump and the subspace have diverged." ) + amplitudes = flat.reshape(n_a, n_b) + _to_alpha_first_signs(amplitudes, solved_strings_a, solved_strings_b, norb) sci_state = SCIState( - amplitudes=flat.reshape(n_a, n_b), + amplitudes=amplitudes, ci_strs_a=solved_strings_a, ci_strs_b=solved_strings_b, norb=norb, @@ -688,6 +690,41 @@ def _ci_strings_to_sbd_dets( return backend.sort_bitarray(dets) +def _to_alpha_first_signs(amplitudes, strings_a, strings_b, norb, block_elements=1 << 24): + """Convert SBD's amplitudes, in place, to the sign convention SCIState assumes. + + The two conventions order the creation operators of a determinant differently. + SBD interleaves the spins orbital by orbital (alpha 0, beta 0, alpha 1, ...); + PySCF, and so ``SCIState.rdm()`` and ``SCIState.spin_square()``, put every alpha + operator before every beta one. Moving each beta operator past the alpha + operators of higher orbitals gives the determinant the phase + + (-1) ** sum over occupied beta orbitals p of (occupied alpha orbitals above p) + + The magnitudes are the same, so energies, SBD's own RDMs, occupancies and any + selection by amplitude magnitude are unaffected, which is why this went unnoticed: + only what PySCF computes *from the amplitudes* was wrong -- e.g. ``spin_square()`` + of 0.036 for an exact singlet, and an energy from ``rdm()`` 89 mHa off. + + Args: + amplitudes: The ``(len(strings_a), len(strings_b))`` matrix, modified in place. + strings_a: Alpha CI strings labelling the rows; bit ``p`` is orbital ``p``. + strings_b: Beta CI strings labelling the columns. + norb: Number of spatial orbitals. + block_elements: Rows are converted in blocks of about this many elements, so + the temporary sign matrix stays bounded at large subspaces. + """ + orbitals = np.arange(norb, dtype=np.int64) + occupied_b = ((np.asarray(strings_b, dtype=np.int64)[:, None] >> orbitals) & 1).astype(np.int32) + occupied_a = (np.asarray(strings_a, dtype=np.int64)[:, None] >> orbitals) & 1 + # alpha electrons in orbitals strictly above p, kept only as a parity + above_a = ((np.cumsum(occupied_a[:, ::-1], axis=1)[:, ::-1] - occupied_a) & 1).astype(np.int32) + rows = max(1, block_elements // max(1, len(strings_b))) + for start in range(0, len(strings_a), rows): + parity = (above_a[start:start + rows] @ occupied_b.T) & 1 + amplitudes[start:start + rows][parity == 1] *= -1.0 + + def _sbd_dets_to_ci_strings( dets: list[list[int]], norb: int, backend, bit_length: int = SBD_DEFAULT_BIT_LENGTH, diff --git a/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml new file mode 100644 index 0000000..c32ec5c --- /dev/null +++ b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml @@ -0,0 +1,15 @@ +--- +fixes: + - | + ``SCIResult.sci_state`` returned by :func:`sbd.sbd_solver.solve_sci` and + :func:`sbd.sbd_solver.solve_sci_batch` now holds amplitudes in the sign + convention ``qiskit-addon-sqd``'s ``SCIState`` assumes. SBD orders a + determinant's creation operators orbital by orbital, PySCF puts all alpha + operators before all beta ones, so some determinants carried the opposite + sign. Everything ``SCIState`` computes through PySCF from the amplitudes + was wrong as a result: ``spin_square()`` returned 0.036 for an exact + singlet, and an energy rebuilt from ``rdm()`` was 89 mHa off. Energies, + ``SCIResult.rdm1``/``rdm2`` (SBD's own RDMs), orbital occupancies and + carryover selection, which uses only amplitude magnitudes, were not + affected. The raw wavefunction files SBD writes keep its own convention; + the :func:`sbd.tpb_diag` documentation gives the conversion. diff --git a/test/test_sci_state.py b/test/test_sci_state.py new file mode 100644 index 0000000..1f39aac --- /dev/null +++ b/test/test_sci_state.py @@ -0,0 +1,126 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""``SCIResult.sci_state`` from ``solve_sci`` must be the state PySCF would return. + +``SCIState`` hands its amplitudes to PySCF for ``rdm()`` and ``spin_square()``, so +they have to follow PySCF's sign convention, in which every alpha creation operator +precedes every beta one. SBD orders them orbital by orbital instead, which gives +some determinants the opposite sign. Energies, SBD's own RDMs (``SCIResult.rdm1`` / +``rdm2``), occupancies and selection by amplitude magnitude never see the +difference, which is why only these comparisons catch it. + +Each case solves one subspace twice -- with SBD through ``solve_sci``, and with +qiskit-addon-sqd's PySCF solver -- and compares the states, not just the energies. +The subspaces are full spaces and seeded random truncations of them: a small +truncation can happen to need no sign change at all, so a lone small case can pass +by accident. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +pyscf = pytest.importorskip("pyscf", reason="pyscf builds the reference Hamiltonians") +fermion = pytest.importorskip("qiskit_addon_sqd.fermion", + reason="qiskit-addon-sqd provides SCIState and the reference solver") + +from pyscf import ao2mo, gto, scf # noqa: E402 + +from sbd.sbd_solver import solve_sci # noqa: E402 + +SYSTEMS = { + # (atom, basis, 2S): a doublet, a singlet, and one with more orbitals + "h3-doublet": ("H 0 0 0; H 0 0 0.8; H 0 0 2.2", "sto-3g", 1), + "h4-singlet": ("H 0 0 0; H 0 0 0.9; H 0 0 1.9; H 0 0 3.0", "sto-3g", 0), + "h6-singlet": ("; ".join(f"H 0 0 {0.95 * i:.2f}" for i in range(6)), "sto-3g", 0), +} + + +def _hamiltonian(name): + atom, basis, spin = SYSTEMS[name] + mol = gto.M(atom=atom, basis=basis, spin=spin, verbose=0) + mf = (scf.ROHF if spin else scf.RHF)(mol) + mf.conv_tol = 1e-12 + mf.run() + c = mf.mo_coeff + norb = mol.nao + h1 = c.T @ scf.hf.get_hcore(mol) @ c + eri = ao2mo.restore(1, ao2mo.kernel(mol, c), norb) + return h1, eri, norb, mol.nelec + + +def _strings(norb, n_electrons, keep=None, seed=0): + full = np.array(sorted(sum(1 << o for o in occ) + for occ in itertools.combinations(range(norb), n_electrons))) + if keep is None or keep >= len(full): + return full + chosen = np.random.default_rng(seed).choice(len(full), size=keep, replace=False) + return np.sort(full[np.union1d(chosen, [0])]) # keep Hartree-Fock + + +CASES = [ + ("h3-doublet", None), + ("h4-singlet", None), + ("h6-singlet", None), + ("h6-singlet", 12), +] + + +def _energy_from_rdms(state, h1, eri): + return (np.einsum("pr,pr->", state.rdm(rank=1, spin_summed=True), h1) + + 0.5 * np.einsum("prqs,prqs->", state.rdm(rank=2, spin_summed=True), eri)) + + +def _check_sci_state_matches_pyscf(device_config, name, keep): + h1, eri, norb, nelec = _hamiltonian(name) + strings = (_strings(norb, nelec[0], keep, seed=1), _strings(norb, nelec[1], keep, seed=2)) + + # PySCF's default tolerance leaves the vector good to ~1e-6; converge both tightly + # so that the elementwise comparison measures the convention, not convergence. + expected = fermion.solve_sci(strings, h1, eri, norb=norb, nelec=nelec, tol=1e-14) + got = solve_sci(strings, h1, eri, norb=norb, nelec=nelec, device_config=device_config, + sbd_config={"eps": 1e-12, "max_it": 300}) + from mpi4py import MPI + if MPI.COMM_WORLD.Get_rank() != 0: + return # the amplitudes are gathered onto rank 0 only + + assert got.energy == pytest.approx(expected.energy, abs=1e-8) + assert np.array_equal(got.sci_state.ci_strs_a, expected.sci_state.ci_strs_a) + assert np.array_equal(got.sci_state.ci_strs_b, expected.sci_state.ci_strs_b) + + # The same vector up to one global sign -- the comparison that catches a + # per-determinant sign, where every magnitude-based one passes. + a, b = got.sci_state.amplitudes, expected.sci_state.amplitudes + overlap = np.vdot(b, a) + assert abs(overlap) == pytest.approx(1.0, abs=1e-7) + np.testing.assert_allclose(a * np.sign(overlap), b, atol=1e-6) + + # And what users compute from it. + assert got.sci_state.spin_square() == pytest.approx(expected.sci_state.spin_square(), abs=1e-7) + assert _energy_from_rdms(got.sci_state, h1, eri) == pytest.approx(expected.energy, abs=1e-7) + + +# One body, two wrappers: pytest-mpi skips unmarked tests under --only-mpi and marked +# ones without it (see test_sqd_integration.py). +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_standalone(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep) + + +@pytest.mark.mpi +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_mpi(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep) From 8a09a91a81d7040ec768778147b43006734b4bf0 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Fri, 16 Oct 2026 09:52:00 -0400 Subject: [PATCH 5/6] Give SCIState its amplitudes in PySCF's sign convention (#62) * 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 * Apply batched suggestions from code review Co-authored-by: Jim Garrison --------- Co-authored-by: Claude Opus 5.5 Co-authored-by: Jim Garrison Co-authored-by: Jim Garrison (cherry picked from commit 73c71a4a09abd77a6ff35da7adab7fbc308594dc) --- python/__init__.py | 10 ++ python/sbd_solver.py | 39 +++++- ...tate-sign-convention-58b57b05535e1ca2.yaml | 15 +++ test/test_sci_state.py | 126 ++++++++++++++++++ 4 files changed, 189 insertions(+), 1 deletion(-) create mode 100644 releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml create mode 100644 test/test_sci_state.py diff --git a/python/__init__.py b/python/__init__.py index d0a9b78..e0f1c7e 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -727,6 +727,16 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, them. ``sort_bitarray`` removes duplicates as well as sorting, so pass lists already in that form if you need to label the rows and columns yourself. + + Both files use SBD's sign convention, which orders a determinant's + creation operators orbital by orbital (alpha 0, beta 0, alpha 1, ...). + PySCF, and so qiskit-addon-sqd's ``SCIState``, puts all alpha operators + before all beta ones, so some determinants differ in sign: multiply each + amplitude by ``(-1) ** n``, where ``n`` counts, over the occupied beta + orbitals ``p``, the occupied alpha orbitals above ``p``, before handing + it to PySCF. Magnitudes, and so selection by weight, are the same either + way. :func:`sbd.sbd_solver.solve_sci` converts the amplitudes in its + ``SCIState`` already. """ _ensure_initialized() backend = get_backend(device) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index b346158..59f975f 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -297,8 +297,10 @@ def _solve_sci_core( "a different size means the dump and the subspace have diverged." ) + amplitudes = flat.reshape(n_a, n_b) + _to_alpha_first_signs(amplitudes, solved_strings_a, solved_strings_b, norb) sci_state = SCIState( - amplitudes=flat.reshape(n_a, n_b), + amplitudes=amplitudes, ci_strs_a=solved_strings_a, ci_strs_b=solved_strings_b, norb=norb, @@ -688,6 +690,41 @@ def _ci_strings_to_sbd_dets( return backend.sort_bitarray(dets) +def _to_alpha_first_signs(amplitudes, strings_a, strings_b, norb, block_elements=1 << 24): + """Convert SBD's amplitudes, in place, to the sign convention SCIState assumes. + + The two conventions order the creation operators of a determinant differently. + SBD interleaves the spins orbital by orbital (alpha 0, beta 0, alpha 1, ...); + PySCF, and so ``SCIState.rdm()`` and ``SCIState.spin_square()``, put every alpha + operator before every beta one. Moving each beta operator past the alpha + operators of higher orbitals gives the determinant the phase + + (-1) ** sum over occupied beta orbitals p of (occupied alpha orbitals above p) + + The magnitudes are the same, so energies, SBD's own RDMs, occupancies and any + selection by amplitude magnitude are unaffected, which is why this went unnoticed: + only what PySCF computes *from the amplitudes* was wrong -- e.g. ``spin_square()`` + of 0.036 for an exact singlet, and an energy from ``rdm()`` 89 mHa off. + + Args: + amplitudes: The ``(len(strings_a), len(strings_b))`` matrix, modified in place. + strings_a: Alpha CI strings labelling the rows; bit ``p`` is orbital ``p``. + strings_b: Beta CI strings labelling the columns. + norb: Number of spatial orbitals. + block_elements: Rows are converted in blocks of about this many elements, so + the temporary sign matrix stays bounded at large subspaces. + """ + orbitals = np.arange(norb, dtype=np.int64) + occupied_b = ((np.asarray(strings_b, dtype=np.int64)[:, None] >> orbitals) & 1).astype(np.int32) + occupied_a = (np.asarray(strings_a, dtype=np.int64)[:, None] >> orbitals) & 1 + # alpha electrons in orbitals strictly above p, kept only as a parity + above_a = ((np.cumsum(occupied_a[:, ::-1], axis=1)[:, ::-1] - occupied_a) & 1).astype(np.int32) + rows = max(1, block_elements // max(1, len(strings_b))) + for start in range(0, len(strings_a), rows): + parity = (above_a[start:start + rows] @ occupied_b.T) & 1 + amplitudes[start:start + rows][parity == 1] *= -1.0 + + def _sbd_dets_to_ci_strings( dets: list[list[int]], norb: int, backend, bit_length: int = SBD_DEFAULT_BIT_LENGTH, diff --git a/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml new file mode 100644 index 0000000..c32ec5c --- /dev/null +++ b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml @@ -0,0 +1,15 @@ +--- +fixes: + - | + ``SCIResult.sci_state`` returned by :func:`sbd.sbd_solver.solve_sci` and + :func:`sbd.sbd_solver.solve_sci_batch` now holds amplitudes in the sign + convention ``qiskit-addon-sqd``'s ``SCIState`` assumes. SBD orders a + determinant's creation operators orbital by orbital, PySCF puts all alpha + operators before all beta ones, so some determinants carried the opposite + sign. Everything ``SCIState`` computes through PySCF from the amplitudes + was wrong as a result: ``spin_square()`` returned 0.036 for an exact + singlet, and an energy rebuilt from ``rdm()`` was 89 mHa off. Energies, + ``SCIResult.rdm1``/``rdm2`` (SBD's own RDMs), orbital occupancies and + carryover selection, which uses only amplitude magnitudes, were not + affected. The raw wavefunction files SBD writes keep its own convention; + the :func:`sbd.tpb_diag` documentation gives the conversion. diff --git a/test/test_sci_state.py b/test/test_sci_state.py new file mode 100644 index 0000000..1f39aac --- /dev/null +++ b/test/test_sci_state.py @@ -0,0 +1,126 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""``SCIResult.sci_state`` from ``solve_sci`` must be the state PySCF would return. + +``SCIState`` hands its amplitudes to PySCF for ``rdm()`` and ``spin_square()``, so +they have to follow PySCF's sign convention, in which every alpha creation operator +precedes every beta one. SBD orders them orbital by orbital instead, which gives +some determinants the opposite sign. Energies, SBD's own RDMs (``SCIResult.rdm1`` / +``rdm2``), occupancies and selection by amplitude magnitude never see the +difference, which is why only these comparisons catch it. + +Each case solves one subspace twice -- with SBD through ``solve_sci``, and with +qiskit-addon-sqd's PySCF solver -- and compares the states, not just the energies. +The subspaces are full spaces and seeded random truncations of them: a small +truncation can happen to need no sign change at all, so a lone small case can pass +by accident. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +pyscf = pytest.importorskip("pyscf", reason="pyscf builds the reference Hamiltonians") +fermion = pytest.importorskip("qiskit_addon_sqd.fermion", + reason="qiskit-addon-sqd provides SCIState and the reference solver") + +from pyscf import ao2mo, gto, scf # noqa: E402 + +from sbd.sbd_solver import solve_sci # noqa: E402 + +SYSTEMS = { + # (atom, basis, 2S): a doublet, a singlet, and one with more orbitals + "h3-doublet": ("H 0 0 0; H 0 0 0.8; H 0 0 2.2", "sto-3g", 1), + "h4-singlet": ("H 0 0 0; H 0 0 0.9; H 0 0 1.9; H 0 0 3.0", "sto-3g", 0), + "h6-singlet": ("; ".join(f"H 0 0 {0.95 * i:.2f}" for i in range(6)), "sto-3g", 0), +} + + +def _hamiltonian(name): + atom, basis, spin = SYSTEMS[name] + mol = gto.M(atom=atom, basis=basis, spin=spin, verbose=0) + mf = (scf.ROHF if spin else scf.RHF)(mol) + mf.conv_tol = 1e-12 + mf.run() + c = mf.mo_coeff + norb = mol.nao + h1 = c.T @ scf.hf.get_hcore(mol) @ c + eri = ao2mo.restore(1, ao2mo.kernel(mol, c), norb) + return h1, eri, norb, mol.nelec + + +def _strings(norb, n_electrons, keep=None, seed=0): + full = np.array(sorted(sum(1 << o for o in occ) + for occ in itertools.combinations(range(norb), n_electrons))) + if keep is None or keep >= len(full): + return full + chosen = np.random.default_rng(seed).choice(len(full), size=keep, replace=False) + return np.sort(full[np.union1d(chosen, [0])]) # keep Hartree-Fock + + +CASES = [ + ("h3-doublet", None), + ("h4-singlet", None), + ("h6-singlet", None), + ("h6-singlet", 12), +] + + +def _energy_from_rdms(state, h1, eri): + return (np.einsum("pr,pr->", state.rdm(rank=1, spin_summed=True), h1) + + 0.5 * np.einsum("prqs,prqs->", state.rdm(rank=2, spin_summed=True), eri)) + + +def _check_sci_state_matches_pyscf(device_config, name, keep): + h1, eri, norb, nelec = _hamiltonian(name) + strings = (_strings(norb, nelec[0], keep, seed=1), _strings(norb, nelec[1], keep, seed=2)) + + # PySCF's default tolerance leaves the vector good to ~1e-6; converge both tightly + # so that the elementwise comparison measures the convention, not convergence. + expected = fermion.solve_sci(strings, h1, eri, norb=norb, nelec=nelec, tol=1e-14) + got = solve_sci(strings, h1, eri, norb=norb, nelec=nelec, device_config=device_config, + sbd_config={"eps": 1e-12, "max_it": 300}) + from mpi4py import MPI + if MPI.COMM_WORLD.Get_rank() != 0: + return # the amplitudes are gathered onto rank 0 only + + assert got.energy == pytest.approx(expected.energy, abs=1e-8) + assert np.array_equal(got.sci_state.ci_strs_a, expected.sci_state.ci_strs_a) + assert np.array_equal(got.sci_state.ci_strs_b, expected.sci_state.ci_strs_b) + + # The same vector up to one global sign -- the comparison that catches a + # per-determinant sign, where every magnitude-based one passes. + a, b = got.sci_state.amplitudes, expected.sci_state.amplitudes + overlap = np.vdot(b, a) + assert abs(overlap) == pytest.approx(1.0, abs=1e-7) + np.testing.assert_allclose(a * np.sign(overlap), b, atol=1e-6) + + # And what users compute from it. + assert got.sci_state.spin_square() == pytest.approx(expected.sci_state.spin_square(), abs=1e-7) + assert _energy_from_rdms(got.sci_state, h1, eri) == pytest.approx(expected.energy, abs=1e-7) + + +# One body, two wrappers: pytest-mpi skips unmarked tests under --only-mpi and marked +# ones without it (see test_sqd_integration.py). +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_standalone(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep) + + +@pytest.mark.mpi +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_mpi(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep) From ca78f4cf138ec147ed97cf455083762e4df2493e Mon Sep 17 00:00:00 2001 From: Jim Garrison Date: Fri, 16 Oct 2026 12:20:13 -0400 Subject: [PATCH 6/6] Give the static TLS block enough headroom for libgomp (#64) 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 --- tox.ini | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/tox.ini b/tox.ini index aaa6118..d149506 100644 --- a/tox.ini +++ b/tox.ini @@ -29,6 +29,25 @@ passenv = OMP_NUM_THREADS PATH LD_LIBRARY_PATH +# Work around Qiskit/qiskit#17213: with qiskit 2.6.0 installed, importing our +# extension fails on Linux with "cannot allocate memory in static TLS block". +# +# glibc sizes the per-thread TLS block at startup and keeps only a small surplus for +# libraries dlopen'ed afterwards that need initial-exec TLS. Every Python extension +# arrives that way, so they all draw on that surplus. qiskit 2.6.0's _accelerate grew +# from 736 to 1544 bytes of it (matrixmultiply 0.3.11 enabling AVX-512 kernels, which +# quadruple one thread-local buffer), which leaves too little for the libgomp our +# extension pulls in. The error names libgomp, but it is whoever asks last that fails. +# +# Raising the surplus restores the headroom. It is reserved per thread and only commits +# the pages actually used, so a generous value costs little. Drop this once 17213 is +# fixed upstream. +# +# Linux-only: this is a glibc mechanism, and macOS links libomp with no equivalent +# limit. Note glibc silently ignores an unrecognized tunable name, so a typo here looks +# like a working config. +setenv = + linux: GLIBC_TUNABLES = glibc.rtld.optional_static_tls=65536 extras = test commands = @@ -52,6 +71,9 @@ allowlist_externals = mpirun setenv = SBD_TEST_NPROCS = {env:SBD_TEST_NPROCS:2} + # setenv is replaced, not merged, in a child environment, so repeat the + # static-TLS headroom [testenv] sets; see the comment there. + linux: GLIBC_TUNABLES = glibc.rtld.optional_static_tls=65536 commands = mpirun -n {env:SBD_TEST_NPROCS} pytest --only-mpi {posargs}