From a887f82f9cd6c5187088c94a2767526d7886f3ad Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 12:58:19 -0400 Subject: [PATCH 01/22] gdb_diag: take a sharded determinant list, and numpy instead of nested lists GDB spans a subspace with an explicit determinant list rather than the Cartesian product TPB uses, which is what makes it the right solver for a sparse subspace. But the binding required every rank to pass the WHOLE list, pinning b_comm_size to 1, and that cascaded into three ceilings: - b_comm is the only dimension that divides the basis. h_comm is a row stride within a block (gdb/mult.h:70) closed by an allreduce (:208) and MakeHelpers ignores it, so at b=1 every rank held the entire basis AND the entire excitation lookup however many ranks were used. Memory did not scale at all. - t_comm_size was forced to 1 too: GDB runs one task per basis-ring station and a single block has one station. - GPU GDB was capped at ONE rank. gdb/mult_thrust.h:310-314 throws unless h_comm_size == 1, and h = ranks/(t*b), so with b=t=1 the helper dimension took every rank. Multi-GPU GDB was unreachable. So b_comm_size now selects the contract: 1 keeps the previous meaning, >1 means each rank passes its own shard. What upstream does not check is checked here, collectively so the verdict is unanimous before any rank raises: shards identical across a b_comm position (the in-memory path never broadcasts the list, unlike the file overload at gdb/sbdiag.h:764), globally sorted and disjoint via the neighbour exchange that load_basis_from_files uses, t <= b, t*b dividing the rank count exactly (upstream's integer division otherwise yields communicators of unequal size), a non-empty shard for the schemes that index config[0] unconditionally, and h == 1 on Thrust. Completeness cannot be checked, so global_dim is returned for the caller to assert on. All six of the app's placement schemes are exposed as one determinant_distribution string. The three legacy booleans are deliberately not: they encode a four-way choice with a priority rule where two are unreachable unless do_redist_alpha_eq is explicitly zeroed. do_shuffle is dropped from GDB_SBD for the same reason h_comm_size was -- upstream parses it and never reads it, so the attribute could only mislead. Determinants now cross as an (n, words) uint64 array, built into det_vector with a single memcpy instead of three copies, and from_strings packs in bulk rather than once per determinant across the boundary. A nested list still works via forcecast. The shape is passed as an explicit std::vector: a braced initializer list is ambiguous under g++ between array_t's ShapeContainer constructor and its copy/move, which nvc++ and clang accept, so the Linux CPU build would otherwise not compile. method 2 and 3 are rejected. They are valid for TPB, where they select Lanczos; GDB has none, and gdb::diag assigns energy only inside its method 0 and 1 branches (gdb/sbdiag.h:418, :501), so passing 2 returned an uninitialized double with no diagonalization having run. Closes half of #22: h_comm_size was already unexposed, and the b_comm_size == 1 restriction it flagged is lifted here using the redistribution APIs the 93ebabec bump brought in. Co-Authored-By: Claude Opus 5 (1M context) --- python/__init__.py | 102 ++++++- python/bindings.cpp | 673 +++++++++++++++++++++++++++++++++++++++----- 2 files changed, 698 insertions(+), 77 deletions(-) diff --git a/python/__init__.py b/python/__init__.py index d74a2fe..75dfd92 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -535,6 +535,28 @@ def from_string(s, bit_length, total_bit_length, device=None): return get_backend(device).from_string(s, bit_length, total_bit_length) +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 + 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. + """ + _ensure_initialized() + return get_backend(device).from_strings(list(strings), bit_length, total_bit_length) + + +def sort_bitarray_array(dets, device=None): + """Sort packed determinants into canonical order, array in and array out. + + The array form of :func:`sort_bitarray`, for a determinant list held as an + ``(n, words)`` array. Deduplicates, like the list form. + """ + _ensure_initialized() + return get_backend(device).sort_bitarray_array(dets) + + def sort_bitarray(dets, device=None): """Sort determinants into canonical order, removing duplicates. @@ -596,7 +618,9 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, def gdb_diag(fcidump, det, sbd_data, - loadname="", savename="", device=None): + loadname="", savename="", device=None, + determinant_distribution="", determinant_grid_a=0, + determinant_grid_b=0): """ Perform GDB diagonalization over an explicit list of determinants. @@ -605,29 +629,79 @@ def gdb_diag(fcidump, det, sbd_data, themselves, so an arbitrary sparse subspace can be diagonalized. Each determinant is a ``2 * norb``-bit configuration packed into words of - ``sbd_data.bit_length`` bits, as returned by :func:`from_string`. Bit ``2 * i`` - is the occupation of spin-alpha orbital ``i`` and bit ``2 * i + 1`` that of - spin-beta orbital ``i``. + ``sbd_data.bit_length`` bits. Bit ``2 * i`` is the occupation of spin-alpha + orbital ``i`` and bit ``2 * i + 1`` that of spin-beta orbital ``i``. + :func:`from_strings` packs a list of bitstrings into the expected + ``(ndets, words)`` array; a nested list is accepted and converted. + + **The shard contract.** ``sbd_data.b_comm_size`` decides what ``det`` means: + + - ``1`` — every rank passes the whole basis. + - ``> 1`` — every rank passes **its own shard**, and the union over b_comm + positions is the basis. This is the only way GDB's memory scales, because + b_comm is the only dimension that divides the basis (and with it the + excitation lookup); the derived helper dimension divides work but not + storage. + + The shard index is ``rank % b_comm_size``, and ranks sharing one must pass + identical determinants — the in-memory path does not broadcast the list. + Shards must be globally sorted and disjoint: shard ``i`` strictly below shard + ``i + 1``. All of this is checked, and a violation raises rather than + silently diagonalizing the wrong subspace. What cannot be checked is + *completeness* — that the union is the basis you meant — so compare the + returned ``global_dim`` against what you expect. + + ``t_comm_size`` must not exceed ``b_comm_size``: GDB runs one task per + basis-ring station and there are exactly ``b_comm_size`` of them. + ``t_comm_size * b_comm_size`` must divide the rank count exactly. On the + Thrust backend the helper dimension must be 1, i.e. every rank goes to + ``t_comm_size * b_comm_size``. Args: fcidump: FCIDump object. - det: Determinants spanning the subspace. Must be distinct; they are - sorted into SBD's canonical order internally, which - :func:`sort_bitarray` reproduces. - sbd_data: GDB_SBD configuration object. ``b_comm_size`` must be 1. - loadname: Path to load initial wavefunction (optional). + det: Determinants spanning the subspace, or this rank's shard of them, as + an ``(ndets, words)`` array. Must be distinct. + sbd_data: GDB_SBD configuration object. ``method`` must be 0 or 1 — GDB + implements Davidson only, so TPB's Lanczos methods 2 and 3 do not + exist here and are rejected. + loadname: Path to load an initial wavefunction (optional). savename: Path prefix to save the final wavefunction to (optional). SBD - writes ``f"{savename}000000.bin"``, holding the determinants in - canonical order and their amplitudes. + 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). + 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 + matvec), ``count``, ``count-sorted``, ``grid-cyclic`` or + ``grid-cyclic-balanced``. Ignored when ``b_comm_size`` is 1. + determinant_grid_a: Grid rows for the grid-cyclic schemes. With + ``determinant_grid_b``, must multiply to ``b_comm_size``; give both or + neither, and the default is the factor pair nearest square. + determinant_grid_b: Grid columns; see ``determinant_grid_a``. Returns: - dict with keys: energy, density, carryover_det, one_p_rdm, two_p_rdm. + dict with keys ``energy``, ``density``, ``carryover_det``, ``one_p_rdm``, + ``two_p_rdm``, ``local_dim``, ``global_dim`` and + ``determinant_distribution``. What is replicated and what is sharded: + + - ``energy`` — replicated, bit-identical on every rank. + - ``density``, ``one_p_rdm``, ``two_p_rdm`` — replicated on every rank. + - ``carryover_det`` — **this rank's shard**, as an ``(n, words)`` array, + not the whole list. For ``carryover_type`` 1 it is split over b_comm and + *duplicated* across the helper dimension, so gathering means taking one + representative per ``rank % b_comm_size``. For types 2 and 3 it is split + over the world communicator with no duplication, so a plain allgather is + correct; those types return the parents together with the new + candidates, i.e. the next subspace rather than only the additions. + - ``local_dim`` / ``global_dim`` — determinants on this rank, and summed + over b_comm. """ _ensure_initialized() backend = get_backend(device) return backend.gdb_diag( - _global_comm, sbd_data, fcidump, det, loadname, savename + _global_comm, sbd_data, fcidump, det, loadname, savename, + determinant_distribution, determinant_grid_a, determinant_grid_b, ) @@ -738,7 +812,9 @@ def print_info(): 'LoadAlphaDets', 'makestring', 'from_string', + 'from_strings', 'sort_bitarray', + 'sort_bitarray_array', 'tpb_diag_from_files', 'tpb_diag', 'gdb_diag', diff --git a/python/bindings.cpp b/python/bindings.cpp index 5aa3c3b..b68f262 100644 --- a/python/bindings.cpp +++ b/python/bindings.cpp @@ -34,8 +34,14 @@ // repo rather than in the vendored upstream submodule. #include +#include +#include +#include +#include + #include #include +#include #include #include #include @@ -134,6 +140,102 @@ static void sbd_pin_offload_device(int mpi_rank) { #define SBD_MODULE_NAME _core #endif +// --------------------------------------------------------------------------- +// GDB distributed-basis helpers +// +// GDB's matvec is systolic over b_comm: each rank owns the bra rows of its basis +// block, and the ket vector plus its index map rotate around the ring +// (gdb/mult.h:39-41 for each rank's initial offset, :196-202 for the hops). Two +// consequences drive everything below. +// +// 1. b_comm is the ONLY dimension that shards memory. h_comm is a row stride +// within a block (mult.h:70) closed by an allreduce (:208), and MakeHelpers +// ignores it, so every h-rank builds the full excitation lookup. +// 2. The ring has exactly mpi_size_b stations and one "task" is one station, so +// t_comm_size <= b_comm_size is structural. +// +// Upstream checks neither, nor that t*b divides the rank count. The failures are +// a segfault and silently ragged communicators respectively, so they are checked +// here instead. +// --------------------------------------------------------------------------- + +/** GDB determinant placement over b_comm. The six schemes of the upstream app. */ +enum class SbdDetDistribution { + input, equal_bra_a, count, count_sorted, grid_cyclic, grid_cyclic_balanced +}; + +/** + * Resolve a scheme name, accepting underscores for dashes as the app does. + * + * Empty means the default. The app reaches its default through three legacy + * booleans with a priority rule (do_redist_alpha_eq, which defaults true, then + * do_sort_det, then do_redist_det -- main.cc:33-56); those are deliberately not + * exposed here, since they encode a four-way choice as three flags where two are + * unreachable unless the first is explicitly zeroed. The name says it directly. + */ +static SbdDetDistribution sbd_resolve_det_distribution(const std::string& raw) { + std::string name = raw; + std::replace(name.begin(), name.end(), '_', '-'); + if (name.empty() || name == "equal-bra-a") return SbdDetDistribution::equal_bra_a; + if (name == "input") return SbdDetDistribution::input; + if (name == "count") return SbdDetDistribution::count; + if (name == "count-sorted") return SbdDetDistribution::count_sorted; + if (name == "grid-cyclic") return SbdDetDistribution::grid_cyclic; + if (name == "grid-cyclic-balanced") return SbdDetDistribution::grid_cyclic_balanced; + throw std::invalid_argument( + "unknown determinant_distribution '" + raw + "'; expected one of: input, " + "equal-bra-a (default), count, count-sorted, grid-cyclic, grid-cyclic-balanced"); +} + +static const char* sbd_det_distribution_name(SbdDetDistribution d) { + switch (d) { + case SbdDetDistribution::input: return "input"; + case SbdDetDistribution::equal_bra_a: return "equal-bra-a"; + case SbdDetDistribution::count: return "count"; + case SbdDetDistribution::count_sorted: return "count-sorted"; + case SbdDetDistribution::grid_cyclic: return "grid-cyclic"; + case SbdDetDistribution::grid_cyclic_balanced: return "grid-cyclic-balanced"; + } + return "unknown"; +} + +/** + * True on every rank iff local_ok holds on every rank of comm. + * + * Every per-rank check below routes through this. Throwing on one rank only + * would leave the others inside the next collective, so the verdict has to be + * unanimous before anyone raises. + */ +static bool sbd_all_ranks_ok(bool local_ok, MPI_Comm comm) { + int local = local_ok ? 1 : 0; + int all = 0; + MPI_Allreduce(&local, &all, 1, MPI_INT, MPI_LAND, comm); + return all != 0; +} + +/** + * Order-insensitive-after-sort 64-bit fingerprint of a shard (FNV-1a). + * + * Used to compare the replicas that ranks sharing a b_comm position must pass: + * the in-memory diag never broadcasts the determinant list (contrast the + * file-based overload's MpiBcast(det,0,h_comm) at gdb/sbdiag.h:764), so a + * divergence there is silently wrong rather than detected. Taken after the local + * sort, so a caller that supplies the same set in a different order still agrees. + */ +static uint64_t sbd_shard_fingerprint(const sbd::det_vector& det) { + uint64_t h = 1469598103934665603ULL; + const auto& flat = det.cflat(); + const unsigned char* bytes = reinterpret_cast(flat.data()); + const size_t n = flat.size() * sizeof(size_t); + for (size_t i = 0; i < n; ++i) { + h ^= static_cast(bytes[i]); + h *= 1099511628211ULL; + } + // Fold in the row count so an empty shard and a zero-filled one differ. + h ^= static_cast(det.size()) + 0x9e3779b97f4a7c15ULL; + return h; +} + PYBIND11_MODULE(SBD_MODULE_NAME, m) { // Set module docstring based on backend #ifdef SBD_THRUST @@ -245,8 +347,13 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { "Task communicator size") .def_readwrite("b_comm_size", &sbd::gdb::SBD::b_comm_size, "Basis communicator size") + // GDB implements Davidson only -- there is no Lanczos anywhere under + // include/sbd/chemistry/gdb/, so TPB's methods 2 and 3 do not exist + // here. gdb_diag rejects them rather than passing them through: see the + // check in gdb_diag for what upstream does with an out-of-range value. .def_readwrite("method", &sbd::gdb::SBD::method, - "Diagonalization method (0=Davidson, 1=Davidson+Ham, 2=Lanczos, 3=Lanczos+Ham)") + "Diagonalization method (0=Davidson, 1=Davidson+Ham). GDB has " + "no Lanczos; TPB's 2 and 3 are not available") .def_readwrite("max_it", &sbd::gdb::SBD::max_it, "Maximum number of iterations") .def_readwrite("max_nb", &sbd::gdb::SBD::max_nb, @@ -259,8 +366,15 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { "Initialization method") .def_readwrite("seed", &sbd::gdb::SBD::seed, "Seed for the initial vector") - .def_readwrite("do_shuffle", &sbd::gdb::SBD::do_shuffle, - "Shuffle determinants flag") + // do_shuffle is deliberately NOT exposed, for the same reason as + // h_comm_size: GDB declares the field (gdb/sbdiag.h:25), parses it from + // argv (:94) and copies it into a local (:219), then never reads it + // again -- the only other "shuffle" under gdb/ is an unrelated + // warp-shuffle comment in mult_thrust.h. TPB does use it + // (tpb/sbdiag.h:888-900), which is why it stays on TPB_SBD. An + // attribute that silently does nothing is worse than an absent one, + // and determinant placement is controlled by + // gdb_diag's determinant_distribution argument instead. .def_readwrite("do_rdm", &sbd::gdb::SBD::do_rdm, "Calculate RDM flag (0=density only, 1=full RDM)") .def_readwrite("carryover_type", &sbd::gdb::SBD::carryover_type, @@ -318,6 +432,85 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { py::arg("bit_length"), py::arg("total_bit_length")); + // Packing one determinant per call from Python is the bottleneck at scale: + // issue #31 measures it at tens of thousands of strings, and a sampled + // subspace is far larger. Loop in C++ and hand back one array instead. + m.def("from_strings", + [](const std::vector& strings, size_t bit_length, + size_t total_bit_length) { + const size_t words = + (total_bit_length + bit_length - 1) / bit_length; + const std::vector shape{ + static_cast(strings.size()), + static_cast(words)}; + py::array_t out(shape); + size_t* data = out.mutable_data(); + std::memset(data, 0, strings.size() * words * sizeof(size_t)); + for (size_t r = 0; r < strings.size(); ++r) { + const std::string& s = strings[r]; + if (s.size() != total_bit_length) { + throw std::invalid_argument( + "from_strings: expected every bitstring to be " + + std::to_string(total_bit_length) + " characters, got " + + std::to_string(s.size()) + " at index " + + std::to_string(r)); + } + size_t* row = data + r * words; + for (size_t i = 0; i < total_bit_length; ++i) { + if (s[total_bit_length - 1 - i] == '1') { + row[i / bit_length] |= + (static_cast(1) << (i % bit_length)); + } + } + } + return out; + }, + "Pack many bitstrings into one (n, words) array, as from_string does " + "for a single determinant", + py::arg("strings"), + py::arg("bit_length"), + py::arg("total_bit_length")); + + // Array in, array out, so a sharded determinant list can be put in canonical + // order without a round trip through Python lists. Deduplicates locally, as + // sbd::sort_bitarray does. + m.def("sort_bitarray_array", + [](py::array_t arr) + -> py::array_t { + if (arr.ndim() != 2) { + // An empty shard is legal and carries no width. + if (arr.ndim() <= 1 && arr.size() == 0) { + const std::vector empty_shape{ + static_cast(0), + static_cast(0)}; + return py::array_t(empty_shape); + } + throw std::invalid_argument( + "sort_bitarray_array expects a 2-D (ndets, words) array"); + } + const size_t n = static_cast(arr.shape(0)); + const size_t words = static_cast(arr.shape(1)); + std::vector> rows(n, std::vector(words)); + const size_t* in = arr.data(); + for (size_t r = 0; r < n; ++r) { + std::memcpy(rows[r].data(), in + r * words, + words * sizeof(size_t)); + } + sbd::sort_bitarray(rows); + const std::vector shape{ + static_cast(rows.size()), + static_cast(words)}; + py::array_t out(shape); + size_t* data = out.mutable_data(); + for (size_t r = 0; r < rows.size(); ++r) { + std::memcpy(data + r * words, rows[r].data(), + words * sizeof(size_t)); + } + return out; + }, + "Sort packed determinants into canonical order, removing duplicates", + py::arg("dets")); + m.def("sort_bitarray", [](std::vector>& dets) { sbd::sort_bitarray(dets); @@ -417,37 +610,385 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { // ======================================================================== // Main GDB diagonalization function (data structure version) // - // The determinant list is the subspace itself: unlike TPB, the subspace is - // not the Cartesian product of two half-determinant lists, so an arbitrary - // sparse set of determinants can be diagonalized. + // The determinant list IS the subspace: unlike TPB, it is not the Cartesian + // product of two half-determinant lists, so an arbitrary sparse set can be + // diagonalized. + // + // THE SHARD CONTRACT. The list may be distributed: + // + // b_comm_size == 1 -> every rank passes the WHOLE basis. + // b_comm_size > 1 -> every rank passes ITS OWN SHARD, and the union over + // b_comm positions is the basis. + // + // Sharding is what makes GDB's memory scale: b_comm is the only dimension + // that divides the basis (and hence the excitation lookup), so at + // b_comm_size == 1 every rank stores everything no matter how many ranks + // there are. Ranks sharing a b_comm position -- i.e. one h_comm, which is + // world_rank % b_comm_size -- must pass IDENTICAL shards, because the + // in-memory gdb::diag never broadcasts the list the way the file-based + // overload does (gdb/sbdiag.h:764). That is checked, not trusted. + // + // Input shards must be globally sorted and disjoint. This is stricter than + // upstream (grid-cyclic tolerates arbitrary shards), but it is what slicing a + // globally sorted list gives you for free, and it lets a caller's partition + // be checked rather than silently producing a wrong subspace. Completeness -- + // that the union is the basis you MEANT -- cannot be checked here, so the + // global dimension is returned for the caller to assert on. // ======================================================================== m.def("gdb_diag", [](py::object py_comm, const sbd::gdb::SBD& sbd_data, const sbd::FCIDump& fcidump, - const std::vector>& det, + py::array_t det_in, const std::string& loadname, - const std::string& savename) { + const std::string& savename, + const std::string& determinant_distribution, + int determinant_grid_a, + int determinant_grid_b) { - // Convert MPI communicator MPI_Comm comm = get_mpi_comm(py_comm); + int mpi_rank; MPI_Comm_rank(comm, &mpi_rank); + int mpi_size; MPI_Comm_size(comm, &mpi_size); - int mpi_rank; - MPI_Comm_rank(comm, &mpi_rank); + const int b_comm_size = sbd_data.b_comm_size; + const int t_comm_size = sbd_data.t_comm_size; - // Every rank passes the whole determinant list, which is what - // sbd::gdb::diag expects only when the basis is not split over - // b_comm. Splitting it would require distributing the determinants - // over b_comm first, as the file-based entry point does. - if (sbd_data.b_comm_size != 1) { + // ---- grid shape ------------------------------------------------- + // These depend only on the config, which every rank passes + // identically, so they can throw directly without a vote. + if (b_comm_size < 1 || t_comm_size < 1) { + throw std::invalid_argument( + "gdb_diag requires b_comm_size >= 1 and t_comm_size >= 1"); + } + // GDB runs one task per basis-ring station and there are exactly + // b_comm_size of them, so a larger t_comm_size starves a rank and + // upstream's MakeHelpers then dereferences an empty lookup + // (gdb/helper.h:736-739 then :761) -- a segfault, not an error. + if (t_comm_size > b_comm_size) { + throw std::invalid_argument( + "gdb_diag requires t_comm_size <= b_comm_size: GDB runs one task " + "per basis-ring station and there are exactly b_comm_size of them, " + "so a larger t_comm_size starves a rank and upstream's MakeHelpers " + "then dereferences an empty lookup (segfault, not an error)"); + } + // gdb::diag dispatches with `if (method == 0) {...} else if (method == 1) + // {...}` and no else, assigning `energy` only inside those branches + // (gdb/sbdiag.h:418, :501). A method of 2 or 3 -- valid for TPB, where + // they select Lanczos -- runs no diagonalization at all and leaves + // `energy` uninitialized. The Thrust build masks the value first + // (`method &= 1`, sbdiag.h:210-211), so it is the CPU and OMP-offload + // backends that would return garbage. Reject it on every backend. + if (sbd_data.method != 0 && sbd_data.method != 1) { throw std::invalid_argument( - "gdb_diag requires b_comm_size == 1; distribute work over " - "t_comm_size, and over the derived helper dimension by " - "changing the rank count (h_comm_size is not settable)"); + "gdb_diag requires method 0 or 1 (Davidson, or Davidson storing " + "the Hamiltonian); GDB has no Lanczos, so TPB's methods 2 and 3 " + "do not exist here and would return an uninitialized energy"); } - if (det.empty()) { - throw std::invalid_argument("gdb_diag requires at least one determinant"); + const long long named_grid = + static_cast(b_comm_size) * static_cast(t_comm_size); + // h_comm_size is derived by integer division upstream with no check, + // so a rank count that is not a multiple of b*t silently yields + // ragged communicators: unequal h_comm sizes in one run, and a rank + // alone in its own b_comm believing the ring has one station. + if (named_grid > mpi_size || mpi_size % named_grid != 0) { + throw std::invalid_argument( + "gdb_diag requires t_comm_size * b_comm_size to divide the rank " + "count exactly (the helper dimension is the quotient); got " + + std::to_string(t_comm_size) + " * " + std::to_string(b_comm_size) + + " against " + std::to_string(mpi_size) + " ranks"); + } + const int h_comm_size = static_cast(mpi_size / named_grid); +#ifdef SBD_THRUST + // gdb/mult_thrust.h:310-314 throws on every kernel launch unless the + // helper dimension is trivial. Say so before any work happens. + if (h_comm_size != 1) { + throw std::invalid_argument( + "GDB on the Thrust backend requires h_comm_size == 1, i.e. " + "t_comm_size * b_comm_size == ranks; got a helper dimension of " + + std::to_string(h_comm_size) + ". Spend every rank on " + "b_comm_size (and t_comm_size <= b_comm_size) instead"); + } +#endif + + // ---- placement scheme ------------------------------------------- + const SbdDetDistribution distribution = + sbd_resolve_det_distribution(determinant_distribution); + const bool uses_grid = + distribution == SbdDetDistribution::grid_cyclic || + distribution == SbdDetDistribution::grid_cyclic_balanced; + int grid_a = determinant_grid_a; + int grid_b = determinant_grid_b; + if (uses_grid) { + if (grid_a < 0 || grid_b < 0) { + throw std::invalid_argument( + "determinant grid dimensions must be positive"); + } + if ((grid_a == 0) != (grid_b == 0)) { + throw std::invalid_argument( + "specify both determinant grid dimensions or neither"); + } + if (grid_a == 0) { + // The factor pair of b_comm_size nearest square, as main.cc does. + grid_a = static_cast(std::sqrt(static_cast(b_comm_size))); + while (grid_a > 1 && b_comm_size % grid_a != 0) --grid_a; + if (grid_a < 1) grid_a = 1; + grid_b = b_comm_size / grid_a; + } + if (static_cast(grid_a) * grid_b != b_comm_size) { + throw std::invalid_argument( + "determinant grid dimensions must multiply to b_comm_size (" + + std::to_string(b_comm_size) + ")"); + } + } else if (grid_a != 0 || grid_b != 0) { + throw std::invalid_argument( + "determinant grid dimensions require a grid-cyclic distribution"); + } + + // ---- orbital count and packing width ---------------------------- + size_t norb = 0; + for (const auto& kv : fcidump.header) { + if (kv.first == std::string("NORB")) { + norb = static_cast(std::atoi(kv.second.c_str())); + } + } + if (norb == 0) { + throw std::invalid_argument( + "gdb_diag could not read NORB from the FCIDUMP header"); + } + const size_t total_bits = 2 * norb; + const size_t bit_length = sbd_data.bit_length; + if (bit_length == 0) { + throw std::invalid_argument("gdb_diag requires bit_length >= 1"); + } + const size_t expected_words = (total_bits + bit_length - 1) / bit_length; + + // ---- the determinant buffer ------------------------------------- + // An empty shard is legal (a rank may own nothing), and numpy gives a + // 1-D zero-length array for an empty list, which carries no width -- + // hence the agreement step below rather than trusting shape(1). + const py::ssize_t ndim = det_in.ndim(); + size_t n_local = 0; + size_t local_words = 0; + bool shape_ok = true; + if (ndim == 2) { + n_local = static_cast(det_in.shape(0)); + local_words = static_cast(det_in.shape(1)); + if (n_local > 0 && local_words == 0) shape_ok = false; + } else if (ndim <= 1 && det_in.size() == 0) { + n_local = 0; + local_words = 0; + } else { + shape_ok = false; + } + if (!sbd_all_ranks_ok(shape_ok, comm)) { + throw std::invalid_argument( + "gdb_diag expects det as a 2-D (ndets, words) array of packed " + "determinants, as from_strings() returns; an empty shard may be " + "an empty array"); + } + + unsigned long long words_local = static_cast(local_words); + unsigned long long words_max = 0; + MPI_Allreduce(&words_local, &words_max, 1, MPI_UNSIGNED_LONG_LONG, + MPI_MAX, comm); + if (words_max == 0) { + throw std::invalid_argument( + "gdb_diag requires at least one determinant somewhere on the " + "communicator"); + } + const size_t words = static_cast(words_max); + if (!sbd_all_ranks_ok(local_words == 0 || local_words == words, comm)) { + throw std::invalid_argument( + "gdb_diag requires every rank's determinants to have the same " + "number of packed words"); + } + if (words != expected_words) { + throw std::invalid_argument( + "gdb_diag got determinants packed into " + std::to_string(words) + + " word(s), but NORB=" + std::to_string(norb) + " with bit_length=" + + std::to_string(bit_length) + " implies " + + std::to_string(expected_words) + + "; pack with the same bit_length the config carries"); + } + + // det_vector's row width is a process-global property fixed by the + // first container built, so set it while the GIL is held and report a + // mismatch before any work is done. The half-determinant width is set + // too, matching the upstream app: the grid-cyclic path builds + // half-determinant containers before any diagonalization runs. + try { + sbd::det_vector::init_elem_size(words); + sbd::det_vector::init_elem_size( + (norb + bit_length - 1) / bit_length); + } catch (const std::length_error&) { + throw std::invalid_argument( + "gdb_diag was already called in this process with a different " + "number of words per determinant, which SBD fixes for the " + "lifetime of the process. Keep norb and bit_length fixed, or " + "run the new problem in a fresh process."); + } + + sbd::det_vector det; + det.resize(n_local); + if (n_local > 0) { + std::memcpy(det.flat().data(), det_in.data(), + n_local * words * sizeof(size_t)); + } + + // SBD indexes the subspace with binary searches, so each shard must be + // in canonical order; an unsorted list silently yields a wrong energy. + // sort_bitarray also removes duplicates, which would leave part of the + // subspace unreachable, so reject a shrink rather than dropping them. + const size_t n_before = det.size(); + sbd::sort_bitarray(det); + if (!sbd_all_ranks_ok(det.size() == n_before, comm)) { + throw std::invalid_argument( + "gdb_diag requires distinct determinants"); + } + + // ---- communicators ---------------------------------------------- + // Built here, before diag builds its own set, because the placement + // schemes redistribute over b_comm. With both named dimensions + // trivial, h_comm is the whole communicator, so skip the split. + const bool split_comms = (b_comm_size > 1 || t_comm_size > 1); + MPI_Comm h_comm = comm, b_comm = MPI_COMM_NULL, t_comm = MPI_COMM_NULL; + if (split_comms) { + sbd::gdb::DetBasisCommunicator(comm, h_comm_size, b_comm_size, + t_comm_size, h_comm, b_comm, t_comm); + } + struct CommGuard { + bool owns; MPI_Comm *h, *b, *t; + ~CommGuard() { + if (!owns) return; + if (*h != MPI_COMM_NULL) MPI_Comm_free(h); + if (*b != MPI_COMM_NULL) MPI_Comm_free(b); + if (*t != MPI_COMM_NULL) MPI_Comm_free(t); + } + } comm_guard{split_comms, &h_comm, &b_comm, &t_comm}; + + // ---- shard agreement across h_comm ------------------------------ + { + int h_size = 1; + MPI_Comm_size(h_comm, &h_size); + if (h_size > 1) { + unsigned long long fp = sbd_shard_fingerprint(det); + unsigned long long lo = 0, hi = 0; + MPI_Allreduce(&fp, &lo, 1, MPI_UNSIGNED_LONG_LONG, MPI_MIN, h_comm); + MPI_Allreduce(&fp, &hi, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, h_comm); + if (!sbd_all_ranks_ok(lo == hi, comm)) { + throw std::invalid_argument( + "gdb_diag requires ranks sharing a b_comm position to pass " + "identical determinants: the shard index is " + "world_rank % b_comm_size, and the in-memory path does not " + "broadcast the list, so a divergence would silently " + "diagonalize different subspaces"); + } + } + } + + // ---- shards are globally sorted and disjoint --------------------- + if (b_comm_size > 1) { + int rank_b = 0, size_b = 0; + MPI_Comm_rank(b_comm, &rank_b); + MPI_Comm_size(b_comm, &size_b); + std::vector my_last(words, 0), prev_last(words, 0); + const bool have_mine = det.size() > 0; + if (have_mine) { + std::memcpy(my_last.data(), &det.cflat()[(det.size() - 1) * words], + words * sizeof(size_t)); + } + int send_n = have_mine ? 1 : 0; + int recv_n = 0; + const int dst = (rank_b + 1 < size_b) ? rank_b + 1 : MPI_PROC_NULL; + const int srcr = (rank_b > 0) ? rank_b - 1 : MPI_PROC_NULL; + MPI_Sendrecv(&send_n, 1, MPI_INT, dst, 91, + &recv_n, 1, MPI_INT, srcr, 91, b_comm, MPI_STATUS_IGNORE); + MPI_Sendrecv(my_last.data(), static_cast(words), + MPI_UNSIGNED_LONG, dst, 92, + prev_last.data(), static_cast(words), + MPI_UNSIGNED_LONG, srcr, 92, b_comm, MPI_STATUS_IGNORE); + // Only comparable when both sides hold something; an empty shard + // in between is skipped rather than treated as a failure. + bool ordered = true; + if (recv_n == 1 && have_mine) { + std::vector my_first(words, 0); + std::memcpy(my_first.data(), det.cflat().data(), + words * sizeof(size_t)); + const bool duplicate = + std::memcmp(my_first.data(), prev_last.data(), + words * sizeof(size_t)) == 0; + ordered = !duplicate && !sbd::less_from_back(my_first, prev_last); + } + if (!sbd_all_ranks_ok(ordered, comm)) { + throw std::invalid_argument( + "gdb_diag requires the input shards to be globally sorted and " + "disjoint: shard i must hold a strictly lower range than shard " + "i+1, where the shard index is world_rank % b_comm_size. Slice " + "a globally sorted determinant list to get this"); + } + } + + // ---- schemes that cannot take an empty shard -------------------- + // redistribution/reordering index config[0] unconditionally + // (framework/bit_manipulation.h:566, :577) and SaveWavefunction does + // basis[0].size() (caop/basic/restart.h:35), all UB when a rank owns + // nothing. redistribution_bitarray and the grid-cyclic path are safe. + const bool needs_nonempty = + distribution == SbdDetDistribution::count || + distribution == SbdDetDistribution::count_sorted || + !savename.empty(); + if (needs_nonempty && !sbd_all_ranks_ok(det.size() > 0, comm)) { + throw std::invalid_argument( + std::string("gdb_diag with determinant_distribution='") + + sbd_det_distribution_name(distribution) + + "'" + (savename.empty() ? "" : " or a savename") + + " requires a non-empty shard on every rank: upstream indexes the " + "first determinant unconditionally on those paths"); + } + + // ---- redistribute over b_comm ----------------------------------- + // Nothing to place with a single basis block, and equal-bra-a there + // would still pay an allgather of every alpha string, so skip it. + if (b_comm_size > 1 && distribution != SbdDetDistribution::input) { + py::gil_scoped_release release; + switch (distribution) { + case SbdDetDistribution::equal_bra_a: + sbd::redistribution_equal_bra_a(det, bit_length, total_bits, + b_comm); + break; + case SbdDetDistribution::count: + sbd::redistribution(det, bit_length, total_bits, b_comm); + break; + case SbdDetDistribution::count_sorted: + sbd::redistribution(det, bit_length, total_bits, b_comm); + sbd::reordering(det, bit_length, total_bits, b_comm); + break; + case SbdDetDistribution::grid_cyclic: + case SbdDetDistribution::grid_cyclic_balanced: + sbd::gdb::redistribution_grid_bra_ab_cyclic( + det, bit_length, total_bits, + static_cast(grid_a), static_cast(grid_b), + b_comm); + if (distribution == SbdDetDistribution::grid_cyclic_balanced) { + sbd::redistribution_bitarray(det, b_comm); + } + // Mandatory: the grid-cyclic output is locally unsorted, + // and MakeHelpers requires each shard in canonical order. + sbd::sort_bitarray(det); + break; + case SbdDetDistribution::input: + break; + } + } + + // ---- global dimension ------------------------------------------- + // Summed over b_comm, since h_comm and t_comm hold replicas. + unsigned long long local_dim = static_cast(det.size()); + unsigned long long global_dim = local_dim; + if (b_comm_size > 1) { + MPI_Allreduce(&local_dim, &global_dim, 1, MPI_UNSIGNED_LONG_LONG, + MPI_SUM, b_comm); } #ifdef SBD_THRUST @@ -468,76 +1009,80 @@ PYBIND11_MODULE(SBD_MODULE_NAME, m) { sbd_pin_offload_device(mpi_rank); #endif - // Output variables - double energy; + double energy = 0.0; std::vector density; sbd::det_vector co_det; std::vector> one_p_rdm; std::vector> two_p_rdm; - // det_vector packs each determinant into a fixed number of words, - // which is a process-global property of the type: it is fixed by the - // first determinant container built in the process and cannot be - // changed afterwards. Set it here, while the GIL is still held, so - // that a mismatch is reported before any work is done. - try { - sbd::det_vector::init_elem_size(det[0].size()); - } catch (const std::length_error&) { - throw std::invalid_argument( - "gdb_diag was already called in this process with a different " - "number of words per determinant, which SBD fixes for the " - "lifetime of the process. Keep norb and bit_length fixed, or " - "run the new problem in a fresh process."); - } - - // SBD indexes the subspace with binary searches, so the determinants - // must be in canonical order; an unsorted list silently yields a - // wrong energy. Sort here rather than asking the caller to; the - // exposed sort_bitarray reproduces this order for a caller that - // needs it. sort_bitarray also removes duplicates, which would leave - // part of the subspace unreachable, so reject them instead of - // dropping them silently. - std::vector> det_sorted(det); - sbd::sort_bitarray(det_sorted); - if (det_sorted.size() != det.size()) { - throw std::invalid_argument( - "gdb_diag requires distinct determinants"); - } - - std::vector> co_det_vvs; - { - // Release GIL for long computation py::gil_scoped_release release; - sbd::det_vector det_packed(det_sorted.begin(), det_sorted.end()); - - sbd::gdb::diag(comm, sbd_data, fcidump, det_packed, + sbd::gdb::diag(comm, sbd_data, fcidump, det, loadname, savename, energy, density, co_det, one_p_rdm, two_p_rdm); - for (const auto& row : co_det) { - co_det_vvs.emplace_back(row.begin(), row.end()); + // With do_rdm == 0 the occupation density is computed only where + // mpi_rank_t == 0 (gdb/sbdiag.h:525) and left untouched + // elsewhere, so ranks off that plane would see an empty list. + // Fill them in so the Python contract does not depend on t. + if (t_comm_size > 1 && t_comm != MPI_COMM_NULL) { + unsigned long long n_den = density.size(), n_den_max = 0; + MPI_Allreduce(&n_den, &n_den_max, 1, MPI_UNSIGNED_LONG_LONG, + MPI_MAX, t_comm); + if (n_den_max > 0) { + density.resize(static_cast(n_den_max), 0.0); + MPI_Bcast(density.data(), static_cast(n_den_max), + MPI_DOUBLE, 0, t_comm); + } } } - // Return results as dictionary + // Carryover comes back as a shard, not a gathered list: for + // carryover_type 1 it is split over b_comm AND duplicated across + // h_comm (gdb/carryover.h), and for 2/3 it is split over the world + // communicator with no duplication (gdb/expansion.h:821-822). Handing + // back the raw shard keeps the memory win; the caller gathers if it + // wants the whole list, taking one representative per b_comm position + // for type 1. + const size_t n_co = co_det.size(); + // Explicit shape vector: a braced list here is ambiguous under g++ + // (array_t's ShapeContainer constructor versus its copy/move), which + // nvc++ and clang both accepted -- so the Linux CPU build broke while + // the macOS and Thrust builds passed. + const std::vector co_shape{ + static_cast(n_co), + static_cast(words)}; + py::array_t co_out(co_shape); + if (n_co > 0) { + std::memcpy(co_out.mutable_data(), co_det.cflat().data(), + n_co * words * sizeof(size_t)); + } + py::dict results; results["energy"] = energy; results["density"] = density; - results["carryover_det"] = co_det_vvs; + results["carryover_det"] = co_out; results["one_p_rdm"] = one_p_rdm; results["two_p_rdm"] = two_p_rdm; + results["local_dim"] = local_dim; + results["global_dim"] = global_dim; + results["determinant_distribution"] = + std::string(sbd_det_distribution_name(distribution)); return results; }, - "Perform GDB diagonalization over an explicit list of determinants", + "Perform GDB diagonalization over an explicit list of determinants, which " + "may be sharded across b_comm", py::arg("comm"), py::arg("sbd_data"), py::arg("fcidump"), py::arg("det"), py::arg("loadname") = "", - py::arg("savename") = ""); + py::arg("savename") = "", + py::arg("determinant_distribution") = "", + py::arg("determinant_grid_a") = 0, + py::arg("determinant_grid_b") = 0); // ======================================================================== // Main TPB diagonalization function (file-based version) From 60868aa1fa4e8d94ff23db361338a55c8c5370f6 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 12:58:19 -0400 Subject: [PATCH 02/22] test: cover GDB, which had no test anywhere GDB was untested in this wrapper and, at the pinned upstream commit, upstream's only GDB test exercises the distribution helpers rather than a diagonalization. The anchor needs no pinned number: TPB's subspace is one GDB can also express, so interleaving every (alpha, beta) pair into a full determinant points both solvers at exactly the same Hilbert space and their energies must agree. On h2o they agree to 1.4e-14, and the full 275^2 interleave reproduces the -76.23594663 published for that alpha list, which test_reference_energies.py asserts through TPB. Also covered: the sharded path against a TPB run of the same subspace; the task dimension becoming usable once b_comm_size allows it; all six placement schemes agreeing, since placement is a load-balancing decision and must not move the answer; and each guard with its own case -- non-disjoint shards, a rank count the grid cannot tile, unknown and misconfigured placement schemes, duplicate determinants, and Lanczos methods. from_strings is checked elementwise against a per-determinant from_string loop, and a nested list is checked to still work. The Fe4S4 case runs in a subprocess: det_vector's row width is fixed per process, and h2o packs into one 64-bit word where Fe4S4 needs two, so they cannot share one. That is the constraint gdb_diag reports, and forking is the workaround it prescribes. Co-Authored-By: Claude Opus 5 (1M context) --- test/test_gdb_equivalence.py | 560 +++++++++++++++++++++++++++++++++++ 1 file changed, 560 insertions(+) create mode 100644 test/test_gdb_equivalence.py diff --git a/test/test_gdb_equivalence.py b/test/test_gdb_equivalence.py new file mode 100644 index 0000000..d270346 --- /dev/null +++ b/test/test_gdb_equivalence.py @@ -0,0 +1,560 @@ +# 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. + +"""Check GDB against TPB on a subspace both can express. + +GDB spans the subspace with the determinants it is handed; TPB spans it with the +Cartesian product of an alpha and a beta determinant list. The product is a +subspace GDB can also express -- interleave every (alpha, beta) pair into one +full determinant -- so the two solvers can be pointed at exactly the same +Hilbert space and their energies must agree. That equality is the anchor here, +and it needs no pinned reference value. + +It also gives GDB a ground-truth number for free. The energies published with +the upstream TPB data are for the full product of the alpha list with itself, so +the full interleave of ``h2o-1em3-alpha.txt`` must reproduce the published +-76.23594663 -- see ``test_reference_energies.py``, which asserts the same value +through TPB. + +Before this file, GDB had no test anywhere: not in the wrapper, and not in +``vendor/sbd-upstream/tests/functionality`` apart from upstream's +``test_gdb_grid_distribution.cc``, which exercises the distribution helpers +rather than a diagonalization. +""" + +from __future__ import annotations + +import itertools +import pathlib +import subprocess +import sys + +import pytest + +import sbd + +# One word per determinant for h2o either way: 24 orbitals is 24 bits of alpha +# (TPB's half determinants) and 48 bits of alpha+beta (GDB's full ones), both +# under 64. ``det_vector::init_elem_size`` fixes the word count process-wide on +# first use, so every case in this file -- and in test_reference_energies.py, +# which shares the process -- must agree on it. +BIT_LENGTH = 64 + +# Enough alpha strings to make a non-trivial subspace, few enough that the +# product stays small: the |A|^2 scaling is the reason this is 24 and not 275. +ALPHA_LIMIT = 24 + +# Published for the full 275 x 275 product of h2o-1em3-alpha.txt. +H2O_1EM3_ENERGY = -76.23594663 + + +def _interleave(alpha: str, beta: str) -> str: + """Interleave two norb-bit strings into one 2*norb-bit GDB determinant. + + Bit ``2 * i`` is alpha orbital ``i`` and bit ``2 * i + 1`` is beta orbital + ``i``, counting bits from the right -- the order the GDB app's README + specifies ("the rightmost bit corresponds to alpha-spin orbital 1, the next + to beta-spin orbital 1") and the order ``from_string`` then packs. + """ + a_rev, b_rev = alpha[::-1], beta[::-1] + return "".join(a_rev[i] + b_rev[i] for i in range(len(a_rev)))[::-1] + + +def _alpha_strings(path, limit=None): + with open(path, encoding="utf-8") as handle: + strings = [line.strip() for line in handle if line.strip()] + return strings[:limit] if limit else strings + + +def _tpb_energy(backend, fcidump, alpha_strings, norb, **overrides): + """Diagonalize the product of ``alpha_strings`` with itself, via TPB.""" + config = backend.TPB_SBD() + config.eps = 1e-10 + config.max_it = 200 + config.bit_length = BIT_LENGTH + for name, value in overrides.items(): + setattr(config, name, value) + half = backend.sort_bitarray( + [backend.from_string(s, BIT_LENGTH, norb) for s in alpha_strings] + ) + return sbd.tpb_diag(fcidump, half, half, config)["energy"] + + +# Placement is a property of the call, not of the solver config, so these are +# kwargs on gdb_diag rather than fields to setattr onto GDB_SBD. +_CALL_KWARGS = ("determinant_distribution", "determinant_grid_a", "determinant_grid_b") + + +def _gdb_result(backend, fcidump, det, norb, **overrides): + """Diagonalize an explicit determinant list (or this rank's shard) via GDB. + + ``det`` may be a list of bitstrings, which is packed here, or an already + packed ``(n, words)`` array. + """ + config = backend.GDB_SBD() + config.eps = 1e-10 + config.max_it = 200 + config.bit_length = BIT_LENGTH + call = {k: overrides.pop(k) for k in list(overrides) if k in _CALL_KWARGS} + for name, value in overrides.items(): + setattr(config, name, value) + # Bitstrings need packing; anything else is already packed, as an array or as + # the nested sequence the binding's forcecast accepts. + if len(det) and isinstance(det[0], str): + det = sbd.from_strings(list(det), BIT_LENGTH, 2 * norb) + return sbd.gdb_diag(fcidump, det, config, **call) + + +def _gdb_energy(backend, fcidump, det, norb, **overrides): + """The energy alone, for the many cases that assert only on it.""" + return _gdb_result(backend, fcidump, det, norb, **overrides)["energy"] + + +def _shard(full, b_comm_size, rank): + """This rank's slice of a globally sorted list. + + Mirrors SBD's own ``q = N/p``, remainder-to-the-low-ranks split + (``balanced_begin``, framework/bit_manipulation.h:1523-1531), and indexes by + ``rank % b_comm_size`` because that is the b_comm position. + """ + if b_comm_size == 1: + return full + index = rank % b_comm_size + quotient, remainder = divmod(len(full), b_comm_size) + begin = index * quotient + min(index, remainder) + end = begin + quotient + (1 if index < remainder else 0) + return full[begin:end] + + +@pytest.fixture(scope="module") +def h2o(data_dir, backend): + """The h2o FCIDUMP, its orbital count, the alpha-list path, and a truncation.""" + molecule_dir = data_dir / "h2o" + fcidump = backend.LoadFCIDump(str(molecule_dir / "fcidump.txt")) + norb = int(fcidump.header["NORB"]) + alpha_path = molecule_dir / "h2o-1em3-alpha.txt" + return fcidump, norb, alpha_path, _alpha_strings(alpha_path, ALPHA_LIMIT) + + +def test_gdb_matches_tpb_on_the_same_subspace(backend, h2o): + """Both solvers give the same energy for the same Hilbert space. + + The strong form of the check: no reference value, so it cannot pass by + coincidence with a wrong Hamiltonian or a misread FCIDUMP. Only the + determinant *representation* differs -- |A| x |A| half determinants against + |A|^2 interleaved full ones. + """ + fcidump, norb, _, alpha = h2o + product = [_interleave(a, b) for a, b in itertools.product(alpha, alpha)] + assert len(product) == len(alpha) ** 2 + + tpb = _tpb_energy(backend, fcidump, alpha, norb) + gdb = _gdb_energy(backend, fcidump, product, norb) + assert gdb == pytest.approx(tpb, abs=1e-9) + + +def test_gdb_rejects_duplicate_determinants(backend, h2o): + """A repeated determinant is an error, not a silently smaller subspace. + + ``sort_bitarray`` deduplicates, so the binding compares sizes afterwards and + raises. Recorded here because TPB takes the opposite branch and drops + duplicates silently -- the divergence tracked in issue #31. If that is ever + reconciled in GDB's favour, this test is what has to change. + """ + fcidump, norb, _, alpha = h2o + doubled = [_interleave(alpha[0], alpha[0])] * 2 + with pytest.raises(ValueError, match="distinct determinants"): + _gdb_energy(backend, fcidump, doubled, norb) + + +def test_gdb_rejects_a_basis_split_the_ranks_cannot_tile(backend, h2o): + """b_comm_size > 1 needs ranks to match it. + + Splitting the basis is supported now, but the grid still has to tile: upstream + derives the helper dimension as ``ranks / (t * b)`` by integer division and + never checks the remainder, so asking for more basis blocks than there are + ranks would silently produce communicators of unequal size. Serially that + means any b_comm_size above 1 is refused. + """ + fcidump, norb, _, alpha = h2o + product = [_interleave(alpha[0], b) for b in alpha] + with pytest.raises(ValueError, match="divide the rank count"): + _gdb_energy(backend, fcidump, product, norb, b_comm_size=2) + + +def test_gdb_rejects_lanczos_methods(backend, h2o): + """Methods 2 and 3 are refused rather than returning an uninitialized energy. + + They are valid for TPB, where they select Lanczos. GDB has no Lanczos: + ``gdb::diag`` handles only ``method == 0`` and ``method == 1`` and assigns + ``energy`` only inside those branches, so passing 2 through would return + whatever was on the stack -- with no diagonalization having run. + """ + fcidump, norb, _, alpha = h2o + product = [_interleave(alpha[0], b) for b in alpha] + with pytest.raises(ValueError, match="method 0 or 1"): + _gdb_energy(backend, fcidump, product, norb, method=2) + + +@pytest.mark.slow +def test_gdb_reproduces_the_published_h2o_energy(backend, h2o): + """The full interleave reproduces the energy published for this alpha list. + + The subspace is 275^2 = 75,625 determinants -- the same Hilbert space + ``test_reference_energies.py`` covers through TPB in its fast tier, but GDB + carries its dimension explicitly rather than as a product, so the determinant + list itself is 75,625 entries. Marked slow for the cost of that rather than + for being intractable: measured at 5.5 s of diagonalization plus 0.2 s to + build the list, on 8 host threads, returning -76.2359466308 against the + published -76.23594663. Promote it to the fast tier if that budget is fine. + """ + fcidump, norb, alpha_path, _ = h2o + # The full list, not the fixture's ALPHA_LIMIT truncation: the published + # energy is for the product of all of it with itself. + alpha = _alpha_strings(alpha_path) + product = [_interleave(a, b) for a, b in itertools.product(alpha, alpha)] + assert len(product) == len(alpha) ** 2 + + energy = _gdb_energy(backend, fcidump, product, norb) + # Quoted to eight decimals upstream, so compared to that rather than to the + # solver's own 1e-10 convergence tolerance. + assert energy == pytest.approx(H2O_1EM3_ENERGY, abs=1e-8) + + +@pytest.mark.mpi +def test_gdb_under_mpi_matches_tpb_on_the_same_subspace(backend, h2o): + """Distributing GDB over ranks does not change the answer. + + GDB decomposes as ``t_comm_size x b_comm_size x helper``. Both named + dimensions are pinned to 1 on the in-memory path -- b_comm_size because every + rank passes the whole determinant list, t_comm_size because GDB enumerates one + task per b_comm rank and so cannot split a single task further (see + ``test_gdb_rejects_a_task_split``). Ranks therefore all land on the derived + helper dimension, and this checks that spending them there is harmless. + + The comparison is against a TPB run of the same subspace in the same process, + which reaches the same Hilbert space through an independent decomposition + (``adet_comm_size``). No pinned reference value is needed, and a broken helper + distribution shows up as disagreement. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + fcidump, norb, _, alpha = h2o + product = [_interleave(a, b) for a, b in itertools.product(alpha, alpha)] + + gdb = _gdb_energy(backend, fcidump, product, norb) + tpb = _tpb_energy(backend, fcidump, alpha, norb, adet_comm_size=size) + + if comm.Get_rank() == 0: + assert gdb == pytest.approx(tpb, abs=1e-9), ( + f"{size} ranks (helper={size}) gave {gdb}, but TPB on the same " + f"subspace gave {tpb}" + ) + + +def test_gdb_rejects_a_task_split(backend, h2o): + """t_comm_size > 1 is refused rather than segfaulting. + + ``t_comm_size <= b_comm_size`` is structural, not incidental. GDB's matvec + rotates the ket around ``b_comm`` as a ring (gdb/mult.h:39-41, :196-202), so + the ring has exactly ``b_comm_size`` stations and one "task" is one station; + ``t_comm`` parallelizes stations, and cannot have more workers than there are + stations. ``MakeHelpers`` encodes this as ``task_end = mpi_size_b`` split + across ``t_comm`` (gdb/helper.h:736-739). + + Upstream does not check it: a starved rank resizes ``exidx`` to zero and then + reads ``exidx[0].slide`` anyway because its ``task_begin`` is nonzero + (helper.h:761), which is a null dereference during helper construction -- + observed as EXC_BAD_ACCESS at 0x0 on 2 and 4 ranks before this guard existed. + With ``b_comm_size`` pinned to 1 in memory, that makes any ``t_comm_size > 1`` + fatal. + + The check runs serially because the guard is a configuration check that + precedes any communication. + """ + fcidump, norb, _, alpha = h2o + product = [_interleave(alpha[0], b) for b in alpha] + with pytest.raises(ValueError, match="t_comm_size <= b_comm_size"): + _gdb_energy(backend, fcidump, product, norb, t_comm_size=2) + + +# --------------------------------------------------------------------------- +# Distributed basis: b_comm_size > 1, where each rank passes only its shard +# --------------------------------------------------------------------------- + + +@pytest.fixture(scope="module") +def h2o_product(backend, h2o): + """The truncated product basis, packed and globally sorted.""" + fcidump, norb, _, alpha = h2o + strings = [_interleave(a, b) for a, b in itertools.product(alpha, alpha)] + packed = sbd.sort_bitarray_array( + sbd.from_strings(strings, BIT_LENGTH, 2 * norb) + ) + assert packed.shape[0] == len(alpha) ** 2 + return packed + + +def test_from_strings_matches_the_per_determinant_form(backend, h2o): + """The bulk packer agrees elementwise with from_string called in a loop. + + from_strings exists because the per-call form costs one boundary crossing per + determinant, which issue #31 measures as a real cost at tens of thousands of + strings. It is only worth having if it is the same function, hence this. + """ + fcidump, norb, _, alpha = h2o + strings = [_interleave(a, b) for a, b in itertools.product(alpha[:6], alpha[:6])] + bulk = sbd.from_strings(strings, BIT_LENGTH, 2 * norb) + one_by_one = [backend.from_string(s, BIT_LENGTH, 2 * norb) for s in strings] + assert bulk.tolist() == one_by_one + + +def test_gdb_accepts_a_nested_list(backend, h2o, h2o_product): + """A list of lists still works, so existing callers are unaffected. + + The binding takes a forcecast array, which converts a nested sequence, so the + numpy contract is an addition rather than a break. + """ + fcidump, norb, _, _ = h2o + as_list = [list(row) for row in h2o_product] + assert _gdb_energy(backend, fcidump, as_list, norb) == pytest.approx( + _gdb_energy(backend, fcidump, h2o_product, norb), abs=1e-12 + ) + + +def test_gdb_reports_the_dimensions_it_diagonalized(backend, h2o, h2o_product): + """global_dim is what the caller needs to check completeness itself. + + Whether the union of shards is the basis the caller *meant* cannot be checked + inside the binding, so the dimension is reported back instead. + """ + fcidump, norb, _, _ = h2o + result = _gdb_result(backend, fcidump, h2o_product, norb) + assert result["global_dim"] == h2o_product.shape[0] + assert result["local_dim"] == h2o_product.shape[0] + assert result["determinant_distribution"] == "equal-bra-a" + + +@pytest.mark.mpi +def test_gdb_sharded_basis_matches_the_whole_basis(backend, h2o, h2o_product): + """Splitting the basis across b_comm does not change the energy. + + The core claim of the distributed path: each rank passes only its slice, and + the answer still matches TPB on the same subspace. Checked against a TPB run + in this same process, which reaches the same Hilbert space through an + independent decomposition, so no pinned value is involved. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + fcidump, norb, _, alpha = h2o + + gdb = _gdb_energy(backend, fcidump, _shard(h2o_product, size, comm.Get_rank()), + norb, b_comm_size=size) + tpb = _tpb_energy(backend, fcidump, alpha, norb, adet_comm_size=size) + if comm.Get_rank() == 0: + assert gdb == pytest.approx(tpb, abs=1e-9) + + +@pytest.mark.mpi +def test_gdb_task_dimension_works_once_the_ring_has_stations(backend, h2o, h2o_product): + """t_comm_size > 1 becomes usable exactly when b_comm_size allows it. + + ``t <= b`` is structural -- one task per basis-ring station -- so with the + basis in one block the task dimension is unreachable. Shard the basis and it + opens up. Needs at least 4 ranks for t=2, b=2. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + if size < 4 or size % 4: + pytest.skip(f"needs a rank count divisible by 4 for t=2 x b=2, got {size}") + fcidump, norb, _, alpha = h2o + + gdb = _gdb_energy(backend, fcidump, _shard(h2o_product, 2, comm.Get_rank()), + norb, b_comm_size=2, t_comm_size=2) + tpb = _tpb_energy(backend, fcidump, alpha, norb, adet_comm_size=size) + if comm.Get_rank() == 0: + assert gdb == pytest.approx(tpb, abs=1e-9) + + +@pytest.mark.mpi +@pytest.mark.parametrize( + "scheme", + ["input", "equal-bra-a", "count", "count-sorted", + "grid-cyclic", "grid-cyclic-balanced"], +) +def test_gdb_placement_does_not_change_the_energy(backend, h2o, h2o_product, scheme): + """All six placement schemes agree. + + Placement is a load-balancing decision -- which rank owns which determinants, + and how the ring is laid out -- so it must not move the answer. ``input`` keeps + the caller's slices, ``equal-bra-a`` equalizes distinct alpha strings (what + actually balances the matvec, whose outer loop is over alpha), ``count`` and + ``count-sorted`` equalize determinant counts, and the two grid-cyclic schemes + spread alpha and beta keys over a 2-D rank grid. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + fcidump, norb, _, alpha = h2o + + gdb = _gdb_energy(backend, fcidump, _shard(h2o_product, size, comm.Get_rank()), + norb, b_comm_size=size, determinant_distribution=scheme) + tpb = _tpb_energy(backend, fcidump, alpha, norb, adet_comm_size=size) + if comm.Get_rank() == 0: + assert gdb == pytest.approx(tpb, abs=1e-9) + + +@pytest.mark.mpi +def test_gdb_rejects_shards_that_are_not_disjoint(backend, h2o, h2o_product): + """Overlapping shards are refused, not silently diagonalized. + + Every rank passing the whole list while claiming b_comm_size > 1 is the + mistake this catches: each rank would believe its copy was a shard, and the + ring would rotate duplicated blocks. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + if size < 2: + pytest.skip("needs at least 2 ranks to have distinct shards") + fcidump, norb, _, _ = h2o + with pytest.raises(ValueError, match="globally sorted and disjoint"): + _gdb_energy(backend, fcidump, h2o_product, norb, b_comm_size=size) + + +@pytest.mark.mpi +def test_gdb_rejects_a_rank_count_the_grid_does_not_tile(backend, h2o, h2o_product): + """t * b must divide the rank count exactly. + + Upstream derives the helper dimension by integer division and never checks the + remainder, which silently yields communicators of unequal size and a rank alone + in its own basis ring. Refuse it instead. + """ + from mpi4py import MPI + + comm = MPI.COMM_WORLD + size = comm.Get_size() + fcidump, norb, _, _ = h2o + bad = size + 1 + with pytest.raises(ValueError, match="divide the rank count"): + _gdb_energy(backend, fcidump, _shard(h2o_product, bad, comm.Get_rank()), + norb, b_comm_size=bad) + + +def test_gdb_rejects_an_unknown_placement_scheme(backend, h2o, h2o_product): + """A misspelled scheme names the alternatives rather than falling back.""" + fcidump, norb, _, _ = h2o + with pytest.raises(ValueError, match="unknown determinant_distribution"): + _gdb_energy(backend, fcidump, h2o_product, norb, + determinant_distribution="equal-bra-alpha") + + +def test_gdb_accepts_underscores_in_the_scheme_name(backend, h2o, h2o_product): + """``grid_cyclic`` and ``grid-cyclic`` are the same scheme, as upstream has it.""" + fcidump, norb, _, _ = h2o + result = _gdb_result(backend, fcidump, h2o_product, norb, + determinant_distribution="grid_cyclic") + assert result["determinant_distribution"] == "grid-cyclic" + + +@pytest.mark.parametrize( + "grid_a,grid_b,scheme,message", + [ + (2, 0, "grid-cyclic", "both determinant grid dimensions or neither"), + (3, 3, "grid-cyclic", "multiply to b_comm_size"), + (1, 1, "equal-bra-a", "require a grid-cyclic distribution"), + ], +) +def test_gdb_validates_the_determinant_grid(backend, h2o, h2o_product, + grid_a, grid_b, scheme, message): + """The grid dimensions are checked against b_comm_size before any work. + + ``redistribution_grid_bra_ab_cyclic`` validates the product itself, but against + the communicator it is handed, so checking here names b_comm_size in the error. + """ + fcidump, norb, _, _ = h2o + with pytest.raises(ValueError, match=message): + _gdb_energy(backend, fcidump, h2o_product, norb, + determinant_distribution=scheme, + determinant_grid_a=grid_a, determinant_grid_b=grid_b) + + +# Upstream's only general-determinant data, and the only case here at a realistic +# size. It lives in the app directory rather than under data/, which is alpha-only. +FE4S4_DIR = ( + "vendor/sbd-upstream/apps/chemistry_gdb_selected_basis_diagonalization" +) +# Measured by this wrapper, NOT published upstream: run.sh ships no expected value +# and data/ has reference tables only for the TPB molecules. So this is a +# regression guard against our own verified result, not an independent check. It +# was identical across b_comm_size 1, 2 and 4 and across t_comm_size 1 and 2. +FE4S4_ENERGY = -326.6982518821 + + +@pytest.mark.slow +def test_gdb_fe4s4_at_realistic_size(): + """59,536 determinants over 36 orbitals, upstream's own GDB input. + + Everything else here is h2o at a few hundred to tens of thousands of + determinants; this is the only case with a real determinant list and a real + orbital count, and the only one whose four input files are a ready-made + four-way partition (equal sizes, disjoint, individually sorted, concatenation + globally sorted). Roughly 21 s on 8 host threads. + + Runs in a **subprocess**, because ``det_vector``'s row width is fixed for the + lifetime of a process: h2o packs into one 64-bit word (48 bits) and Fe4S4 needs + two (72 bits), so the two cannot share a process. That is the constraint + ``gdb_diag`` reports rather than a problem with this case, and forking is the + workaround it tells callers to use. + """ + root = pathlib.Path(__file__).resolve().parents[1] / FE4S4_DIR + if not root.is_dir(): + pytest.skip(f"upstream GDB app data not found at {root}") + + script = f""" +import sbd +root = {str(root)!r} +backend = sbd.get_backend() +fcidump = backend.LoadFCIDump(root + "/fcidump_Fe4S4.txt") +norb = int(fcidump.header["NORB"]) +strings = [] +for index in range(4): + with open(root + "/det%d.txt" % index) as handle: + strings.extend(line.strip() for line in handle if line.strip()) +assert norb == 36, norb +assert len(strings) == 59536, len(strings) +config = backend.GDB_SBD() +config.eps, config.max_it, config.bit_length = 1e-6, 30, {BIT_LENGTH} +det = sbd.from_strings(strings, {BIT_LENGTH}, 2 * norb) +result = sbd.gdb_diag(fcidump, det, config) +print("ENERGY", repr(result["energy"]), "DIM", result["global_dim"]) +""" + completed = subprocess.run( + [sys.executable, "-c", script], + capture_output=True, text=True, timeout=1800, + ) + assert completed.returncode == 0, ( + f"subprocess failed:\n{completed.stdout[-2000:]}\n{completed.stderr[-2000:]}" + ) + line = [l for l in completed.stdout.splitlines() if l.startswith("ENERGY")] + assert line, f"no energy reported:\n{completed.stdout[-2000:]}" + _, energy, _, dim = line[-1].split() + assert int(dim) == 59536 + assert float(energy) == pytest.approx(FE4S4_ENERGY, abs=1e-8) From dbca48dced13ef37606e7ddddfd95fa84950ac83 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 12:58:46 -0400 Subject: [PATCH 03/22] examples/gdb: drivers for GDB diagonalization and heatbath expansion Two drivers, both staying on the in-memory entry point -- determinant text is read in Python and handed to the binding as an array, never through SBD's file-based overload. run_gdb_diag.py runs a single diagonalization over an explicit determinant list. It takes full-determinant files (defaulting to upstream's four Fe4S4 files) or forms the product basis from an alpha list, which is the subspace TPB would build and so the cross-check against it. With --b_comm_size it shards: when the number of files is any multiple of b, each rank reads its own contiguous block of files, which has to be contiguous rather than strided because shard i must hold a strictly lower range than shard i+1. run_gdb_heatbath.py iterates: diagonalize, let SBD expand the subspace from the resulting wavefunction, diagonalize the larger subspace, repeat. carryover_type 2 and 3 return the parents together with the new candidates, so one round's result is the next round's subspace and the loop needs nothing extra. It is structured as a cutoff LADDER because the expansion reaches a self-consistent size for a fixed heatbath_cutoff and then stops growing -- rounds run at one cutoff until the energy settles or the subspace stops growing, then the next rung starts, with --max_dim capping the whole run. --seed picks the starting subspace (determinant files, the Hartree-Fock determinant alone, an interleaved alpha list, or an arbitrary bitstring file) so one driver produces every row of a seed comparison. The README documents what each path does and the constraints that are easy to get wrong: b_comm_size is the only dimension that divides memory; t <= b because there is one task per basis-ring station; t*b must divide the rank count; the helper dimension must be 1 on Thrust, which is why multi-GPU GDB needs b > 1; expansion and carryover run on the host even in a GPU build, so OMP_NUM_THREADS should stay generous; and heatbath_truncation prunes parents rather than admitting candidates, so leaving it at 0 is almost always right. It also records that upstream's Fe4S4 "GDB" data is the full 244x244 product of AlphaDets.txt rather than a sparse subspace, which is worth knowing before reading it as a GDB benchmark, and points at apps/gen_dets for producing pre-split shard files rather than shipping a second tool for the same job. Deliberately no timing or hardware figures: the point is which paths exist and what each costs in memory and constraints, and users should measure on their own systems rather than inherit ours. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 3 + examples/gdb/README.md | 322 +++++++++++++++++++++ examples/gdb/run_gdb_diag.py | 452 +++++++++++++++++++++++++++++ examples/gdb/run_gdb_heatbath.py | 476 +++++++++++++++++++++++++++++++ examples/tpb/README.md | 2 + 5 files changed, 1255 insertions(+) create mode 100644 examples/gdb/README.md create mode 100755 examples/gdb/run_gdb_diag.py create mode 100755 examples/gdb/run_gdb_heatbath.py diff --git a/README.md b/README.md index e091a44..5f0ad5b 100644 --- a/README.md +++ b/README.md @@ -58,6 +58,9 @@ matter for it. - [`examples/tpb/`](examples/tpb/README.md) — tensor-product basis: standalone TPB diagonalization, the SQD loops, and the subspace-enlargement driver. +- [`examples/gdb/`](examples/gdb/README.md) — general determinant basis: standalone + GDB diagonalization over an explicit determinant list, and an iterative + heatbath-expansion driver. ## Integration with qiskit-addon-sqd diff --git a/examples/gdb/README.md b/examples/gdb/README.md new file mode 100644 index 0000000..7a22a20 --- /dev/null +++ b/examples/gdb/README.md @@ -0,0 +1,322 @@ +# GDB examples — general determinant basis + +GDB spans the subspace with the determinants it is given, rather than with the +Cartesian product of an alpha and a beta list that TPB uses. TPB's dimension is +`|adet| x |bdet|`; GDB's is exactly the number of determinants passed. That is what +makes it the right solver for an arbitrary sparse subspace — a set of sampled +bitstrings used *as sampled*, with no product completion. + +For backend selection, `--device` values and bundled test data, see +[`../README.md`](../README.md). For TPB and the SQD loops, see +[`../tpb/README.md`](../tpb/README.md). + +## run_gdb_diag.py — standalone GDB diagonalization + +Stays on the in-memory entry point (`sbd.gdb_diag`) throughout: determinant text is +read in Python and handed to the binding as a list. SBD's own file-based entry +point is deliberately not used. + +```bash +# Fe4S4, upstream's own GDB data: 4 files, 59,536 determinants, 36 orbitals. +# Upstream publishes no reference energy for this case. +python run_gdb_diag.py + +# Shard the basis: Fe4S4's four files, one per rank. For this data that is +# already the balanced globally-sorted split, so no redistribution is needed. +mpirun -np 4 -x OMP_NUM_THREADS=8 python run_gdb_diag.py --b_comm_size 4 + +# Spend ranks on both named dimensions (t <= b, and t*b must divide the ranks) +mpirun -np 8 python run_gdb_diag.py --b_comm_size 4 --t_comm_size 2 + +# Choose how determinants are placed across b_comm +mpirun -np 4 python run_gdb_diag.py --b_comm_size 4 \ + --determinant_distribution grid-cyclic + +# A subspace TPB can also express: interleave an alpha list with itself into the +# full |A|^2 product basis. This is the cross-check against tpb_diag -- same +# subspace, two independent solvers, energies must agree. +python run_gdb_diag.py \ + --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \ + --from-alpha ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt \ + --alpha-limit 30 + +# Grow the subspace with SBD's own heatbath expansion (one round; the expanded +# list comes back as carryover_det) +python run_gdb_diag.py --carryover_type 2 --heatbath_cutoff 1e-4 +``` + +The determinant bit order is the one the upstream app documents: reading from the +right, alpha orbital 1, beta orbital 1, alpha orbital 2, and so on — i.e. bit +`2*i` is alpha orbital `i` and bit `2*i + 1` is beta orbital `i`. + +## run_gdb_heatbath.py — grow the subspace, round after round + +Selected CI (HCI) driven from Python: diagonalize, let SBD expand the subspace from +the resulting wavefunction, diagonalize the larger subspace, repeat. The loop needs +no extra machinery — `carryover_type` 2 and 3 return the parents *together with* the +new candidates, so one round's result **is** the next round's subspace. + +```bash +# Fe4S4 from upstream's shipped subspace, one cutoff +python run_gdb_heatbath.py --cutoffs 1e-3 + +# A ladder, stopping before the subspace passes 2M determinants +python run_gdb_heatbath.py --cutoffs 1e-3,1e-4,1e-5 --max_dim 2000000 + +# The no-input null: start from the Hartree-Fock determinant alone +python run_gdb_heatbath.py --seed hf --cutoffs 1e-3,1e-4 + +# Sharded. At --b_comm_size == ranks with t=1 the expansion comes back already +# sharded for the next round, so nothing has to be gathered. +mpirun -np 4 python run_gdb_heatbath.py --b_comm_size 4 --cutoffs 1e-3 +mpirun -np 4 python run_gdb_heatbath.py --b_comm_size 4 --device gpu --cutoffs 1e-4 + +# Record the (dimension, energy) series for a comparison table +python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log ladder.json +``` + +### Why a cutoff *ladder* + +`heatbath_cutoff` admits a candidate when its estimated contribution exceeds the +threshold, so **for a fixed cutoff the subspace reaches a self-consistent size and +stops growing** — more rounds then buy nothing. Going deeper means lowering the +cutoff. So `--cutoffs` takes a ladder: rounds run at one cutoff until the energy +stops moving (`--energy_tol`, the standard HCI criterion) or the subspace stops +growing, then the next rung begins. `--max_dim` caps the whole run, which is what +keeps rows of a comparison cost-matched rather than cutoff-matched. + +On Fe4S4 (36 orbitals, 54 electrons) starting from upstream's four shipped files, +a `1e-3,1e-4` ladder walks the subspace out like this: + +| rung | dimension | energy | +|---|---|---| +| seed (upstream's four files) | 59,536 | −326.6982518821 | +| `1e-3` round 1 | 63,569 | −326.7298282871 | +| `1e-3` settled | 63,887 | −326.7302438694 | +| `1e-4` round 1 | 568,538 | −326.7773323983 | +| `1e-4` round 3 | 755,107 | −326.7839718082 | + +`1e-3` gains ~32 mHa for 7% more determinants and then saturates — that is the fixed +point for that cutoff, not convergence of the method. One step to `1e-4` multiplies the +subspace by roughly nine. Use `--max_dim` to stop before a rung outgrows your memory. + +The energy is variational, so it must fall monotonically; the driver flags a rise, +which would mean the subspace shrank or a round failed to converge. For h2o and n2 the +FCI limits in that basis (−76.24377680 and −109.04874199) are hard ceilings a correct +run can never cross, which makes a sparse run self-checking even without a reference. + +### Seeds + +`--seed` selects where the starting subspace comes from, so one driver produces +every row of a seed comparison: + +| `--seed` | starting subspace | +|---|---| +| `files` *(default)* | determinant files, defaulting to upstream's four Fe4S4 files | +| `hf` | the Hartree-Fock determinant alone, built from the FCIDUMP header — the no-input null | +| `from-alpha` | an alpha list interleaved with itself, i.e. a TPB-shaped product space | +| `strings` | an arbitrary bitstring file, e.g. sampled configurations | + +Which seed is best is system-dependent, and worth measuring rather than assuming. +The `hf` seed is the useful baseline: it starts from a single determinant and needs no +input subspace at all, so it shows what the classical expansion achieves on its own. +On Fe4S4 it climbs from −326.5243550095 (the HF energy) into the same range the +file-seeded run reaches, at a comparable dimension — so on that system the expansion +is doing most of the work and a pre-existing subspace adds little. Whether that holds +for your system is exactly the kind of thing to check with `--log` and a dimension +sweep. + +Compare seeds **at matched dimension**, not at matched cutoff. Different seeds reach +different sizes from the same cutoff, so a cutoff-matched comparison mostly reports +subspace size rather than seed quality. `--max_dim` and the `--log` series are there +for that. + +### `--heatbath_truncation` is not the cutoff + +It discards **parents** by weight *before* expansion starts, and its default of 0 +(keep every parent) is almost always what you want. Setting it to `1e-4` on the +59,536-determinant Fe4S4 wavefunction cut the subspace to 506 rather than growing it. + +## Expected results + +| case | energy | note | +|---|---|---| +| h2o, first 24 alpha interleaved into a 576-determinant product basis | **-76.0588897208** | matches `tpb_diag` on the same subspace to 1.4e-14 | +| h2o, full 275² = 75,625-determinant interleave | **-76.2359466308** | against the **-76.23594663** published for that alpha list | +| Fe4S4, 59,536 determinants, 36 orbitals | **-326.6982518821** | measured here; upstream publishes no reference for this case | + +The h2o product-basis energy is unchanged across `b_comm_size` 1, 2 and 4, across +`t_comm_size` 1 and 2, and across all six placement schemes. Fe4S4 is likewise +identical at `b=1`, `b=4` (one file per rank) and `b=2, t=2`. + +A note on what sharding buys: `b_comm_size` divides **memory**, not necessarily time. +At modest sizes the ring communication can offset the extra parallelism, so the reason +to shard is to hold a subspace that would not fit on one rank — and, on GPUs, because +it is the only way to use more than one card at all (see [On GPUs](#on-gpus)). Measure +on your own hardware rather than extrapolating from anyone else's. + +## Preparing pre-split determinant files (optional) + +Neither driver needs pre-split files: both read a list and slice it deterministically, +so one big file is always correct. Splitting is an **I/O and memory optimization** — +with one file every rank reads the whole thing and keeps only its slice, whereas with +one file per rank the read is parallel and peak memory is a single shard. Both drivers +take that path automatically when the number of `--detfiles` equals `--b_comm_size`. + +Upstream already ships the tool for producing them: **`apps/gen_dets`** (built as +`gdet`). It takes an alpha determinant list, forms the full alpha × beta product, +sorts it, and writes it as N shards: + +```bash +cd vendor/sbd-upstream/apps/gen_dets +# edit Configuration for your compiler, then: +make +mpirun -np 1 ./gdet --adetfile AlphaDets.txt \ + --detfiles det0.txt,det1.txt,det2.txt,det3.txt \ + --bit_length 20 --norb 36 +``` + +The number of output files comes from `--detfiles`, **not** from the rank count — rank +0 writes all of them itself, splitting with `get_mpi_range`, so `-np 1` is enough +(upstream's own `run.sh` passes `-np 4`, which works but is not required). Set +`--norb` to the orbital count and give one file per rank you intend to run on, e.g. +eight files for eight GPUs. + +The result satisfies what `gdb_diag` requires of input shards — each shard sorted, +shards disjoint, and the concatenation globally sorted so shard *i* sits strictly +below shard *i+1* — because `get_mpi_range` is the same `q = N/p` +remainder-to-the-low-ranks split the drivers use. + +Two gaps worth knowing. `gdet` only takes an **alpha list** and can only produce the +alpha × beta *product*, so there is no upstream way to shard an arbitrary sparse or +heatbath-expanded determinant list into files; feed those to the drivers as a single +file (or in memory, which is what `run_gdb_heatbath.py` does between rounds). And +upstream's error messages point at a `scripts/sort-basis-shards.py` that is not +shipped — there is no `scripts/` directory in the tree. + +### Note on Fe4S4's shipped basis + +Worth knowing when reading benchmark numbers: `AlphaDets.txt` next to `gen_dets` +holds exactly **244** alpha determinants, 244² = **59,536**, and the four det files +total 59,536 and are byte-identical to the GDB app's. So upstream's Fe4S4 "GDB" +benchmark subspace is the **full Cartesian product** of 244 alpha determinants with +itself — a TPB-shaped space written out as an explicit determinant list, not a sparse +one. TPB fed `AlphaDets.txt` returns the same −326.6982518821 that GDB returns on the +four files. + +## MPI decomposition + +GDB decomposes as `t_comm_size × b_comm_size × helper`, with its own field names +rather than TPB's `adet`/`bdet`/`task`. Only two are settable; `helper` is derived: + +``` +helper = ranks / (t_comm_size × b_comm_size) +``` + +**`t_comm_size × b_comm_size` must divide the rank count exactly.** Upstream takes +that quotient by integer division and never checks the remainder, which silently +produces communicators of unequal size and a rank alone in its own basis ring, so +`gdb_diag` refuses it. + +### What each dimension buys + +| | parallelizes | shards memory? | constraint | +|---|---|---|---| +| **`b_comm_size`** | the basis itself — each rank owns a block, and the blocks form the ring the ket rotates around | **yes — the only one that does** | ≥ 1 | +| **`t_comm_size`** | the ring's stations, each task rank starting at its own offset | no — costs `t` replicas of the ket | **`t ≤ b`** | +| **`helper`** | rows within a block (`idet % helper`), reduced by an allreduce | no — every helper rank builds the full excitation lookup | none on CPU; **must be 1 on Thrust** | + +`t ≤ b` is structural, not a limitation of this wrapper: the matvec rotates the ket +around `b_comm` as a ring (`gdb/mult.h:39-41`, `:196-202`), so there are exactly +`b_comm_size` stations and one task is one station. Upstream does not check it and a +starved task rank dereferences an empty lookup (`gdb/helper.h:761`), so `gdb_diag` +rejects it rather than letting it segfault. + +Because only `b` shards memory, `--b_comm_size 1` on many ranks means every rank +holds the whole basis *and* the whole excitation lookup no matter how many ranks +there are. That is the wall to watch at scale. + +### The shard contract + +With `--b_comm_size R > 1` each rank passes **its own shard**, not the whole list: + +- the shard index is `rank % R`, and ranks sharing one must pass **identical** + determinants (the in-memory path does not broadcast the list); +- shards must be **globally sorted and disjoint** — shard `i` strictly below shard + `i + 1`, which slicing a globally sorted list gives you for free; +- the union must be the basis you meant. That part cannot be checked here, so + compare the returned `global_dim` against what you expect. + +Everything else is checked and raises rather than quietly diagonalizing the wrong +subspace. + +### Placement across `b_comm` + +`--determinant_distribution` chooses who owns what. All six agree on the energy — +placement is a load-balancing decision — so pick by balance, not by result. + +| scheme | equalizes | leaves uneven | +|---|---|---| +| `input` | nothing; your slices are kept | whatever imbalance you passed | +| `equal-bra-a` *(default)* | distinct **alpha** strings per rank, which is what the matvec's outer loop runs over | determinant counts, when one alpha has many betas | +| `count` | determinant counts | one alpha's determinants can straddle ranks | +| `count-sorted` | as `count`, plus a reorder within each rank | same | +| `grid-cyclic` | spreads alpha and beta keys cyclically over a `grid_a × grid_b` grid | no count rebalance after | +| `grid-cyclic-balanced` | `grid-cyclic`, then a count rebalance | — | + +The grid schemes exist because a one-dimensional alpha split can be badly skewed for +an irregular subspace — exactly the sampled-bitstring case. `--determinant_grid_a` +and `--determinant_grid_b` must multiply to `b_comm_size`; give both or neither, and +the default is the factor pair nearest square. + +For Fe4S4 the four shipped files are already the balanced globally-sorted split +(14,884 each, disjoint, individually sorted, concatenation in order), so at +`--b_comm_size 4` the driver hands file *i* to shard *i* and `input` placement is +already count-balanced. + +### On GPUs + +GDB has **Thrust kernels only** (`--device gpu`, NVIDIA). There is no `omp target` +code under `include/sbd/chemistry/gdb/` at all, so under `--device gpu-omp` a GDB run +pins a device and then diagonalizes on the host; the driver warns when it resolves +there. AMD has no Thrust build, so AMD GDB is CPU-only. + +The Thrust path additionally requires **`helper == 1`** (`gdb/mult_thrust.h:310-314`, +checked on every kernel launch), which means every rank must go to +`t_comm_size × b_comm_size`. Since `helper = ranks / (t × b)`, leaving the basis in one +block caps GPU GDB at a single rank — so **multi-GPU GDB requires `--b_comm_size`**. + +Because `h = ranks / (t · b)`, leaving the basis in a single block caps GPU GDB at one +rank — the helper dimension takes every rank you add. So **multi-GPU GDB requires +`--b_comm_size`**, e.g. `-np 4 --b_comm_size 4`. Asking for more ranks than `t · b` is +refused up front, naming the helper dimension, rather than throwing from inside a +kernel launch. + +Only the Davidson runs on the device. `gdb/expansion.h` and `gdb/carryover.h` contain +no `thrust::` code and are included outside any `SBD_THRUST` guard (`inc_all.h:23`), so +heatbath expansion and carryover selection run on the host with OpenMP even in a GPU +build. Keep `OMP_NUM_THREADS` generous on GPU runs — one thread per rank throttles that +half of every iteration — and expect less than full GPU utilisation during a ladder for +the same reason. + +### `method` + +Only 0 and 1. TPB's 2 and 3 select Lanczos, which GDB does not implement; `gdb_diag` +rejects them rather than letting `gdb::diag` fall through both of its branches and +return an uninitialized energy. + +### Outputs, and what is replicated + +| key | distribution | +|---|---| +| `energy` | replicated, bit-identical on every rank | +| `density`, `one_p_rdm`, `two_p_rdm` | replicated on every rank | +| `carryover_det` | **this rank's shard**, as an `(n, words)` array. For `carryover_type` 1 it is split over `b_comm` *and duplicated* across the helper dimension, so gathering means one representative per `rank % b_comm_size`; for types 2 and 3 it is split over the world communicator with no duplication, so a plain allgather is correct | +| `local_dim` / `global_dim` | determinants on this rank, and summed over `b_comm` | +| `savename` | one file per shard, `f"{savename}{rank_b:06d}.bin"` — rank 0's file is not the whole wavefunction, and GDB has no combined matrix-form dump | + +## See Also + +- [`../README.md`](../README.md) — backend selection, test data, performance tips +- [`../tpb/README.md`](../tpb/README.md) — TPB and the SQD loops +- [Repository README](../../README.md) — installation, API reference diff --git a/examples/gdb/run_gdb_diag.py b/examples/gdb/run_gdb_diag.py new file mode 100755 index 0000000..eaa83af --- /dev/null +++ b/examples/gdb/run_gdb_diag.py @@ -0,0 +1,452 @@ +#!/usr/bin/env python3 + +# 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. + +""" +Standalone GDB diagonalization over an explicit determinant list, in memory. + +GDB (general determinant basis) spans the subspace with the determinants it is +given. TPB spans it with the Cartesian product of an alpha and a beta list, so +TPB's dimension is |adet| x |bdet| while GDB's is exactly the number of +determinants passed. That is the whole point of GDB: an arbitrary sparse +subspace -- for example a set of sampled bitstrings used as sampled, with no +product completion. + +This driver stays on the in-memory entry point (``sbd.gdb_diag``) throughout. +Determinant text is read in Python and handed to the binding as a list; SBD's +own file-based entry point is deliberately not used. + +The determinant list can be sharded: with ``--b_comm_size R`` each rank passes only +its own slice, which is the only way GDB's memory scales, because b_comm is the one +dimension that divides the basis (and with it the excitation lookup). The derived +helper dimension divides work but not storage. ``--t_comm_size`` must not exceed +``--b_comm_size`` (one task per basis-ring station), and their product must divide +the rank count. + +Usage: + # Fe4S4, upstream's own GDB data: 4 files, 59,536 determinants, 36 orbitals. + # Note that upstream publishes no reference energy for this case. + python run_gdb_diag.py + + # Shard Fe4S4's four files one per rank: file i goes to b_comm position i, + # which for this data is already the balanced globally-sorted split + mpirun -np 4 -x OMP_NUM_THREADS=8 python run_gdb_diag.py --b_comm_size 4 + + # Spend ranks on both named dimensions (t <= b, and t*b must divide the ranks) + mpirun -np 8 python run_gdb_diag.py --b_comm_size 4 --t_comm_size 2 + + # Choose how determinants are placed across b_comm + mpirun -np 4 python run_gdb_diag.py --b_comm_size 4 \ + --determinant_distribution grid-cyclic + + # A subspace TPB can also express: interleave an alpha list with itself to + # form the full |A|^2 product basis. This is the cross-check against + # tpb_diag -- same subspace, two independent solvers, energies must agree. + python run_gdb_diag.py \ + --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \ + --from-alpha ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt \ + --alpha-limit 30 + + # Grow the subspace with SBD's own heatbath expansion (one round; the + # expanded list comes back as carryover_det) + python run_gdb_diag.py --carryover_type 2 --heatbath_cutoff 1e-4 +""" + +import argparse +import sys +import time + +import numpy as np + +# Fe4S4 is the only general-determinant data upstream ships, and it lives in the +# app directory rather than under data/ (which is alpha-only, i.e. TPB). +_UPSTREAM = "../../vendor/sbd-upstream" +_GDB_APP = f"{_UPSTREAM}/apps/chemistry_gdb_selected_basis_diagonalization" +_DEFAULT_FCIDUMP = f"{_GDB_APP}/fcidump_Fe4S4.txt" +_DEFAULT_DETFILES = ",".join(f"{_GDB_APP}/det{i}.txt" for i in range(4)) + + +def parse_args(): + """Parse command line arguments for all GDB_SBD parameters.""" + parser = argparse.ArgumentParser( + description="GDB diagonalization over an explicit determinant list", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + + parser.add_argument('--device', default='cpu', + choices=['cpu', 'gpu', 'gpu-omp', 'auto'], + help="Compute backend. 'gpu' is NVHPC Thrust, which is the " + "only GPU backend with GDB kernels; 'gpu-omp' has none, " + "so GDB runs on the host there") + + # --- input ------------------------------------------------------------ + parser.add_argument('--fcidump', default=_DEFAULT_FCIDUMP, + help='FCIDUMP file defining the Hamiltonian') + parser.add_argument('--detfiles', default=_DEFAULT_DETFILES, + help='Comma-separated files of full determinants, one ' + '2*norb-bit string per line. Concatenated in Python ' + 'and passed as a single in-memory list') + parser.add_argument('--from-alpha', default='', metavar='FILE', + dest='from_alpha', + help='Instead of --detfiles, read a norb-bit alpha list and ' + 'form the full |A|^2 product basis by interleaving it ' + 'with itself. This is the subspace TPB would build from ' + 'the same file, so the two solvers are comparable') + parser.add_argument('--alpha-limit', type=int, default=0, dest='alpha_limit', + help='With --from-alpha, keep only the first N alpha strings. ' + '0 means all, which costs |A|^2 determinants') + + # --- MPI decomposition ------------------------------------------------ + parser.add_argument('--b_comm_size', type=int, default=1, + help='Basis communicator size, i.e. how many shards the ' + 'determinant list is split into. The only dimension that ' + 'divides memory') + parser.add_argument('--t_comm_size', type=int, default=1, + help='Task communicator size. Must not exceed b_comm_size: GDB ' + 'runs one task per basis-ring station and there are ' + 'exactly b_comm_size of them') + parser.add_argument('--determinant_distribution', default='', + choices=['', 'input', 'equal-bra-a', 'count', 'count-sorted', + 'grid-cyclic', 'grid-cyclic-balanced'], + help="How to place determinants across b_comm. Default " + "equal-bra-a, which equalizes distinct alpha strings per " + "rank -- the quantity that balances the matvec") + parser.add_argument('--determinant_grid_a', type=int, default=0, + help='Grid rows for the grid-cyclic schemes; with ' + '--determinant_grid_b must multiply to b_comm_size') + parser.add_argument('--determinant_grid_b', type=int, default=0, + help='Grid columns; see --determinant_grid_a') + + # --- diagonalization -------------------------------------------------- + # GDB implements Davidson only: there is no Lanczos anywhere under + # include/sbd/chemistry/gdb/. gdb::diag is `if(method==0){} else if(method==1){}` + # with no else, and `energy` is assigned only inside those branches + # (gdb/sbdiag.h:418, :501), so method 2 or 3 returns an uninitialized double on + # the CPU backend. The Thrust path is accidentally safe -- sbdiag.h:210-211 does + # `method &= 1`. Hence choices=[0, 1], unlike run_sbd_diag.py's [0, 1, 2, 3]. + parser.add_argument('--method', type=int, default=0, choices=[0, 1], + help='0=Davidson, 1=Davidson storing the Hamiltonian. GDB has ' + 'no Lanczos, so TPB methods 2 and 3 do not exist here') + parser.add_argument('--iteration', '--max_it', type=int, default=100, + dest='max_it', help='Maximum Davidson iterations') + parser.add_argument('--block', '--max_nb', type=int, default=10, + dest='max_nb', help='Maximum number of basis vectors') + parser.add_argument('--tolerance', '--eps', type=float, default=1e-6, + dest='eps', help='Convergence tolerance') + parser.add_argument('--max_time', type=float, default=1e10, + help='Maximum wall time in seconds') + parser.add_argument('--init', type=int, default=0, + help='Initial vector policy') + parser.add_argument('--seed', type=int, default=1729, + help='Seed for a random initial vector (init != 0)') + parser.add_argument('--bit_length', type=int, default=64, + help='Bits per packed word. Words per determinant is ' + 'ceil(2*norb / bit_length)') + + # --- RDMs ------------------------------------------------------------- + parser.add_argument('--rdm_output', default='', metavar='FILE', + help='Save spin-summed rdm1 and rdm2 to one .npz file. ' + 'Implies do_rdm=1') + + # --- carryover / expansion ------------------------------------------- + parser.add_argument('--carryover_type', type=int, default=0, + choices=[0, 1, 2, 3], + help='0=off, 1=keep the top --ratio fraction by weight, ' + '2/3=weight-truncate then heatbath-expand (2 and 3 are ' + 'the two heatbath variants). 2/3 GROW the subspace') + parser.add_argument('--carryover_ratio', '--ratio', type=float, default=0.0, + dest='ratio', + help='Fraction of determinants kept by carryover_type 1') + parser.add_argument('--carryover_threshold', '--threshold', type=float, + default=0.01, dest='threshold', + help='Weight threshold for carryover_type 1') + parser.add_argument('--heatbath_cutoff', type=float, default=1e-4, + help='Integral-magnitude cutoff admitting a candidate ' + 'determinant (carryover_type 2/3)') + parser.add_argument('--heatbath_truncation', type=float, default=0.0, + help='Weight threshold applied to parents before expanding') + parser.add_argument('--heatbath_batch_size', type=int, default=1000000, + help='Heatbath expansion batch size per rank') + + # --- wavefunction I/O ------------------------------------------------- + parser.add_argument('--loadname', default='', + help='Load an initial wavefunction from this file') + parser.add_argument('--savename', default='', + help='Save the final wavefunction under this prefix') + + return parser.parse_args() + + +def read_det_strings(paths): + """Read full-determinant bitstrings from one or more text files. + + Returns them in file order, unfiltered: gdb_diag sorts its own input and + rejects duplicates, so both properties are left for it to enforce rather + than silently repaired here. + """ + strings = [] + for path in paths: + with open(path, encoding='utf-8') as handle: + for line in handle: + line = line.strip() + if line: + strings.append(line) + return strings + + +def interleave_spin_strings(alpha, beta): + """Interleave an alpha and a beta bitstring into one GDB determinant string. + + Both inputs are norb-bit strings written most-significant-first, as in + SBD's alpha determinant files. The result is 2*norb bits in the order the + GDB app documents: reading from the right, alpha orbital 1, beta orbital 1, + alpha orbital 2, and so on -- i.e. bit 2*i is alpha orbital i and bit + 2*i + 1 is beta orbital i, which is what from_string() then packs. + """ + a_rev = alpha[::-1] + b_rev = beta[::-1] + return ''.join(a_rev[i] + b_rev[i] for i in range(len(a_rev)))[::-1] + + +def shard_bounds(total, b_comm_size, index): + """This rank's slice of a globally sorted list of ``total`` determinants. + + Mirrors SBD's own ``q = N/p``, remainder-to-the-low-ranks split + (``balanced_begin``, framework/bit_manipulation.h:1523-1531), so an ``input`` + placement is already the count-balanced one. + """ + quotient, remainder = divmod(total, b_comm_size) + begin = index * quotient + min(index, remainder) + return begin, begin + quotient + (1 if index < remainder else 0) + + +def product_basis(alpha_strings): + """Every (alpha, beta) pair from one list: the subspace TPB spans from it.""" + return [interleave_spin_strings(a, b) + for a in alpha_strings for b in alpha_strings] + + +def main(): + args = parse_args() + + import sbd + from sbd.sbd_solver import assemble_rdms + + sbd.init(device=args.device) + + rank = sbd.get_rank() + size = sbd.get_world_size() + device = sbd.get_device() + + if rank == 0: + print("=" * 70) + print("SBD GDB - general determinant basis, in memory") + print("=" * 70) + sbd.print_info() + print() + + # Every rank builds the same list: that is what the in-memory entry point + # expects, and it is why b_comm_size must be 1. + fcidump = sbd.LoadFCIDump(args.fcidump) + norb = int(fcidump.header["NORB"]) + nelec = int(fcidump.header["NELEC"]) + total_bits = 2 * norb + + b_size = args.b_comm_size + shard_index = rank % b_size if b_size > 1 else 0 + + if args.from_alpha: + alpha = read_det_strings([args.from_alpha]) + if args.alpha_limit: + alpha = alpha[:args.alpha_limit] + bad = [s for s in alpha if len(s) != norb] + if bad: + print(f"ERROR: --from-alpha strings must be {norb} bits " + f"(NORB from the FCIDUMP); found one of length {len(bad[0])}", + file=sys.stderr) + return 1 + det_strings = product_basis(alpha) + source = (f"{args.from_alpha} ({len(alpha)} alpha strings -> " + f"{len(alpha)}^2 product determinants)") + whole = sbd.sort_bitarray_array( + sbd.from_strings(det_strings, args.bit_length, total_bits)) + begin, end = shard_bounds(whole.shape[0], b_size, shard_index) + det = whole[begin:end] if b_size > 1 else whole + placement = (f"sliced from the globally sorted list " + f"[{begin}:{end}]") if b_size > 1 else "whole basis" + else: + paths = [p for p in args.detfiles.split(',') if p] + if not paths: + print("ERROR: no determinant files given", file=sys.stderr) + return 1 + # When the files already partition the basis -- each sorted, disjoint, and + # in globally sorted order, as gdet emits and Fe4S4's four files are -- + # a rank can read only its own share and nothing is read twice. + # + # The share must be a CONTIGUOUS BLOCK of files, not a stride: shard i has + # to hold a strictly lower range than shard i+1, so with 24 files over 8 + # ranks it is files 0-2, 3-5, ... and NOT 0,8,16. A strided assignment + # would interleave the ranges and fail the sorted-disjoint check. + took_file_block = b_size > 1 and len(paths) % b_size == 0 + if took_file_block: + per_rank = len(paths) // b_size + first = shard_index * per_rank + mine = paths[first:first + per_rank] + placement = (f"file {first}" if per_rank == 1 else + f"files {first}-{first + per_rank - 1}") + placement += f" of {len(paths)}" + elif b_size > 1: + mine = paths + placement = "sliced from the globally sorted concatenation" + else: + mine = paths + placement = "whole basis" + det_strings = read_det_strings(mine) + bad = [s for s in det_strings if len(s) != total_bits] + if bad: + print(f"ERROR: determinant strings must be {total_bits} bits " + f"(2 * NORB); found one of length {len(bad[0])}", file=sys.stderr) + return 1 + if not det_strings: + print("ERROR: no determinants read", file=sys.stderr) + return 1 + det = sbd.sort_bitarray_array( + sbd.from_strings(det_strings, args.bit_length, total_bits)) + if b_size > 1 and not took_file_block: + begin, end = shard_bounds(det.shape[0], b_size, shard_index) + det = det[begin:end] + placement += f" [{begin}:{end}]" + source = f"{len(paths)} file(s): {', '.join(paths)}" + + if det.shape[0] == 0: + print("ERROR: no determinants read", file=sys.stderr) + return 1 + words = det.shape[1] + + config = sbd.GDB_SBD() + config.t_comm_size = args.t_comm_size + config.b_comm_size = args.b_comm_size + config.method = args.method + config.max_it = args.max_it + config.max_nb = args.max_nb + config.eps = args.eps + config.max_time = args.max_time + config.init = args.init + config.seed = args.seed + config.do_rdm = 1 if args.rdm_output else 0 + config.bit_length = args.bit_length + config.carryover_type = args.carryover_type + config.ratio = args.ratio + config.threshold = args.threshold + config.heatbath_cutoff = args.heatbath_cutoff + config.heatbath_truncation = args.heatbath_truncation + config.heatbath_batch_size = args.heatbath_batch_size + + if rank == 0: + grid = args.t_comm_size * args.b_comm_size + helper = size // grid if grid and size % grid == 0 else 0 + print("Configuration:") + print(f" Device: {device}") + print(f" Method: {'Davidson' if args.method == 0 else 'Davidson + stored H'}") + print(f" Max iterations: {config.max_it}") + print(f" Tolerance: {config.eps}") + print(f" MPI ranks: {size}") + print(f" Decomposition: t_comm_size={args.t_comm_size} " + f"b_comm_size={args.b_comm_size} helper={helper}") + print(f" Placement: {args.determinant_distribution or 'equal-bra-a (default)'}") + if args.b_comm_size == 1 and size > 1: + print(" NOTE: b_comm_size=1 means every rank holds the whole basis and " + "the whole excitation lookup. Only b_comm_size shards memory -- the " + "helper dimension divides work but not storage.") + if device == 'gpu' and helper != 1: + print(" ERROR: GDB on the Thrust backend requires helper == 1, i.e. " + "t_comm_size * b_comm_size == ranks. Give every rank to " + "--b_comm_size (with --t_comm_size <= it).") + if device == 'gpu-omp': + print(" WARNING: the OMP-offload backend has no GDB kernels -- there is " + "no `omp target` code under include/sbd/chemistry/gdb/ at all. " + "This run pins a device and then diagonalizes on the host. Use " + "--device gpu (Thrust, NVIDIA only) for GDB on a GPU.") + print() + print("Subspace:") + print(f" Source: {source}") + print(f" This rank: {placement}") + print(f" Determinants on this rank: {det.shape[0]}") + print(f" Orbitals: {norb} ({total_bits} spin orbitals), electrons: {nelec}") + print(f" Packing: {words} word(s) of {args.bit_length} bits") + print() + print("Running GDB diagonalization...") + print() + + start = time.perf_counter() + results = sbd.gdb_diag( + fcidump, det, config, + loadname=args.loadname, savename=args.savename, + determinant_distribution=args.determinant_distribution, + determinant_grid_a=args.determinant_grid_a, + determinant_grid_b=args.determinant_grid_b, + ) + elapsed = time.perf_counter() - start + + if rank != 0: + return 0 + + print("=" * 70) + print("Results") + print("=" * 70) + print(f"Device: {device.upper()}") + print(f"Ground state energy: {results['energy']:.10f} Hartree") + print(f"Subspace dimension: {results['global_dim']} " + f"({results['local_dim']} on this rank)") + print(f"Placement applied: {results['determinant_distribution']}") + print(f"Wall time: {elapsed:.2f} s") + + density = results['density'] + combined = [density[2 * i] + density[2 * i + 1] for i in range(len(density) // 2)] + print(f"Density: {np.round(combined, 6).tolist()}") + print(f" (sums to {sum(combined):.6f}; should equal the electron count {nelec})") + + carryover = results['carryover_det'] + if args.carryover_type: + print(f"Carryover determinants on this rank: {carryover.shape[0]}") + if size > 1: + print(" NOTE: the carryover list is this rank's shard, not the whole " + "list -- " + "HeatbathExpansion ends with sort_global_bitarray + " + "redistribution_bitarray over the WORLD communicator " + "(gdb/expansion.h), so each rank holds a different slice and a " + "driver that feeds it back must allgather first.") + if args.carryover_type in (2, 3): + print(" The expanded list includes its parents: " + "local_heatbath_expansion starts from `edet = det` " + "(gdb/expansion.h:543), so this is the next subspace, not just " + "the new candidates.") + + if args.rdm_output: + # assemble_rdms' layout was verified against PySCF on the TPB path; the + # keys and flat-index convention are shared with GDB, but that + # verification has not been repeated here. + rdm1, rdm2 = assemble_rdms(results, norb) + if rdm1 is None: + print("No RDMs returned (do_rdm was 0)") + else: + np.savez(args.rdm_output, rdm1=rdm1, rdm2=rdm2) + print(f"\nRDMs saved to {args.rdm_output}") + print(f"1-RDM trace: {np.trace(rdm1):.6f} " + f"(should equal the electron count {nelec})") + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py new file mode 100755 index 0000000..3066a99 --- /dev/null +++ b/examples/gdb/run_gdb_heatbath.py @@ -0,0 +1,476 @@ +#!/usr/bin/env python3 + +# 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. + +""" +Grow a GDB subspace with SBD's own heatbath expansion, round after round. + +This is selected CI (HCI) driven from Python: diagonalize, let SBD expand the +subspace from the resulting wavefunction, diagonalize the larger subspace, repeat. +No new machinery is needed for the loop itself -- ``carryover_type`` 2 and 3 return +the parents *together with* the new candidates (``local_heatbath_expansion`` opens +with ``edet = det``), so the result of one round IS the next round's subspace. + +**The cutoff is the dial, and the loop converges per cutoff.** Expansion admits a +candidate when its estimated contribution exceeds ``heatbath_cutoff``, so for a +fixed cutoff the subspace reaches a self-consistent size and stops growing -- +further rounds buy nothing. Going deeper means lowering the cutoff. Hence a +*ladder*: spend rounds at one cutoff until the energy stops moving, then step to +the next. ``--cutoffs`` takes the whole ladder. + +Do not confuse ``--heatbath_truncation`` with the cutoff: it discards *parents* +by weight before expansion even starts, and its default of 0 (keep everything) is +almost always what you want. Setting it to 1e-4 on a 59,536-determinant Fe4S4 +wavefunction cut the subspace to 506. + +Seeds, so the same driver produces every row of a seed comparison: + + --seed files determinant files (default: upstream's four Fe4S4 files) + --seed hf the Hartree-Fock determinant alone, the no-input null + --seed from-alpha an alpha list interleaved with itself (a TPB-shaped space) + --seed strings an arbitrary bitstring file, e.g. sampled configurations + +Usage: + # Fe4S4 from upstream's shipped subspace, one cutoff + python run_gdb_heatbath.py --cutoffs 1e-3 + + # A ladder, stopping if the subspace would pass 2M determinants + python run_gdb_heatbath.py --cutoffs 1e-3,1e-4,1e-5 --max_dim 2000000 + + # The null hypothesis: no input subspace at all, just Hartree-Fock + python run_gdb_heatbath.py --seed hf --cutoffs 1e-3,1e-4 + + # Sharded across 4 ranks/GPUs. At --b_comm_size == ranks with t=1 the + # expansion comes back already sharded for the next round, no gather needed. + mpirun -np 4 python run_gdb_heatbath.py --b_comm_size 4 --cutoffs 1e-3,1e-4 + mpirun -np 4 python run_gdb_heatbath.py --b_comm_size 4 --device gpu --cutoffs 1e-4 + + # Record the (dimension, energy) series for a comparison table + python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log fe4s4_ladder.json +""" + +import argparse +import json +import sys +import time + +import numpy as np + +_UPSTREAM = "../../vendor/sbd-upstream" +_GDB_APP = f"{_UPSTREAM}/apps/chemistry_gdb_selected_basis_diagonalization" +_DEFAULT_FCIDUMP = f"{_GDB_APP}/fcidump_Fe4S4.txt" +_DEFAULT_DETFILES = ",".join(f"{_GDB_APP}/det{i}.txt" for i in range(4)) + + +def parse_args(): + """Parse command line arguments.""" + parser = argparse.ArgumentParser( + description="Iteratively grow a GDB subspace by heatbath expansion", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + + parser.add_argument('--device', default='cpu', + choices=['cpu', 'gpu', 'gpu-omp', 'auto'], + help="Compute backend. 'gpu' (Thrust) is the only GPU backend " + "with GDB kernels, and it requires the helper dimension " + "to be 1, i.e. t_comm_size * b_comm_size == ranks") + parser.add_argument('--fcidump', default=_DEFAULT_FCIDUMP, + help='FCIDUMP file defining the Hamiltonian') + + # --- the seed --------------------------------------------------------- + parser.add_argument('--seed', default='files', + choices=['files', 'hf', 'from-alpha', 'strings'], + help='Where the starting subspace comes from') + parser.add_argument('--detfiles', default=_DEFAULT_DETFILES, + help='--seed files: comma-separated files of 2*norb-bit ' + 'determinant strings') + parser.add_argument('--alpha-file', default='', dest='alpha_file', + help='--seed from-alpha: a norb-bit alpha determinant list, ' + 'interleaved with itself into the full product basis') + parser.add_argument('--alpha-limit', type=int, default=0, dest='alpha_limit', + help='--seed from-alpha: keep only the first N alpha strings ' + '(the product costs N^2 determinants)') + parser.add_argument('--strings-file', default='', dest='strings_file', + help='--seed strings: a file of 2*norb-bit determinant ' + 'strings, e.g. sampled configurations') + + # --- the ladder ------------------------------------------------------- + parser.add_argument('--cutoffs', default='1e-3', + help='Comma-separated heatbath cutoffs, smallest step last. ' + 'Each is a rung: rounds run at that cutoff until the ' + 'energy stops moving, then the next rung begins') + parser.add_argument('--max_rounds', type=int, default=8, + help='Maximum rounds per rung') + parser.add_argument('--energy_tol', type=float, default=1e-5, + help='Advance to the next rung once |dE| between rounds falls ' + 'below this (Hartree)') + parser.add_argument('--max_dim', type=int, default=0, + help='Stop before diagonalizing a subspace larger than this ' + '(0 = no cap). The cap is what keeps rows of a ' + 'comparison cost-matched') + parser.add_argument('--heatbath_truncation', type=float, default=0.0, + help='Weight threshold applied to PARENTS before expanding. ' + 'Leave at 0 unless you mean to expand only from the ' + 'dominant determinants') + parser.add_argument('--heatbath_batch_size', type=int, default=1000000, + help='Heatbath expansion batch size per rank') + + # --- diagonalization -------------------------------------------------- + parser.add_argument('--method', type=int, default=0, choices=[0, 1], + help='0=Davidson, 1=Davidson storing the Hamiltonian. GDB has ' + 'no Lanczos, so TPB methods 2 and 3 do not exist here') + parser.add_argument('--tolerance', '--eps', type=float, default=1e-6, + dest='eps', help='Davidson convergence tolerance') + parser.add_argument('--iteration', '--max_it', type=int, default=30, + dest='max_it', help='Maximum Davidson iterations per round') + parser.add_argument('--block', '--max_nb', type=int, default=10, + dest='max_nb', help='Maximum number of basis vectors') + parser.add_argument('--bit_length', type=int, default=64, + help='Bits per packed word') + parser.add_argument('--carryover_type', type=int, default=2, choices=[2, 3], + help='Heatbath variant: 2 or 3. Types 0 and 1 do not expand, ' + 'so they cannot drive this loop') + + # --- decomposition ---------------------------------------------------- + parser.add_argument('--b_comm_size', type=int, default=1, + help='Basis shards. The only dimension that divides memory') + parser.add_argument('--t_comm_size', type=int, default=1, + help='Task communicator size; must not exceed b_comm_size') + parser.add_argument('--determinant_distribution', default='', + choices=['', 'input', 'equal-bra-a', 'count', 'count-sorted', + 'grid-cyclic', 'grid-cyclic-balanced'], + help='How determinants are placed across b_comm') + + parser.add_argument('--log', default='', metavar='FILE', + help='Write the per-round (dimension, energy) series as JSON') + + return parser.parse_args() + + +def read_strings(paths): + """Read bitstrings from one or more text files, in file order.""" + out = [] + for path in paths: + with open(path, encoding='utf-8') as handle: + out.extend(line.strip() for line in handle if line.strip()) + return out + + +def interleave(alpha, beta): + """Interleave a norb-bit alpha and beta string into one 2*norb-bit determinant. + + Bit ``2 * i`` is alpha orbital ``i`` and bit ``2 * i + 1`` is beta orbital ``i``, + counting from the right -- the order the GDB app documents and ``from_strings`` + packs. + """ + a_rev, b_rev = alpha[::-1], beta[::-1] + return ''.join(a_rev[i] + b_rev[i] for i in range(len(a_rev)))[::-1] + + +def hartree_fock_string(norb, nelec, ms2): + """The Hartree-Fock determinant: the lowest orbitals doubly occupied. + + Built from the FCIDUMP header alone, so ``--seed hf`` needs no input subspace at + all. That is the point of it: if expansion from HF reaches the same place as + expansion from a sampled subspace, the sampling added nothing. + """ + n_alpha = (nelec + ms2) // 2 + n_beta = nelec - n_alpha + if n_alpha > norb or n_beta > norb: + raise ValueError( + f"cannot place {n_alpha} alpha and {n_beta} beta electrons in {norb} " + f"orbitals") + alpha = ''.join('1' if i < n_alpha else '0' for i in range(norb))[::-1] + beta = ''.join('1' if i < n_beta else '0' for i in range(norb))[::-1] + return interleave(alpha, beta) + + +def shard_bounds(total, b_comm_size, index): + """This rank's slice of a globally sorted list. + + Mirrors SBD's ``q = N/p``, remainder-to-the-low-ranks split + (``balanced_begin``, framework/bit_manipulation.h:1523-1531). + """ + quotient, remainder = divmod(total, b_comm_size) + begin = index * quotient + min(index, remainder) + return begin, begin + quotient + (1 if index < remainder else 0) + + +def build_seed(args, sbd, norb, nelec, ms2, rank): + """The starting determinant array for this rank, plus a description.""" + total_bits = 2 * norb + + if args.seed == 'hf': + strings = [hartree_fock_string(norb, nelec, ms2)] + source = f"Hartree-Fock determinant ({nelec} electrons, MS2={ms2})" + elif args.seed == 'from-alpha': + if not args.alpha_file: + raise ValueError("--seed from-alpha requires --alpha-file") + alpha = read_strings([args.alpha_file]) + if args.alpha_limit: + alpha = alpha[:args.alpha_limit] + bad = [s for s in alpha if len(s) != norb] + if bad: + raise ValueError( + f"--alpha-file strings must be {norb} bits, found {len(bad[0])}") + strings = [interleave(a, b) for a in alpha for b in alpha] + source = f"{args.alpha_file}: {len(alpha)} alpha -> {len(alpha)}^2 product" + else: + paths = ([args.strings_file] if args.seed == 'strings' + else [p for p in args.detfiles.split(',') if p]) + if not paths or not all(paths): + raise ValueError(f"--seed {args.seed} requires input file(s)") + strings = read_strings(paths) + source = f"{len(paths)} file(s): {', '.join(paths)}" + + bad = [s for s in strings if len(s) != total_bits] + if bad: + raise ValueError( + f"determinant strings must be {total_bits} bits (2 * NORB), " + f"found one of length {len(bad[0])}") + if not strings: + raise ValueError("the seed is empty") + + # Sorted globally here so that slicing gives the disjoint, ordered shards + # gdb_diag requires. + whole = sbd.sort_bitarray_array( + sbd.from_strings(strings, args.bit_length, total_bits, device=args.device), + device=args.device) + if args.b_comm_size == 1: + return whole, source + begin, end = shard_bounds(whole.shape[0], args.b_comm_size, + rank % args.b_comm_size) + return whole[begin:end], source + + +def reshard(det, args, comm, rank, size): + """Put the expanded list back into the shape the next round expects. + + ``HeatbathExpansion`` ends with a global sort and a count rebalance over the + WORLD communicator (gdb/expansion.h:821-822), so the result is already + globally sorted, balanced and disjoint across world ranks. When + ``b_comm_size == ranks`` and ``t_comm_size == 1`` the b_comm position equals the + world rank, so that is *exactly* a valid next-round shard and nothing needs to + move -- which is the configuration a GPU run uses anyway, since Thrust requires + the helper dimension to be 1. + + Any other grid means world position and b_comm position disagree, and the only + portable fix from Python is to gather and re-slice. That materializes the whole + list on every rank, so it is the slow path and says so. + """ + if args.b_comm_size == 1 or size == 1: + if size == 1: + return det + # Every rank must hold the whole basis; the shards are disjoint pieces of it. + gathered = comm.allgather(det) + return np.concatenate([g for g in gathered if g.size], axis=0) + + if args.b_comm_size == size and args.t_comm_size == 1: + return det # already aligned: world rank == b_comm position + + if rank == 0: + print(" NOTE: b_comm_size != ranks (or t_comm_size > 1), so the expanded " + "list has to be gathered and re-sliced, which puts the whole list on " + "every rank. Use --b_comm_size == ranks with --t_comm_size 1 to avoid " + "it.", flush=True) + gathered = comm.allgather(det) + whole = np.concatenate([g for g in gathered if g.size], axis=0) + begin, end = shard_bounds(whole.shape[0], args.b_comm_size, + rank % args.b_comm_size) + return whole[begin:end] + + +def main(): + args = parse_args() + + import sbd + from mpi4py import MPI + + sbd.init(device=args.device) + backend = sbd.get_backend(args.device) + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + size = comm.Get_size() + + fcidump = backend.LoadFCIDump(args.fcidump) + header = fcidump.header + norb = int(header["NORB"]) + nelec = int(header["NELEC"]) + ms2 = int(header.get("MS2", 0)) + + cutoffs = [float(c) for c in args.cutoffs.split(',') if c.strip()] + if not cutoffs: + print("ERROR: --cutoffs is empty", file=sys.stderr) + return 1 + + try: + det, source = build_seed(args, sbd, norb, nelec, ms2, rank) + except ValueError as exc: + if rank == 0: + print(f"ERROR: {exc}", file=sys.stderr) + return 1 + + def global_dim(local): + return comm.allreduce(int(local.shape[0]), op=MPI.SUM) \ + if args.b_comm_size > 1 else int(local.shape[0]) + + # Collective: every rank must call it, so it cannot live inside a rank guard. + # This deadlocks otherwise -- rank 0 waits in the allreduce while the others walk + # into gdb_diag's own collectives. It hides in serial testing because the + # b_comm_size == 1 path short-circuits without communicating. + starting_dim = global_dim(det) + if args.b_comm_size > 1 and det.shape[0] == 0: + # Printed by whichever rank is starved, which is the informative one. + print(f" NOTE: rank {rank}'s shard is empty -- the seed has fewer " + f"determinants than b_comm_size. 'count' and 'count-sorted' placement " + f"cannot take an empty shard; 'input' and the grid schemes can.", + flush=True) + + if rank == 0: + grid = args.t_comm_size * args.b_comm_size + helper = size // grid if grid and size % grid == 0 else 0 + print("=" * 78) + print("GDB + heatbath expansion - iterative subspace growth") + print("=" * 78) + print(f" Device: {sbd.get_device()} ranks: {size}") + print(f" Decomposition: t={args.t_comm_size} b={args.b_comm_size} " + f"helper={helper}") + print(f" System: {norb} orbitals ({2 * norb} spin orbitals), " + f"{nelec} electrons") + print(f" Seed: {source}") + print(f" Starting dimension: {starting_dim}") + print(f" Cutoff ladder: {cutoffs}") + print(f" Per rung: up to {args.max_rounds} rounds, advance when " + f"|dE| < {args.energy_tol:g}") + if args.max_dim: + print(f" Dimension cap: {args.max_dim}") + if args.heatbath_truncation: + print(f" WARNING: --heatbath_truncation {args.heatbath_truncation:g} " + f"discards parents BEFORE expanding, which can shrink the " + f"subspace rather than grow it.") + print() + print(f" {'rung':>4} {'round':>5} {'dimension':>12} {'energy':>18} " + f"{'dE':>12} {'secs':>7}") + + history = [] + previous_energy = None # for the displayed dE, continuous across rungs + stop_reason = "ladder complete" + + for rung, cutoff in enumerate(cutoffs): + rung_converged = False + # Deliberately NOT seeded from the previous rung's energy. Round 0 of a new + # rung re-diagonalizes the subspace the previous rung already expanded, so + # its energy is essentially unchanged -- testing convergence against that + # would declare the rung done before its smaller cutoff had expanded + # anything at all. + rung_baseline = None + for round_index in range(args.max_rounds): + dim = global_dim(det) + if args.max_dim and dim > args.max_dim: + stop_reason = f"dimension {dim} exceeds --max_dim {args.max_dim}" + if rank == 0: + print(f" stopping: {stop_reason}") + return _finish(args, history, stop_reason, rank) + + config = backend.GDB_SBD() + config.method = args.method + config.eps = args.eps + config.max_it = args.max_it + config.max_nb = args.max_nb + config.bit_length = args.bit_length + config.b_comm_size = args.b_comm_size + config.t_comm_size = args.t_comm_size + config.carryover_type = args.carryover_type + config.heatbath_cutoff = cutoff + config.heatbath_truncation = args.heatbath_truncation + config.heatbath_batch_size = args.heatbath_batch_size + + start = time.perf_counter() + result = sbd.gdb_diag( + fcidump, det, config, device=args.device, + determinant_distribution=args.determinant_distribution, + ) + elapsed = time.perf_counter() - start + + energy = result["energy"] + delta = None if previous_energy is None else energy - previous_energy + if rank == 0: + shown = "" if delta is None else f"{delta:+.6f}" + print(f" {rung:>4} {round_index:>5} {dim:>12} {energy:>18.10f} " + f"{shown:>12} {elapsed:>7.1f}", flush=True) + history.append({ + "rung": rung, "cutoff": cutoff, "round": round_index, + "dimension": dim, "energy": energy, + "delta_energy": delta, "seconds": elapsed, + }) + + grown = result["carryover_det"] + det = reshard(grown, args, comm, rank, size) + new_dim = global_dim(det) + if rank == 0 and new_dim < dim: + print(f" NOTE: the expansion returned fewer determinants " + f"({new_dim} < {dim}). With --heatbath_truncation > 0 the " + f"parents are pruned before expanding.") + + # Energy-targeted advance, which is the standard HCI protocol: a rung is + # done when another round at the same cutoff no longer moves the answer. + rung_delta = None if rung_baseline is None else energy - rung_baseline + rung_baseline = energy + previous_energy = energy + if rung_delta is not None and abs(rung_delta) < args.energy_tol: + rung_converged = True + break + # A rung is equally done if the subspace stopped growing: the expansion + # has reached its fixed point for this cutoff. + if new_dim == dim: + rung_converged = True + if rank == 0: + print(f" rung {rung} at cutoff {cutoff:g}: subspace stopped " + f"growing at {dim} -- its fixed point for this cutoff") + break + + if rank == 0 and not rung_converged: + print(f" rung {rung} at cutoff {cutoff:g}: hit --max_rounds " + f"{args.max_rounds} without converging") + + return _finish(args, history, stop_reason, rank) + + +def _finish(args, history, stop_reason, rank): + """Print the summary and optionally write the series.""" + if rank != 0: + return 0 + if not history: + print("no rounds completed") + return 1 + first, last = history[0], history[-1] + print() + print("=" * 78) + print(f" {stop_reason}") + print(f" Start: dim {first['dimension']}, E = {first['energy']:.10f}") + print(f" End: dim {last['dimension']}, E = {last['energy']:.10f}") + print(f" Gain: {last['energy'] - first['energy']:+.6f} Hartree over " + f"{len(history)} round(s), " + f"{sum(h['seconds'] for h in history):.1f} s total") + print(" The energy is variational, so it must not rise: a positive gain means " + "the subspace shrank or a round failed to converge.") + if args.log: + with open(args.log, "w", encoding="utf-8") as handle: + json.dump({"stop_reason": stop_reason, "rounds": history}, handle, + indent=2) + print(f" Series written to {args.log}") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/examples/tpb/README.md b/examples/tpb/README.md index 3efe7d3..1794ea7 100644 --- a/examples/tpb/README.md +++ b/examples/tpb/README.md @@ -408,4 +408,6 @@ phases, not a fault. ## See Also +- [`../gdb/README.md`](../gdb/README.md) — the general determinant basis, which + decomposes over MPI differently - [Repository README](../../README.md) — Installation, API reference From f62c479b26438cf7016a5257d227fc09078f4c14 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 13:13:02 -0400 Subject: [PATCH 04/22] test: make the backend fixture and the package default agree On an install with more than one backend built, every test that builds objects from the fixture and then calls a module-level entry point failed: TypeError: tpb_diag(): incompatible function arguments. Invoked with: ..., , , ... The fixture calls get_backend(None), which auto-resolves and PREFERS Thrust when it is built and a GPU is present. The module-level tpb_diag/gdb_diag then reach _ensure_initialized(), whose default is 'cpu'. So the FCIDump and config came from one extension module and the diagonalization was dispatched in another. Pre-existing and not GDB-specific: test_reference_energies.py fails 2/2 the same way. It stayed invisible because a CPU-only build has nothing to disagree with, and that is what CI and a laptop give you. It shows up on any machine where the default build produced cpu + gpu + gpu-omp, i.e. exactly the hardware the GPU paths are meant for. Found while confirming the g++ CPU build on a box with both cpu and Thrust modules present. Pinning the default to whatever the fixture hands out fixes both test files and keeps SBD_TEST_DEVICE working as before. Co-Authored-By: Claude Opus 5 (1M context) --- test/conftest.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/conftest.py b/test/conftest.py index 49040d2..ad666dc 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -100,7 +100,16 @@ def backend(): """The SBD backend module to test, for tests calling the extension directly.""" import sbd - return sbd.get_backend(_requested_device()) + device = _requested_device() + # Pin the package default to the same backend this fixture hands out, or the two + # disagree on an install with more than one backend built. get_backend(None) here + # auto-resolves and PREFERS Thrust, while the module-level entry points reach + # _ensure_initialized() whose default is 'cpu' -- so a test would build its + # FCIDump and config from one module and call tpb_diag/gdb_diag in another, + # failing with "incompatible function arguments". Invisible on a CPU-only build, + # which is why it went unnoticed. + sbd.init(device=device or sbd.get_device()) + return sbd.get_backend(device) @pytest.fixture(scope="session") From 6fa41431f42d71441a567f7039f5749a3035b20f Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 13:20:16 -0400 Subject: [PATCH 05/22] README: frame both bases, and stop implying GDB plugs into the SQD addon Three gaps a reader would hit, now that GDB is more than a footnote. The Overview presented the bindings as TPB's, with GDB in a trailing sentence. It now names both methods by the property that decides which you want -- TPB's subspace is a product of two half-determinant lists, GDB's is the explicit list you pass and so can be sparse -- and says outright to prefer TPB for a product subspace, since it represents that case with the two lists rather than every product element and carries no extra constraints. The qiskit-addon-sqd section read as though either solver could be plugged in. It cannot, and not for want of plumbing: the addon's interface is a product subspace by construction (ci_strings is a (strings_a, strings_b) pair, SCIState.amplitudes an |a| x |b| matrix), so a sparse determinant list cannot be expressed without padding back to the full product and discarding the reason to use GDB. Stated plainly, with a pointer to calling gdb_diag directly. The API reference mentioned gdb_diag without its data contract. It now says det may be a whole basis or this rank's shard depending on b_comm_size, lists the constraints a sharded run carries, and notes they are checked rather than assumed -- then defers to examples/gdb/README.md for the decomposition, the placement schemes and which returned values are replicated. Deliberately no copy of that detail here: duplicating it guarantees the two drift. Also corrects the savename description, which this branch had made wrong: with the basis split it is one file per b_comm position, not a single ...000000.bin. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 52 +++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index 5f0ad5b..ffbfde8 100644 --- a/README.md +++ b/README.md @@ -4,10 +4,24 @@ Python bindings for the Selected Basis Diagonalization (SBD) library, with CPU a ## Overview -SBD (Selected Basis Diagonalization) is a high-performance library for quantum chemistry calculations. The Python bindings provide access to SBD's **Tensor-Product Basis (TPB)** diagonalization method on CPU and GPU. +SBD (Selected Basis Diagonalization) is a high-performance library for quantum +chemistry calculations. These bindings expose two of its diagonalization methods on CPU +and GPU. They differ in the shape of the subspace they span, which is what decides +which one you want: + +- **Tensor-Product Basis (TPB)** — the subspace is the Cartesian product of an alpha + and a beta determinant list, so its dimension is `|adet| × |bdet|`. The mature path, + and the one the qiskit-addon-sqd integration uses. +- **General-Determinant Basis (GDB)** — the subspace is the explicit determinant list + you pass, so it can be an arbitrary *sparse* set rather than a product. Experimental: + newer, a smaller tested surface, and still evolving as upstream SBD does. + +If your subspace is a product, prefer TPB — it represents that case with two +half-determinant lists instead of every product element, and needs no extra +constraints. Reach for GDB when the subspace is not a product. **Key Features:** -- **TPB diagonalization** for quantum chemistry Hamiltonians +- **TPB and GDB diagonalization** for quantum chemistry Hamiltonians - Three backends, selected per call at runtime via `device=`: `'cpu'` (host OpenMP), `'gpu'` (NVHPC Thrust/CUDA, NVIDIA only) and `'gpu-omp'` (OpenMP target offload, **NVIDIA or AMD**). All the backends your @@ -16,7 +30,9 @@ SBD (Selected Basis Diagonalization) is a high-performance library for quantum c - MPI parallelization - Integration with [qiskit-addon-sqd](https://github.com/Qiskit/qiskit-addon-sqd) for SQD workflows -In addition to TPB, this package also contains experimental support for SBD's **General-Determinant Basis (GDB)** method. However, SBD's **Creation/Annihilation operator (CAOP)** method is currently not supported by this wrapper; users who need it should reference and use the C++ CLI apps in the upstream submodule (`vendor/sbd-upstream/apps/`). +SBD's **Creation/Annihilation operator (CAOP)** method is not supported by this +wrapper; users who need it should reference and use the C++ CLI apps in the upstream +submodule (`vendor/sbd-upstream/apps/`). > [!NOTE] > This package is newly open-sourced. The Python API follows semantic versioning, but the build configuration and GPU backends have been exercised on a limited set of platforms — please report issues. @@ -66,6 +82,15 @@ matter for it. SBD can serve as the eigensolver backend for qiskit-addon-sqd's SQD workflow. +> [!NOTE] +> This integration is **TPB-only**, and not merely for want of plumbing. +> `solve_sci`/`solve_sci_batch` call `tpb_diag`, and the addon's interface describes a +> product subspace by construction: `ci_strings` is a `(strings_a, strings_b)` pair and +> `SCIState.amplitudes` is an `|a| × |b|` matrix. A sparse determinant list cannot be +> expressed that way without padding back up to the full product, which discards the +> reason to use GDB. For GDB, call `gdb_diag` directly or use the drivers in +> [`examples/gdb/`](examples/gdb/README.md). + **Note:** Requires [qiskit-addon-sqd](https://github.com/Qiskit/qiskit-addon-sqd) with distributed (SPMD) support — `diagonalize_fermionic_hamiltonian` calling `sci_solver` on every MPI rank. This is available in `qiskit-addon-sqd` version `0.13.1` or higher. ### Plain SQD @@ -253,11 +278,24 @@ results = sbd.gdb_diag(fcidump, det, sbd_data, spin-beta orbital `i`. The determinants must be distinct; they are sorted into SBD's canonical order internally, which `sort_bitarray` reproduces. +`det` may be a single `(ndets, words)` array or **this rank's shard of one**: +`sbd_data.b_comm_size` decides which. At `1` every rank passes the whole basis; above +`1` every rank passes its own shard and the union over b_comm positions is the basis, +which is the only way GDB's memory scales. Sharded runs carry further constraints — +`t_comm_size ≤ b_comm_size`, their product dividing the rank count, a helper dimension +of 1 on the Thrust backend, and shards that are globally sorted and disjoint. All are +checked and raise rather than silently diagonalizing the wrong subspace. See +[`examples/gdb/README.md`](examples/gdb/README.md) for the decomposition, the six +determinant-placement schemes, and which returned values are replicated versus +sharded. + `gdb_diag` does not return the wavefunction amplitudes, because SBD's `gdb::diag` -has no in-memory output for them. Passing `savename` makes SBD write them to -`f"{savename}000000.bin"` instead: two `size_t` headers -(`n_dets`, `words_per_det`), then `n_dets × words_per_det` `size_t` determinant -words in canonical order, then `n_dets` `float64` amplitudes. +has no in-memory output for them. Passing `savename` makes SBD write them instead, as +one file per b_comm position — `f"{savename}{rank_b:06d}.bin"`, so just +`…000000.bin` when `b_comm_size` is 1 and `b_comm_size` files otherwise, each holding +only that shard. Each file is two `size_t` headers (`n_dets`, `words_per_det`), then +`n_dets × words_per_det` `size_t` determinant words in canonical order, then `n_dets` +`float64` amplitudes. The optional `device` parameter overrides the default set by `init()`. From c1e411b88af6b1896d61dfb69530ae913a8f7f4a Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 13:40:36 -0400 Subject: [PATCH 06/22] examples/gdb/README: merge a duplicated paragraph in the GPU section Stripping the benchmark figures replaced the passage from "Verified on 8x H100..." onward, but the paragraph before it already made the same point, so the section stated "multi-GPU GDB requires --b_comm_size" twice in consecutive paragraphs. Merged into one, keeping the mult_thrust.h reference and the worked example. The substance was right and is unchanged: GDB has Thrust kernels only, so GPU GDB is NVIDIA-only and AMD GDB is CPU-only; the Thrust path needs helper == 1, which is why more than one GPU requires splitting the basis. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 7a22a20..8c52a5f 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -282,15 +282,12 @@ pins a device and then diagonalizes on the host; the driver warns when it resolv there. AMD has no Thrust build, so AMD GDB is CPU-only. The Thrust path additionally requires **`helper == 1`** (`gdb/mult_thrust.h:310-314`, -checked on every kernel launch), which means every rank must go to -`t_comm_size × b_comm_size`. Since `helper = ranks / (t × b)`, leaving the basis in one -block caps GPU GDB at a single rank — so **multi-GPU GDB requires `--b_comm_size`**. - -Because `h = ranks / (t · b)`, leaving the basis in a single block caps GPU GDB at one -rank — the helper dimension takes every rank you add. So **multi-GPU GDB requires -`--b_comm_size`**, e.g. `-np 4 --b_comm_size 4`. Asking for more ranks than `t · b` is -refused up front, naming the helper dimension, rather than throwing from inside a -kernel launch. +checked on every kernel launch), so every rank must go to `t_comm_size × b_comm_size`. +Since `helper = ranks / (t × b)`, leaving the basis in one block caps GPU GDB at a +single rank — the helper dimension absorbs every rank you add. So **multi-GPU GDB +requires `--b_comm_size`**, e.g. `-np 4 --b_comm_size 4`. Asking for more ranks than +`t × b` is refused up front, naming the helper dimension, rather than throwing from +inside a kernel launch. Only the Davidson runs on the device. `gdb/expansion.h` and `gdb/carryover.h` contain no `thrust::` code and are included outside any `SBD_THRUST` guard (`inc_all.h:23`), so From 4fd5639da5377d6c4487737531d9582da2691b4c Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 13:49:00 -0400 Subject: [PATCH 07/22] README: say which backends each solver actually has The backend bullet listed cpu, gpu and gpu-omp without qualification, which is right for TPB and misleading for GDB: GDB has Thrust kernels only, so under 'gpu-omp' it pins a device and then runs on the host. A reader of the landing page alone would have concluded AMD GDB has GPU support. Adds the distinction and a pointer, without pulling the GPU detail up from examples/gdb/README.md. That was the last place in the docs where the claim was ambiguous; the GDB and shared example READMEs already stated it, and the driver warns at runtime when --device resolves to gpu-omp. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index ffbfde8..94f2de7 100644 --- a/README.md +++ b/README.md @@ -26,7 +26,9 @@ constraints. Reach for GDB when the subspace is not a product. `'cpu'` (host OpenMP), `'gpu'` (NVHPC Thrust/CUDA, NVIDIA only) and `'gpu-omp'` (OpenMP target offload, **NVIDIA or AMD**). All the backends your toolchain supports can be built into one install; each is imported only when - first used + first used. **TPB runs on all three; GDB has Thrust kernels only**, so GDB on a + GPU is NVIDIA-only and under `'gpu-omp'` it falls back to the host — see + [`examples/gdb/README.md`](examples/gdb/README.md) - MPI parallelization - Integration with [qiskit-addon-sqd](https://github.com/Qiskit/qiskit-addon-sqd) for SQD workflows From 36fdd6264f71e338c54d09784ba8c9d1bbd416fa Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 14:04:39 -0400 Subject: [PATCH 08/22] run_gdb_heatbath: stop over-counting the subspace when shards are replicated global_dim() summed local determinant counts over MPI_COMM_WORLD. Ranks sharing a b_comm position hold identical shards by contract, so with t_comm_size or the helper dimension above 1 every shard was counted once per replica: Fe4S4's 59,536-determinant basis reported as 119,072 at --b_comm_size 2 --t_comm_size 2 on four ranks. Energies were unaffected -- the binding computes its own global_dim over b_comm, which is what run_gdb_diag reports -- but the driver's dimension column was wrong and --max_dim compared against the inflated figure, so a ladder would have stopped early. It hid because every run until now used b == ranks with t = 1, where the world sum and the b_comm sum coincide. Found by deliberately exercising the b != ranks reshard path. Fixed by contributing only from the ranks that are one-per-b-position: with the layout rank = h*(b*t) + t_index*b + b_index those are exactly ranks 0..b_comm_size-1, so the b_comm sum needs no sub-communicator. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/run_gdb_heatbath.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py index 3066a99..3b4cc5f 100755 --- a/examples/gdb/run_gdb_heatbath.py +++ b/examples/gdb/run_gdb_heatbath.py @@ -321,8 +321,20 @@ def main(): return 1 def global_dim(local): - return comm.allreduce(int(local.shape[0]), op=MPI.SUM) \ - if args.b_comm_size > 1 else int(local.shape[0]) + """Determinants across the whole basis, counted once each. + + Summing over the world communicator would over-count: ranks sharing a b_comm + position hold identical shards, so with t_comm_size or the helper dimension + above 1 each shard is counted once per replica. Contribute only from the ranks + that are one-per-b-position -- with the layout + ``rank = h*(b*t) + t_index*b + b_index``, those are exactly ranks + ``0 .. b_comm_size - 1`` -- which gives the b_comm sum without building a + sub-communicator. + """ + if args.b_comm_size == 1: + return int(local.shape[0]) + mine = int(local.shape[0]) if rank < args.b_comm_size else 0 + return comm.allreduce(mine, op=MPI.SUM) # Collective: every rank must call it, so it cannot live inside a rank guard. # This deadlocks otherwise -- rank 0 waits in the allreduce while the others walk From bf7f3607a20ca9b8cb4d233305b78bc9afc40a4b Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 14:49:43 -0400 Subject: [PATCH 09/22] README: the GDB config table still claimed b_comm_size must be 1 This branch lifts that restriction, and the prose further down already described the shard contract, but the configuration table was left saying "must be 1 for gdb_diag" -- actively wrong for the branch that changes it. I had edited the surrounding prose and missed the table. Replaced with what the three dimensions now are and whose constraint each one is, since they are different in kind: b_comm_size is free and is the only dimension that divides memory; t_comm_size <= b_comm_size is upstream's algorithm (one task per basis-ring station, unchecked there, so gdb_diag rejects it up front); and their product must divide the rank count, with the derived helper dimension taking the quotient and required to be 1 on Thrust, which is why more than one GPU needs a split basis. Depth stays in examples/gdb/README.md. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 38 +++++++++++++++++++++----------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/README.md b/README.md index 94f2de7..3b6e9e5 100644 --- a/README.md +++ b/README.md @@ -228,31 +228,35 @@ and replaces the determinant communicators with a single basis communicator: | Attribute | Default | Description | |-----------|---------|-------------| -| `b_comm_size` | 1 | Basis communicator size — must be 1 for `gdb_diag`, see below | -| `t_comm_size` | 1 | Task communicator size — must be 1 while `b_comm_size` is, see below | +| `b_comm_size` | 1 | Basis communicator size, i.e. how many shards the determinant list is split into — the only dimension that divides memory | +| `t_comm_size` | 1 | Task communicator size — must not exceed `b_comm_size`, see below | | `seed` | 1729 | Seed for the initial vector | | `heatbath_cutoff` | 1e-4 | Heatbath expansion cutoff | | `heatbath_truncation` | 0.0 | Weight truncation applied before heatbath expansion | | `heatbath_batch_size` | 200000000 | Heatbath expansion batch size | -**Why both must be 1**, since the two constraints have different owners: +**How the three GDB dimensions relate**, and whose constraint each one is: -`b_comm_size == 1` is a limitation of *this wrapper*, not of SBD. Upstream's in-memory -`gdb::diag` expects each rank to pass **its own shard** of the determinant list — that -is what upstream's file-based entry point hands it, after distributing determinant -files across `b_comm`. This wrapper passes the whole list from every rank, which is -only consistent with a single basis block, so it rejects anything else rather than -have each rank diagonalize the full subspace while believing it held a shard. Upstream -itself runs with a split basis: its own `run.sh` for the GDB app passes -`--b_comm_size 2`. +`b_comm_size` splits the determinant list across ranks and is the only dimension that +divides memory — the other two divide work. Above 1, every rank passes its own shard +rather than the whole list (see `gdb_diag` below). It was pinned to 1 before this +wrapper could shard an in-memory list; upstream never required it, and its own `run.sh` +for the GDB app runs `--b_comm_size 2` through the file-based path. -`t_comm_size == 1` then follows from *upstream's* algorithm rather than from us. GDB's +`t_comm_size ≤ b_comm_size` is *upstream's* algorithm, not a wrapper choice. GDB's matrix-vector product rotates the ket around `b_comm` as a ring, so there are exactly -`b_comm_size` ring stations and one "task" is one station — meaning `t_comm_size` -cannot exceed `b_comm_size`. With the basis in a single block there is a single task. - -Ranks are not wasted in the meantime: the derived helper dimension, -`ranks / (t_comm_size × b_comm_size)`, absorbs them and does not change the energy. +`b_comm_size` ring stations and one "task" is one station — there cannot be more task +ranks than stations. Upstream does not check this, and exceeding it faults inside +helper construction, so `gdb_diag` rejects it up front. + +`t_comm_size × b_comm_size` must divide the rank count exactly. The helper dimension is +the quotient, `ranks / (t_comm_size × b_comm_size)`; it is derived rather than settable +and divides work without dividing memory. On the Thrust backend it must be 1, because +the GPU kernels never implemented that dimension — which is why more than one GPU +requires `b_comm_size > 1`. + +See [`examples/gdb/README.md`](examples/gdb/README.md) for the shard contract, the six +determinant-placement schemes and worked rank layouts. ### Diagonalization From fe142753184ba1fa984baa32482ceb578b9bd49a Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 15:04:25 -0400 Subject: [PATCH 10/22] README: the GDB GPU troubleshooting entry is a fix on this branch, not a limitation The entry inherited from the reorg branch says multi-rank GPU GDB is impossible because b_comm_size is pinned to 1. This branch lifts that, so the same symptom now has a resolution: give every rank to the basis, e.g. -np 4 with b_comm_size 4, leaving a helper dimension of 1. Also notes that gdb_diag checks it up front rather than letting the Thrust kernel throw mid-launch. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 3b6e9e5..0d17b88 100644 --- a/README.md +++ b/README.md @@ -341,13 +341,13 @@ prefer it. Tracked as **Ranks die with `Bus error` or `SIGSEGV` inside the MPI's own copy path** (`MPIR_Localcopy`, `ucp_worker_progress`, ...) **on a GPU backend:** the MPI is not GPU-aware and was handed a device pointer. Rebuild UCX `--with-cuda` / `--with-rocm`, and confirm with `ucx_info -d | grep -i 'Transport: cuda'` (or `rocm`). Two things mislead here. A partly GPU-aware stack fails in only one place: an MPICH with GPU support *disabled* over a CUDA-aware UCX ran OMP-offload fine and crashed only in Thrust, because the inter-rank path went through UCX while the local-copy path did not. And on AMD a non-ROCm-aware MPI does not crash at all — ROCm maps device memory into the process address space, so the host copy succeeds and merely stages everything through the host aperture (measured on MI250X, XNACK off, 8 ranks) — so a working AMD run is not evidence that the MPI is ROCm-aware. -**`GDB Thrust mult does not support h_comm_size > 1` from `gdb_diag` on more than one -rank:** GDB's Thrust kernels never implemented the helper dimension, and the helper -dimension is `ranks / (t_comm_size × b_comm_size)`. Since `gdb_diag` requires -`b_comm_size == 1` (see [Configuration](#configuration)), which forces `t_comm_size` to -1, every rank you add lands in the helper dimension — so GPU GDB is limited to a single -rank in this release. Run GDB on one GPU, or on the CPU backend, where the helper -dimension is unconstrained. TPB is unaffected and shards over `adet_comm_size` / -`bdet_comm_size` as usual. +**GDB on more than one GPU refuses to start, naming the helper dimension:** GDB's Thrust +kernels never implemented that dimension, and it is the quotient +`ranks / (t_comm_size × b_comm_size)` — so any rank you do not assign to `t` or `b` lands +there. Give every rank to the basis: `-np 4` with `b_comm_size = 4` leaves a helper +dimension of 1. Leaving `b_comm_size` at 1 puts *all* ranks in the helper dimension, +which is why multi-GPU GDB requires a split basis. `gdb_diag` checks this before doing +any work rather than letting the kernel throw mid-launch. The CPU backend has no such +restriction, and TPB is unaffected. **Repository:** https://github.com/Qiskit/sbd-eigensolver-python From 74835dd69579615e4952575ffbc558f6bcc34162 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 15:29:38 -0400 Subject: [PATCH 11/22] examples/gdb/README: stand alone now that the shared README is gone The intro and See Also deferred to examples/README.md for backend selection, which no longer exists. Replaced with what a GDB user needs in place: the --device values, with 'gpu' flagged as the only GPU backend that has GDB kernels; available_backends() and loaded_backends(); and the 'cpu'/'gpu-omp' interaction, where the CPU module initializes the shared OpenMP runtime host-only and the offload backend then runs silently on the host. Links to the TPB README for the bundled test data rather than restating it. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 29 ++++++++++++++++++++++++----- 1 file changed, 24 insertions(+), 5 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 8c52a5f..f243ae0 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -6,9 +6,7 @@ Cartesian product of an alpha and a beta list that TPB uses. TPB's dimension is makes it the right solver for an arbitrary sparse subspace — a set of sampled bitstrings used *as sampled*, with no product completion. -For backend selection, `--device` values and bundled test data, see -[`../README.md`](../README.md). For TPB and the SQD loops, see -[`../tpb/README.md`](../tpb/README.md). +For TPB and the SQD loops, see [`../tpb/README.md`](../tpb/README.md). ## run_gdb_diag.py — standalone GDB diagonalization @@ -274,6 +272,28 @@ For Fe4S4 the four shipped files are already the balanced globally-sorted split `--b_comm_size 4` the driver hands file *i* to shard *i* and `input` placement is already count-balanced. +### Choosing a backend + +`--device` selects the compute backend per run: + +```bash +--device cpu # host OpenMP (default) +--device gpu # NVHPC Thrust, NVIDIA only -- the only GPU backend with GDB kernels +--device gpu-omp # OpenMP target offload; compiles for GDB but has no GDB kernels +--device auto # GPU if one is available, else CPU +``` + +`sbd.available_backends()` reports what this install actually built — a static scan, so +it is safe to call outside `mpirun` — and `sbd.loaded_backends()` reports what the +current process has imported. + +One interaction to know about: **`'cpu'` and `'gpu-omp'` must not be used in the same +process.** Both link the same OpenMP runtime and the CPU module is built without offload +support, so whichever loads first initializes that runtime; if it is the CPU backend, the +offload backend can no longer acquire a device and silently runs on the host — right +answers, exit 0, idle GPU. `'cpu'` and `'gpu'` (Thrust) coexist fine, since Thrust does +not route device work through OpenMP. + ### On GPUs GDB has **Thrust kernels only** (`--device gpu`, NVIDIA). There is no `omp target` @@ -314,6 +334,5 @@ return an uninitialized energy. ## See Also -- [`../README.md`](../README.md) — backend selection, test data, performance tips -- [`../tpb/README.md`](../tpb/README.md) — TPB and the SQD loops +- [`../tpb/README.md`](../tpb/README.md) — TPB, the SQD loops, and the bundled test data - [Repository README](../../README.md) — installation, API reference From 59b95ff266af56d55d971982717a1bf1d93f88e8 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 15:47:31 -0400 Subject: [PATCH 12/22] README: keep the GDB API section brief, matching TPB; detail lives in examples/gdb The Configuration section had grown three paragraphs on how the GDB communicators relate, all of which examples/gdb/README.md already covers -- its dimension table, the ring explanation for t <= b, "only b shards memory", and the Thrust helper restriction. Replaced with the same shape TPB uses: the rank arithmetic in one line, plus the one API fact a caller cannot do without (above b_comm_size 1, det is this rank's shard), and a pointer for the rest. The gdb_diag prose further down is trimmed the same way. Two stale claims fixed while there. The "Shares ... with TPB_SBD" line listed method, whose range differs -- GDB has only 0 and 1 -- so method now appears in the GDB table explicitly, as it does in TPB's. It also listed do_shuffle, which this branch removed from GDB_SBD because upstream parses it and never reads it. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 51 +++++++++++++++++---------------------------------- 1 file changed, 17 insertions(+), 34 deletions(-) diff --git a/README.md b/README.md index 0d17b88..469a49f 100644 --- a/README.md +++ b/README.md @@ -222,41 +222,27 @@ Total MPI ranks = `task_comm_size × adet_comm_size × bdet_comm_size`. config = sbd.GDB_SBD() ``` -Shares `method`, `max_it`, `max_nb`, `eps`, `max_time`, `init`, `do_shuffle`, -`do_rdm`, `carryover_type`, `ratio`, `threshold` and `bit_length` with `TPB_SBD`, -and replaces the determinant communicators with a single basis communicator: +Shares `max_it`, `max_nb`, `eps`, `max_time`, `init`, `do_rdm`, `carryover_type`, +`ratio`, `threshold` and `bit_length` with `TPB_SBD`, and replaces the determinant +communicators with a single basis communicator: | Attribute | Default | Description | |-----------|---------|-------------| -| `b_comm_size` | 1 | Basis communicator size, i.e. how many shards the determinant list is split into — the only dimension that divides memory | -| `t_comm_size` | 1 | Task communicator size — must not exceed `b_comm_size`, see below | +| `method` | 0 | 0=Davidson, 1=Davidson+Ham. GDB has no Lanczos, so TPB's 2 and 3 are rejected | +| `b_comm_size` | 1 | Basis communicator size — how many shards the determinant list is split into | +| `t_comm_size` | 1 | Task communicator size — must not exceed `b_comm_size` | | `seed` | 1729 | Seed for the initial vector | | `heatbath_cutoff` | 1e-4 | Heatbath expansion cutoff | | `heatbath_truncation` | 0.0 | Weight truncation applied before heatbath expansion | | `heatbath_batch_size` | 200000000 | Heatbath expansion batch size | -**How the three GDB dimensions relate**, and whose constraint each one is: +Total MPI ranks = `t_comm_size × b_comm_size × helper`, where the helper dimension is +the derived quotient and is not settable. With `b_comm_size > 1` each rank passes its +own shard of the determinant list rather than the whole list. -`b_comm_size` splits the determinant list across ranks and is the only dimension that -divides memory — the other two divide work. Above 1, every rank passes its own shard -rather than the whole list (see `gdb_diag` below). It was pinned to 1 before this -wrapper could shard an in-memory list; upstream never required it, and its own `run.sh` -for the GDB app runs `--b_comm_size 2` through the file-based path. - -`t_comm_size ≤ b_comm_size` is *upstream's* algorithm, not a wrapper choice. GDB's -matrix-vector product rotates the ket around `b_comm` as a ring, so there are exactly -`b_comm_size` ring stations and one "task" is one station — there cannot be more task -ranks than stations. Upstream does not check this, and exceeding it faults inside -helper construction, so `gdb_diag` rejects it up front. - -`t_comm_size × b_comm_size` must divide the rank count exactly. The helper dimension is -the quotient, `ranks / (t_comm_size × b_comm_size)`; it is derived rather than settable -and divides work without dividing memory. On the Thrust backend it must be 1, because -the GPU kernels never implemented that dimension — which is why more than one GPU -requires `b_comm_size > 1`. - -See [`examples/gdb/README.md`](examples/gdb/README.md) for the shard contract, the six -determinant-placement schemes and worked rank layouts. +See [`examples/gdb/README.md`](examples/gdb/README.md) for the decomposition, the shard +contract, the six determinant-placement schemes, and why the helper dimension must be 1 +on the Thrust backend. ### Diagonalization @@ -286,14 +272,11 @@ SBD's canonical order internally, which `sort_bitarray` reproduces. `det` may be a single `(ndets, words)` array or **this rank's shard of one**: `sbd_data.b_comm_size` decides which. At `1` every rank passes the whole basis; above -`1` every rank passes its own shard and the union over b_comm positions is the basis, -which is the only way GDB's memory scales. Sharded runs carry further constraints — -`t_comm_size ≤ b_comm_size`, their product dividing the rank count, a helper dimension -of 1 on the Thrust backend, and shards that are globally sorted and disjoint. All are -checked and raise rather than silently diagonalizing the wrong subspace. See -[`examples/gdb/README.md`](examples/gdb/README.md) for the decomposition, the six -determinant-placement schemes, and which returned values are replicated versus -sharded. +`1` every rank passes its own shard, the union over b_comm positions being the basis. +Sharded input must be globally sorted and disjoint; that and the rank-layout constraints +are checked, and raise rather than silently diagonalizing the wrong subspace. +[`examples/gdb/README.md`](examples/gdb/README.md) has the shard contract, the +placement schemes, and which returned values are replicated versus sharded. `gdb_diag` does not return the wavefunction amplitudes, because SBD's `gdb::diag` has no in-memory output for them. Passing `savename` makes SBD write them instead, as From da123aeed6008c232da8ba4dd48175831a73c1cb Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 15:54:56 -0400 Subject: [PATCH 13/22] examples/gdb/README: note that the heatbath loop needs no amplitudes gdb_diag does not return the wavefunction, which invites the question of how an iterative driver can work without it. It can because the selection happens inside SBD: the amplitudes drive it -- weight truncation keeps determinants by |c| and heatbath scoring is essentially |c_i . H_ij| -- but WeightTruncation and HeatbathExpansion consume them in C++ and return only the expanded determinant list, which is the next subspace. Nothing but determinants crosses the Python boundary. Also says where amplitudes would be needed: Python-side selection, as the TPB enlarge-subspace driver does, which for GDB would mean reading them back from the per-shard savename files. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index f243ae0..97d2816 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -129,6 +129,21 @@ different sizes from the same cutoff, so a cutoff-matched comparison mostly repo subspace size rather than seed quality. `--max_dim` and the `--log` series are there for that. +### The loop needs no amplitudes + +A natural question, since `gdb_diag` does not return the wavefunction: the ladder never +needs it. The amplitudes are what drive the selection — weight truncation keeps +determinants by `|c|`, and heatbath scoring is essentially `|c_i · H_ij|` — but SBD +consumes them internally (`WeightTruncation` then `HeatbathExpansion`, which takes the +coefficients as an input) and hands back only the expanded determinant list. That list +is the next subspace, so the loop closes with nothing but determinants crossing the +Python boundary. + +Where amplitudes *would* be needed is Python-side selection — deciding yourself which +determinants to expand, as `../tpb/run_sqd_enlarge_subspace_sbd.py` does for TPB. For +GDB that means reading them back from the per-shard `savename` files, since +`gdb::diag` has no in-memory amplitude output. + ### `--heatbath_truncation` is not the cutoff It discards **parents** by weight *before* expansion starts, and its default of 0 From d4f4b132e8374211f8919123c116b0c7855f5474 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 16:04:41 -0400 Subject: [PATCH 14/22] docs: amplitudes are optional; move the savename layout into examples/gdb The top README carried the per-shard file layout, which is GDB detail and belongs with the rest of it. It also read as a requirement when it is not: nothing in the normal flow needs the wavefunction -- energy, density and RDMs come back directly, a heatbath ladder iterates on carryover_det, and GDB is not wired into qiskit-addon-sqd, whose SCIState would be the usual consumer. Top README now says amplitudes are not returned, that this rarely matters, and where to look if you do want them. examples/gdb/README.md gains a short "Getting the amplitudes, if you want them" section holding the file layout, and its outputs table stops repeating it. Co-Authored-By: Claude Opus 5 (1M context) --- README.md | 12 +++++------- examples/gdb/README.md | 18 +++++++++++++++++- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 469a49f..361320d 100644 --- a/README.md +++ b/README.md @@ -278,13 +278,11 @@ are checked, and raise rather than silently diagonalizing the wrong subspace. [`examples/gdb/README.md`](examples/gdb/README.md) has the shard contract, the placement schemes, and which returned values are replicated versus sharded. -`gdb_diag` does not return the wavefunction amplitudes, because SBD's `gdb::diag` -has no in-memory output for them. Passing `savename` makes SBD write them instead, as -one file per b_comm position — `f"{savename}{rank_b:06d}.bin"`, so just -`…000000.bin` when `b_comm_size` is 1 and `b_comm_size` files otherwise, each holding -only that shard. Each file is two `size_t` headers (`n_dets`, `words_per_det`), then -`n_dets × words_per_det` `size_t` determinant words in canonical order, then `n_dets` -`float64` amplitudes. +`gdb_diag` returns no wavefunction amplitudes — `gdb::diag` has no in-memory output for +them — and for most uses none are needed: the energy, density and RDMs come back +directly, and an iterative heatbath run gets its next subspace from `carryover_det`. If +you do want the amplitudes, `savename` writes them to disk; the file layout is in +[`examples/gdb/README.md`](examples/gdb/README.md). The optional `device` parameter overrides the default set by `init()`. diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 97d2816..9c50b00 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -345,7 +345,23 @@ return an uninitialized energy. | `density`, `one_p_rdm`, `two_p_rdm` | replicated on every rank | | `carryover_det` | **this rank's shard**, as an `(n, words)` array. For `carryover_type` 1 it is split over `b_comm` *and duplicated* across the helper dimension, so gathering means one representative per `rank % b_comm_size`; for types 2 and 3 it is split over the world communicator with no duplication, so a plain allgather is correct | | `local_dim` / `global_dim` | determinants on this rank, and summed over `b_comm` | -| `savename` | one file per shard, `f"{savename}{rank_b:06d}.bin"` — rank 0's file is not the whole wavefunction, and GDB has no combined matrix-form dump | +| `savename` | optional, and off by default: writes the amplitudes to one file per shard rather than returning them — see below | + +### Getting the amplitudes, if you want them + +Usually you do not. `gdb::diag` has no in-memory output for the wavefunction, and +nothing in the normal flow needs it: the energy, density and RDMs are returned directly, +and a heatbath ladder takes its next subspace from `carryover_det` (see [The loop needs +no amplitudes](#the-loop-needs-no-amplitudes)). GDB is also not wired into +qiskit-addon-sqd, whose `SCIState` would be the usual consumer — that path is TPB-only. + +When you do want them — your own analysis, or a selection step written in Python — pass +`savename` and read the files back. SBD writes one per b_comm position, +`f"{savename}{rank_b:06d}.bin"`: a single `…000000.bin` at `b_comm_size 1`, otherwise +`b_comm_size` files each holding only that shard, so no single file is the whole +wavefunction. Each is two `size_t` headers (`n_dets`, `words_per_det`), then +`n_dets × words_per_det` `size_t` determinant words in canonical order, then `n_dets` +`float64` amplitudes. ## See Also From 7cc5d93e16a5193922f88d51942df4888a003b02 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 16:49:42 -0400 Subject: [PATCH 15/22] examples/gdb: document the default data and parameters, and test the drivers Three gaps, all found by reading the README as a new user would. The commands that pass no input flags looked like they ran on no data at all. They run upstream's Fe4S4 case, because --fcidump and --detfiles default to it; say so in a new "The default data" section, and show one command with the paths written out so the --detfiles syntax is copyable. Add --device gpu examples. GDB's only GPU backend is Thrust, which requires helper == 1, so multi-GPU GDB cannot omit --b_comm_size -- worth showing next to the CPU runs rather than only in the GPU section further down. Document run_gdb_heatbath.py's parameters as tables, grouped seed / ladder / solver, matching examples/tpb/README.md. Cross-checked against argparse: no invented flags, every option covered, aliases named. test/test_gdb_drivers.py is new: the drivers had no test, so argument parsing, determinant reading, the sharding arithmetic and the default paths could all break with every library test still green. Runs each driver as a subprocess, so __main__ and the defaults are exercised rather than bypassed. The h2o case asserts the same -76.0588897208 that test_gdb_equivalence.py anchors against TPB, tying driver to solver; the slow case asserts the flagless run really is Fe4S4 at 59,536 determinants, which is what keeps the README's claim honest. Both guardrails are covered, t > b especially -- unguarded it segfaults instead of erroring. 6 passed in 16 s by default, 7 in 37 s with --run-slow. Verified non-vacuous: perturbing the anchor by 1e-8 fails the test. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 67 ++++++++++- test/test_gdb_drivers.py | 232 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 296 insertions(+), 3 deletions(-) create mode 100644 test/test_gdb_drivers.py diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 9c50b00..0fceb8c 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -8,6 +8,16 @@ bitstrings used *as sampled*, with no product completion. For TPB and the SQD loops, see [`../tpb/README.md`](../tpb/README.md). +## The default data + +Both drivers default `--fcidump` and `--detfiles` to upstream's own GDB app data under +`vendor/sbd-upstream/apps/chemistry_gdb_selected_basis_diagonalization/`: the Fe4S4 +FCIDUMP at 36 orbitals, and its four `det0.txt`-`det3.txt` holding 14,884 determinants +each, 59,536 in total. **Any command below that passes neither flag runs exactly that +case** — the determinants come from those four files, not from nowhere. Upstream +publishes no reference energy for it, so treat it as a self-consistency benchmark +rather than a validation target. + ## run_gdb_diag.py — standalone GDB diagonalization Stays on the in-memory entry point (`sbd.gdb_diag`) throughout: determinant text is @@ -15,10 +25,16 @@ read in Python and handed to the binding as a list. SBD's own file-based entry point is deliberately not used. ```bash -# Fe4S4, upstream's own GDB data: 4 files, 59,536 determinants, 36 orbitals. -# Upstream publishes no reference energy for this case. +# The default case: Fe4S4, upstream's four det files, 59,536 determinants. python run_gdb_diag.py +# The same run with the defaults written out -- this is the --detfiles syntax to +# copy for your own data. Files are comma-separated and concatenated in Python +# into one in-memory list, so their combined order must be sorted and disjoint. +GDB=../../vendor/sbd-upstream/apps/chemistry_gdb_selected_basis_diagonalization +python run_gdb_diag.py --fcidump $GDB/fcidump_Fe4S4.txt \ + --detfiles $GDB/det0.txt,$GDB/det1.txt,$GDB/det2.txt,$GDB/det3.txt + # Shard the basis: Fe4S4's four files, one per rank. For this data that is # already the balanced globally-sorted split, so no redistribution is needed. mpirun -np 4 -x OMP_NUM_THREADS=8 python run_gdb_diag.py --b_comm_size 4 @@ -41,6 +57,13 @@ python run_gdb_diag.py \ # Grow the subspace with SBD's own heatbath expansion (one round; the expanded # list comes back as carryover_det) python run_gdb_diag.py --carryover_type 2 --heatbath_cutoff 1e-4 + +# On GPUs, use --device gpu (Thrust). It requires helper == 1, so every rank has +# to go to t*b -- which means --b_comm_size is not optional for multi-GPU GDB. +mpirun -np 4 python run_gdb_diag.py --device gpu --b_comm_size 4 + +# One rank, one GPU is the exception: b=1 already leaves helper == 1. +python run_gdb_diag.py --device gpu ``` The determinant bit order is the one the upstream app documents: reading from the @@ -55,7 +78,7 @@ no extra machinery — `carryover_type` 2 and 3 return the parents *together wit new candidates, so one round's result **is** the next round's subspace. ```bash -# Fe4S4 from upstream's shipped subspace, one cutoff +# Fe4S4 from upstream's shipped subspace (the default data above), one cutoff python run_gdb_heatbath.py --cutoffs 1e-3 # A ladder, stopping before the subspace passes 2M determinants @@ -73,6 +96,44 @@ mpirun -np 4 python run_gdb_heatbath.py --b_comm_size 4 --device gpu --cutoffs 1 python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log ladder.json ``` +### Parameters + +Seed — where the starting subspace comes from: + +| Parameter | What it controls | Default | +|---|---|---| +| `--seed` | `files` reads `--detfiles`; `hf` starts from the single Hartree-Fock determinant; `from-alpha` builds the `\|A\|^2` product of an alpha list; `strings` reads full determinants, e.g. sampled configurations | `files` | +| `--fcidump` | FCIDUMP defining the Hamiltonian | Fe4S4, see [The default data](#the-default-data) | +| `--detfiles` | `--seed files`: comma-separated determinant files, concatenated in Python. Their combined order must be sorted and disjoint | upstream's four Fe4S4 files | +| `--alpha-file` / `--alpha-limit` | `--seed from-alpha`: the alpha list, and a cap on how many of its strings to keep. The product costs `N^2` determinants, so this is the size dial | none / `0` (all) | +| `--strings-file` | `--seed strings`: a file of `2*norb`-bit determinants | none | + +Ladder — how far the expansion is pushed. A rung runs rounds at one cutoff until the +energy stops moving, then the next cutoff begins: + +| Parameter | What it controls | Default | +|---|---|---| +| `--cutoffs` | The ladder itself: comma-separated `heatbath_cutoff` values, smallest step last. Each rung admits candidates whose estimated contribution exceeds it | `1e-3` | +| `--max_rounds` | Cap on rounds **per rung**, so a rung that never converges cannot run forever | `8` | +| `--energy_tol` | Advance to the next rung once `\|dE\|` between rounds falls below this (Hartree) | `1e-5` | +| `--max_dim` | Stop before diagonalizing a subspace larger than this. The brake that keeps rows of a comparison cost-matched | `0` (no cap) | +| `--carryover_type` | Heatbath variant, 2 or 3. Types 0 and 1 do not expand, so they cannot drive the loop | `2` | +| `--heatbath_truncation` | Weight threshold applied to **parents**, before expanding — not the cutoff. See the warning below; leave it at 0 | `0.0` | +| `--heatbath_batch_size` | Expansion batch size per rank | `1000000` | + +Solver, MPI and output — the same meanings as in `run_gdb_diag.py`: + +| Parameter | What it controls | Default | +|---|---|---| +| `--device` | `cpu`, `gpu` (Thrust, the only GPU backend with GDB kernels), `gpu-omp` (no GDB kernels — runs on the host), `auto` | `cpu` | +| `--method` | 0=Davidson, 1=Davidson storing the Hamiltonian. GDB has no Lanczos | `0` | +| `--tolerance` / `--iteration` / `--block` | Davidson residual tolerance, iteration cap, and basis-vector count. Also accepted as `--eps` / `--max_it` / `--max_nb`, matching SBD's own names | `1e-6` / `30` / `10` | +| `--bit_length` | Bits per packed word | `64` | +| `--b_comm_size` | Basis shards — the only dimension that divides memory | `1` | +| `--t_comm_size` | Tasks per ring station; must not exceed `--b_comm_size` | `1` | +| `--determinant_distribution` | Placement across `b_comm`; see [Placement across `b_comm`](#placement-across-b_comm) | `equal-bra-a` | +| `--log FILE` | Write the per-round `(rung, cutoff, round, dimension, energy, delta_energy, seconds)` series as JSON — the machine-readable form of the table the driver prints | none | + ### Why a cutoff *ladder* `heatbath_cutoff` admits a candidate when its estimated contribution exceeds the diff --git a/test/test_gdb_drivers.py b/test/test_gdb_drivers.py new file mode 100644 index 0000000..cef5fad --- /dev/null +++ b/test/test_gdb_drivers.py @@ -0,0 +1,232 @@ +# 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. + +"""Check the GDB example drivers end to end, as a user invokes them. + +``test_gdb_equivalence.py`` covers the binding: it calls ``sbd.gdb_diag`` +directly and proves the solver right. Nothing covered the layer a user actually +touches -- the argument parsing, the determinant text reading, the sharding +arithmetic, the default file paths -- so a driver could break while every +library test stayed green. + +These run the scripts in a **subprocess**, the same device the Fe4S4 case in +``test_gdb_equivalence.py`` uses and for the same reason: ``det_vector``'s row +width is fixed process-wide on first use, so h2o (one word) and Fe4S4 (two) +cannot share a process. Invoking the driver as a subprocess also means the test +exercises ``__main__``, argument defaults included, rather than importing a +function and bypassing them. + +The energies asserted here are not new reference values. The small h2o case +reproduces the same -76.0588897208 that ``test_gdb_equivalence.py`` anchors +against TPB, which ties the driver to the solver: if they ever disagree, the +driver's reading or sharding is at fault, not the diagonalization. +""" + +from __future__ import annotations + +import json +import pathlib +import re +import subprocess +import sys + +import pytest + +REPO_ROOT = pathlib.Path(__file__).resolve().parents[1] +DRIVER_DIR = REPO_ROOT / "examples" / "gdb" +H2O_DIR = REPO_ROOT / "vendor" / "sbd-upstream" / "data" / "h2o" +FE4S4_DIR = ( + REPO_ROOT + / "vendor" + / "sbd-upstream" + / "apps" + / "chemistry_gdb_selected_basis_diagonalization" +) + +# The 24x24 product of the h2o alpha list, interleaved into 576 full +# determinants. Asserted through the binding in test_gdb_equivalence.py; the +# driver must land on the same value from the same input files. +H2O_576_ENERGY = -76.0588897208 + +# Upstream's shipped GDB subspace: four files of 14,884 determinants over 36 +# orbitals. This is what both drivers use when neither --fcidump nor --detfiles +# is given, so a test of the defaults is also a test that those paths resolve. +FE4S4_DIM = 59536 +FE4S4_ENERGY = -326.6982518821 + +# Small enough to keep the guardrail cases near-instant: 12 alphas is a +# 144-determinant subspace, and these never reach a diagonalization anyway. +TINY_ALPHA_LIMIT = 12 + + +def _h2o_args(alpha_limit: int) -> list[str]: + """Point a driver at the bundled h2o data as an |A|^2 product subspace.""" + return [ + "--fcidump", str(H2O_DIR / "fcidump.txt"), + "--from-alpha", str(H2O_DIR / "h2o-1em3-alpha.txt"), + "--alpha-limit", str(alpha_limit), + ] + + +def _run(script: str, args: list[str], timeout: int = 900): + """Run a driver from its own directory, the way the README shows it.""" + if not (DRIVER_DIR / script).is_file(): + pytest.skip(f"driver not found: {DRIVER_DIR / script}") + return subprocess.run( + [sys.executable, script, *args], + cwd=DRIVER_DIR, + capture_output=True, + text=True, + timeout=timeout, + ) + + +def _require(path: pathlib.Path) -> None: + if not path.exists(): + pytest.skip(f"vendored upstream data not found at {path} (submodule not checked out?)") + + +def _assert_ok(completed) -> str: + assert completed.returncode == 0, ( + "driver exited non-zero:\n" + f"--- stdout ---\n{completed.stdout[-3000:]}\n" + f"--- stderr ---\n{completed.stderr[-3000:]}" + ) + return completed.stdout + + +def _energy(stdout: str) -> float: + match = re.search(r"Ground state energy:\s*(-?\d+\.\d+)", stdout) + assert match, f"no energy line in driver output:\n{stdout[-3000:]}" + return float(match.group(1)) + + +def _dimension(stdout: str) -> int: + match = re.search(r"Subspace dimension:\s*(\d+)", stdout) + assert match, f"no dimension line in driver output:\n{stdout[-3000:]}" + return int(match.group(1)) + + +# --- run_gdb_diag.py ---------------------------------------------------------- + + +def test_diag_reproduces_the_library_anchor(): + """The driver on h2o must equal what the binding gives on the same subspace.""" + _require(H2O_DIR / "h2o-1em3-alpha.txt") + stdout = _assert_ok(_run("run_gdb_diag.py", _h2o_args(24))) + assert _dimension(stdout) == 576 + assert _energy(stdout) == pytest.approx(H2O_576_ENERGY, abs=1e-8) + + +def test_diag_reports_a_sane_electron_count(): + """The occupation density must sum to the electron count, not merely exist. + + A subspace built with the wrong bit order still diagonalizes and still + returns a plausible-looking energy; the electron count is what catches it, + which is why the driver prints the sum rather than the vector alone. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + stdout = _assert_ok(_run("run_gdb_diag.py", _h2o_args(24))) + match = re.search(r"sums to (\d+\.\d+); should equal the electron count (\d+)", stdout) + assert match, f"no electron-count check in output:\n{stdout[-3000:]}" + assert float(match.group(1)) == pytest.approx(float(match.group(2)), abs=1e-6) + + +def test_diag_carryover_returns_parents_with_the_candidates(): + """``carryover_type 2`` must come back larger than the subspace it expanded. + + The expansion starts from ``edet = det`` (gdb/expansion.h:543), so the + returned list is the next subspace rather than only the additions -- the + property the heatbath driver's loop depends on. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + stdout = _assert_ok(_run("run_gdb_diag.py", [ + *_h2o_args(TINY_ALPHA_LIMIT), "--carryover_type", "2", + "--heatbath_cutoff", "1e-3", + ])) + match = re.search(r"Carryover determinants on this rank:\s*(\d+)", stdout) + assert match, f"no carryover line in output:\n{stdout[-3000:]}" + assert int(match.group(1)) > _dimension(stdout) + + +@pytest.mark.parametrize( + "extra, expected", + [ + (["--t_comm_size", "2", "--b_comm_size", "1"], "t_comm_size <= b_comm_size"), + (["--b_comm_size", "2"], "divide the rank count"), + ], + ids=["t-exceeds-b", "decomposition-does-not-divide-ranks"], +) +def test_diag_rejects_an_impossible_decomposition(extra, expected): + """Both guardrails must fail loudly on one rank rather than segfault. + + ``t > b`` is the dangerous one: upstream does not check it, and MakeHelpers + then dereferences an empty lookup, so without this guard the failure mode is + a segfault rather than a message. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + completed = _run("run_gdb_diag.py", [*_h2o_args(TINY_ALPHA_LIMIT), *extra]) + assert completed.returncode != 0, f"expected a refusal, got:\n{completed.stdout[-2000:]}" + combined = completed.stdout + completed.stderr + assert expected in combined, f"guardrail message missing:\n{combined[-3000:]}" + + +@pytest.mark.slow +def test_diag_default_files_are_the_fe4s4_subspace(): + """With no input flags at all, the driver must run upstream's Fe4S4 data. + + This is the only test of the default paths, and the README documents them as + the case every flagless command runs -- so if the vendored layout moves, this + is what says so. + """ + _require(FE4S4_DIR / "det0.txt") + stdout = _assert_ok(_run("run_gdb_diag.py", [], timeout=1800)) + assert _dimension(stdout) == FE4S4_DIM + assert _energy(stdout) == pytest.approx(FE4S4_ENERGY, abs=1e-8) + + +# --- run_gdb_heatbath.py ------------------------------------------------------ + + +def test_heatbath_ladder_grows_and_lowers_the_energy(tmp_path): + """Two rounds on h2o: the subspace must grow and the energy must not rise. + + The energy is variational in the subspace, and each round's subspace + contains the previous one, so a rise means the expansion dropped parents or + a round failed to converge. Round 0 diagonalizes the seed untouched, so it + must equal the same anchor the diagonalization driver reports. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + log = tmp_path / "ladder.json" + _assert_ok(_run("run_gdb_heatbath.py", [ + "--fcidump", str(H2O_DIR / "fcidump.txt"), + "--seed", "from-alpha", + "--alpha-file", str(H2O_DIR / "h2o-1em3-alpha.txt"), + "--alpha-limit", "24", + "--cutoffs", "1e-3", + "--max_rounds", "2", + "--log", str(log), + ])) + + rounds = json.loads(log.read_text())["rounds"] + assert len(rounds) >= 2, f"expected at least two rounds, got {rounds}" + + assert rounds[0]["dimension"] == 576 + assert rounds[0]["energy"] == pytest.approx(H2O_576_ENERGY, abs=1e-8) + assert rounds[0]["delta_energy"] is None, "the seed round has nothing to compare against" + + dims = [entry["dimension"] for entry in rounds] + energies = [entry["energy"] for entry in rounds] + assert dims == sorted(dims), f"subspace shrank across rounds: {dims}" + assert dims[-1] > dims[0], f"expansion added nothing: {dims}" + for before, after in zip(energies, energies[1:]): + assert after <= before + 1e-10, f"energy rose across a round: {energies}" From a37b6faaa185b9be182621c1f51f2bad4e55b98e Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 17:08:24 -0400 Subject: [PATCH 16/22] examples/gdb/README: list the two backends that work, and say which input flags differ The backend block listed all four --device values as if they were alternatives. For GDB only two are: gpu-omp has no GDB kernels and runs on the host, and auto can resolve to gpu-omp and do the same. List cpu and gpu, then one sentence on why the other two are accepted but not GPU paths. Also state plainly that --fcidump and --detfiles are identical across the two drivers, and that two input flags are not: the alpha list is --alpha-file in the heatbath driver but --from-alpha in the diagonalization one, and --seed selects the subspace source here while there it is an integer RNG seed. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 0fceb8c..430da4a 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -98,6 +98,12 @@ python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log ladder.json ### Parameters +`--fcidump` and `--detfiles` are spelled and defaulted exactly as in `run_gdb_diag.py`, +so a command that feeds one driver its data feeds the other. Two input flags do **not** +carry over: the alpha list is `--alpha-file` here but `--from-alpha` there, and `--seed` +means different things in the two drivers — here it selects where the subspace comes +from, while in `run_gdb_diag.py` it is the integer RNG seed for a random initial vector. + Seed — where the starting subspace comes from: | Parameter | What it controls | Default | @@ -125,7 +131,7 @@ Solver, MPI and output — the same meanings as in `run_gdb_diag.py`: | Parameter | What it controls | Default | |---|---|---| -| `--device` | `cpu`, `gpu` (Thrust, the only GPU backend with GDB kernels), `gpu-omp` (no GDB kernels — runs on the host), `auto` | `cpu` | +| `--device` | `cpu` or `gpu` (Thrust) are the two real choices; `gpu-omp` and `auto` run GDB on the host, see [Choosing a backend](#choosing-a-backend) | `cpu` | | `--method` | 0=Davidson, 1=Davidson storing the Hamiltonian. GDB has no Lanczos | `0` | | `--tolerance` / `--iteration` / `--block` | Davidson residual tolerance, iteration cap, and basis-vector count. Also accepted as `--eps` / `--max_it` / `--max_nb`, matching SBD's own names | `1e-6` / `30` / `10` | | `--bit_length` | Bits per packed word | `64` | @@ -355,10 +361,13 @@ already count-balanced. ```bash --device cpu # host OpenMP (default) --device gpu # NVHPC Thrust, NVIDIA only -- the only GPU backend with GDB kernels ---device gpu-omp # OpenMP target offload; compiles for GDB but has no GDB kernels ---device auto # GPU if one is available, else CPU ``` +`gpu-omp` and `auto` are accepted too, but for GDB neither is a GPU path: `gpu-omp` has +no GDB kernels and silently diagonalizes on the host (the drivers warn), and `auto` +prefers Thrust but falls back to `gpu-omp` where Thrust is not built — landing on the +host as well. **For GDB, pick `cpu` or `gpu` explicitly.** + `sbd.available_backends()` reports what this install actually built — a static scan, so it is safe to call outside `mpirun` — and `sbd.loaded_backends()` reports what the current process has imported. From 141afd797017f7e10d5ab3a27e2fec3ed63c0d72 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 17:18:13 -0400 Subject: [PATCH 17/22] examples/gdb: rename --seed to --subspace-from, and accept either alpha-file spelling The two drivers had grown a flag collision and a gratuitous difference. --seed took a string here (files|hf|from-alpha|strings) and an int in run_gdb_diag.py, where it is the RNG seed for a random initial vector. Same flag, incompatible types, so a command could not be moved between the drivers and `--seed hf` failed confusingly in one of them. Renamed to --subspace-from, which says what it does. --seed still works and prints a deprecation notice to stderr; argparse's own deprecated= needs 3.13 and this package supports 3.10, so the old spelling is detected in argv instead. The alpha determinant list was --from-alpha in one driver and --alpha-file in the other, for no reason. Both drivers now accept both spellings. Tests cover all three paths: either alpha spelling gives the same 576-determinant energy, and the deprecated --seed still runs while emitting the notice. 9 passed / 1 slow-skipped in the driver file, 26 passed in the serial suite, 13 under mpirun -n 2. README parameter tables re-checked against argparse: no invented flags, none undocumented. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 25 ++++++++--------- examples/gdb/run_gdb_diag.py | 2 +- examples/gdb/run_gdb_heatbath.py | 36 +++++++++++++++++-------- test/test_gdb_drivers.py | 46 +++++++++++++++++++++++++++++++- 4 files changed, 84 insertions(+), 25 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 430da4a..25bf3c9 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -85,7 +85,7 @@ python run_gdb_heatbath.py --cutoffs 1e-3 python run_gdb_heatbath.py --cutoffs 1e-3,1e-4,1e-5 --max_dim 2000000 # The no-input null: start from the Hartree-Fock determinant alone -python run_gdb_heatbath.py --seed hf --cutoffs 1e-3,1e-4 +python run_gdb_heatbath.py --subspace-from hf --cutoffs 1e-3,1e-4 # Sharded. At --b_comm_size == ranks with t=1 the expansion comes back already # sharded for the next round, so nothing has to be gathered. @@ -98,21 +98,22 @@ python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log ladder.json ### Parameters -`--fcidump` and `--detfiles` are spelled and defaulted exactly as in `run_gdb_diag.py`, -so a command that feeds one driver its data feeds the other. Two input flags do **not** -carry over: the alpha list is `--alpha-file` here but `--from-alpha` there, and `--seed` -means different things in the two drivers — here it selects where the subspace comes -from, while in `run_gdb_diag.py` it is the integer RNG seed for a random initial vector. +Input flags match `run_gdb_diag.py`: `--fcidump`, `--detfiles` and `--alpha-limit` are +spelled and defaulted identically, and the alpha list is accepted as either +`--alpha-file` or `--from-alpha` by both drivers, so a command that feeds one its data +feeds the other. The one flag that is *not* shared is `--seed`, which in +`run_gdb_diag.py` is the integer RNG seed for a random initial vector; this driver's +equivalent is `--subspace-from` (`--seed` still works here but warns, and will go). Seed — where the starting subspace comes from: | Parameter | What it controls | Default | |---|---|---| -| `--seed` | `files` reads `--detfiles`; `hf` starts from the single Hartree-Fock determinant; `from-alpha` builds the `\|A\|^2` product of an alpha list; `strings` reads full determinants, e.g. sampled configurations | `files` | +| `--subspace-from` | `files` reads `--detfiles`; `hf` starts from the single Hartree-Fock determinant; `from-alpha` builds the `\|A\|^2` product of an alpha list; `strings` reads full determinants, e.g. sampled configurations | `files` | | `--fcidump` | FCIDUMP defining the Hamiltonian | Fe4S4, see [The default data](#the-default-data) | -| `--detfiles` | `--seed files`: comma-separated determinant files, concatenated in Python. Their combined order must be sorted and disjoint | upstream's four Fe4S4 files | -| `--alpha-file` / `--alpha-limit` | `--seed from-alpha`: the alpha list, and a cap on how many of its strings to keep. The product costs `N^2` determinants, so this is the size dial | none / `0` (all) | -| `--strings-file` | `--seed strings`: a file of `2*norb`-bit determinants | none | +| `--detfiles` | `--subspace-from files`: comma-separated determinant files, concatenated in Python. Their combined order must be sorted and disjoint | upstream's four Fe4S4 files | +| `--alpha-file` (or `--from-alpha`) / `--alpha-limit` | `--subspace-from from-alpha`: the alpha list, and a cap on how many of its strings to keep. The product costs `N^2` determinants, so this is the size dial | none / `0` (all) | +| `--strings-file` | `--subspace-from strings`: a file of `2*norb`-bit determinants | none | Ladder — how far the expansion is pushed. A rung runs rounds at one cutoff until the energy stops moving, then the next cutoff begins: @@ -172,10 +173,10 @@ run can never cross, which makes a sparse run self-checking even without a refer ### Seeds -`--seed` selects where the starting subspace comes from, so one driver produces +`--subspace-from` selects where the starting subspace comes from, so one driver produces every row of a seed comparison: -| `--seed` | starting subspace | +| `--subspace-from` | starting subspace | |---|---| | `files` *(default)* | determinant files, defaulting to upstream's four Fe4S4 files | | `hf` | the Hartree-Fock determinant alone, built from the FCIDUMP header — the no-input null | diff --git a/examples/gdb/run_gdb_diag.py b/examples/gdb/run_gdb_diag.py index eaa83af..89b6fe3 100755 --- a/examples/gdb/run_gdb_diag.py +++ b/examples/gdb/run_gdb_diag.py @@ -96,7 +96,7 @@ def parse_args(): help='Comma-separated files of full determinants, one ' '2*norb-bit string per line. Concatenated in Python ' 'and passed as a single in-memory list') - parser.add_argument('--from-alpha', default='', metavar='FILE', + parser.add_argument('--from-alpha', '--alpha-file', default='', metavar='FILE', dest='from_alpha', help='Instead of --detfiles, read a norb-bit alpha list and ' 'form the full |A|^2 product basis by interleaving it ' diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py index 3b4cc5f..bce60fc 100755 --- a/examples/gdb/run_gdb_heatbath.py +++ b/examples/gdb/run_gdb_heatbath.py @@ -88,20 +88,24 @@ def parse_args(): help='FCIDUMP file defining the Hamiltonian') # --- the seed --------------------------------------------------------- - parser.add_argument('--seed', default='files', + # --seed is the old spelling, kept working but deprecated: run_gdb_diag.py uses + # --seed for the integer RNG seed of its random initial vector, so the same flag + # meant two incompatible things across the two drivers. + parser.add_argument('--subspace-from', '--seed', default='files', + dest='subspace_from', choices=['files', 'hf', 'from-alpha', 'strings'], help='Where the starting subspace comes from') parser.add_argument('--detfiles', default=_DEFAULT_DETFILES, - help='--seed files: comma-separated files of 2*norb-bit ' + help='--subspace-from files: comma-separated files of 2*norb-bit ' 'determinant strings') - parser.add_argument('--alpha-file', default='', dest='alpha_file', - help='--seed from-alpha: a norb-bit alpha determinant list, ' + parser.add_argument('--alpha-file', '--from-alpha', default='', dest='alpha_file', + help='--subspace-from from-alpha: a norb-bit alpha determinant list, ' 'interleaved with itself into the full product basis') parser.add_argument('--alpha-limit', type=int, default=0, dest='alpha_limit', - help='--seed from-alpha: keep only the first N alpha strings ' + help='--subspace-from from-alpha: keep only the first N alpha strings ' '(the product costs N^2 determinants)') parser.add_argument('--strings-file', default='', dest='strings_file', - help='--seed strings: a file of 2*norb-bit determinant ' + help='--subspace-from strings: a file of 2*norb-bit determinant ' 'strings, e.g. sampled configurations') # --- the ladder ------------------------------------------------------- @@ -154,7 +158,17 @@ def parse_args(): parser.add_argument('--log', default='', metavar='FILE', help='Write the per-round (dimension, energy) series as JSON') - return parser.parse_args() + args = parser.parse_args() + + # argparse's own `deprecated=` needs Python 3.13, and this supports 3.10, so + # detect the old spelling in argv instead. Warn rather than fail: --seed still + # works, it is only ambiguous across the two drivers. + if any(a == '--seed' or a.startswith('--seed=') for a in sys.argv[1:]): + print("NOTE: --seed is deprecated here and will be removed; use " + "--subspace-from. run_gdb_diag.py uses --seed for its integer RNG " + "seed, so the same flag meant two different things.", file=sys.stderr) + + return args def read_strings(paths): @@ -210,10 +224,10 @@ def build_seed(args, sbd, norb, nelec, ms2, rank): """The starting determinant array for this rank, plus a description.""" total_bits = 2 * norb - if args.seed == 'hf': + if args.subspace_from == 'hf': strings = [hartree_fock_string(norb, nelec, ms2)] source = f"Hartree-Fock determinant ({nelec} electrons, MS2={ms2})" - elif args.seed == 'from-alpha': + elif args.subspace_from == 'from-alpha': if not args.alpha_file: raise ValueError("--seed from-alpha requires --alpha-file") alpha = read_strings([args.alpha_file]) @@ -226,10 +240,10 @@ def build_seed(args, sbd, norb, nelec, ms2, rank): strings = [interleave(a, b) for a in alpha for b in alpha] source = f"{args.alpha_file}: {len(alpha)} alpha -> {len(alpha)}^2 product" else: - paths = ([args.strings_file] if args.seed == 'strings' + paths = ([args.strings_file] if args.subspace_from == 'strings' else [p for p in args.detfiles.split(',') if p]) if not paths or not all(paths): - raise ValueError(f"--seed {args.seed} requires input file(s)") + raise ValueError(f"--subspace-from {args.subspace_from} requires input file(s)") strings = read_strings(paths) source = f"{len(paths)} file(s): {', '.join(paths)}" diff --git a/test/test_gdb_drivers.py b/test/test_gdb_drivers.py index cef5fad..605b0c9 100644 --- a/test/test_gdb_drivers.py +++ b/test/test_gdb_drivers.py @@ -209,7 +209,7 @@ def test_heatbath_ladder_grows_and_lowers_the_energy(tmp_path): log = tmp_path / "ladder.json" _assert_ok(_run("run_gdb_heatbath.py", [ "--fcidump", str(H2O_DIR / "fcidump.txt"), - "--seed", "from-alpha", + "--subspace-from", "from-alpha", "--alpha-file", str(H2O_DIR / "h2o-1em3-alpha.txt"), "--alpha-limit", "24", "--cutoffs", "1e-3", @@ -230,3 +230,47 @@ def test_heatbath_ladder_grows_and_lowers_the_energy(tmp_path): assert dims[-1] > dims[0], f"expansion added nothing: {dims}" for before, after in zip(energies, energies[1:]): assert after <= before + 1e-10, f"energy rose across a round: {energies}" + + +@pytest.mark.parametrize( + "alpha_flag", ["--from-alpha", "--alpha-file"], ids=["from-alpha", "alpha-file"] +) +def test_both_drivers_accept_either_alpha_flag_spelling(alpha_flag): + """The two drivers grew different names for the same input; both now take both. + + Asserted on the diagonalization driver, where ``--from-alpha`` is the primary + name and ``--alpha-file`` the alias added for parity with the heatbath driver. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + stdout = _assert_ok(_run("run_gdb_diag.py", [ + "--fcidump", str(H2O_DIR / "fcidump.txt"), + alpha_flag, str(H2O_DIR / "h2o-1em3-alpha.txt"), + "--alpha-limit", "24", + ])) + assert _dimension(stdout) == 576 + assert _energy(stdout) == pytest.approx(H2O_576_ENERGY, abs=1e-8) + + +def test_heatbath_still_accepts_the_deprecated_seed_flag(tmp_path): + """``--seed`` keeps working, and says it is on the way out. + + It is the old spelling of ``--subspace-from``, renamed because + ``run_gdb_diag.py`` uses ``--seed`` for an integer RNG seed -- the same flag + taking a string here and an int there. Kept working so existing commands do + not break; warned about so they get updated. + """ + _require(H2O_DIR / "h2o-1em3-alpha.txt") + log = tmp_path / "ladder.json" + completed = _run("run_gdb_heatbath.py", [ + "--fcidump", str(H2O_DIR / "fcidump.txt"), + "--seed", "from-alpha", + "--alpha-file", str(H2O_DIR / "h2o-1em3-alpha.txt"), + "--alpha-limit", str(TINY_ALPHA_LIMIT), + "--cutoffs", "1e-3", "--max_rounds", "1", + "--log", str(log), + ]) + _assert_ok(completed) + assert "--seed is deprecated" in completed.stderr, ( + f"no deprecation notice:\n{completed.stderr[-2000:]}" + ) + assert json.loads(log.read_text())["rounds"][0]["dimension"] == TINY_ALPHA_LIMIT ** 2 From 78cab43490504595f5240e90a45cc7f5ff936f69 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 17:23:19 -0400 Subject: [PATCH 18/22] examples/gdb: drop the --seed alias, and fix the mentions the rename missed The branch is unmerged, so nothing depends on the old spelling: --seed is gone rather than deprecated, along with its argv-sniffing warning and the test for it. `--seed hf` now exits 2 with "unrecognized arguments". The previous commit's rename was pattern-by-pattern and missed four lines of the module docstring plus two messages, so `--help` still advertised --seed while argparse had moved on. Swept the whole file this time and asserted no mention survives; the only remaining --seed in examples/gdb is the README line contrasting it with run_gdb_diag.py's integer RNG seed, which is deliberate. 9 passed with --run-slow, 25 in the serial suite, 13 under mpirun -n 2. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 6 +++--- examples/gdb/run_gdb_heatbath.py | 31 +++++++++---------------------- test/test_gdb_drivers.py | 25 ------------------------- 3 files changed, 12 insertions(+), 50 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 25bf3c9..a50fc58 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -101,9 +101,9 @@ python run_gdb_heatbath.py --cutoffs 1e-3,1e-4 --log ladder.json Input flags match `run_gdb_diag.py`: `--fcidump`, `--detfiles` and `--alpha-limit` are spelled and defaulted identically, and the alpha list is accepted as either `--alpha-file` or `--from-alpha` by both drivers, so a command that feeds one its data -feeds the other. The one flag that is *not* shared is `--seed`, which in -`run_gdb_diag.py` is the integer RNG seed for a random initial vector; this driver's -equivalent is `--subspace-from` (`--seed` still works here but warns, and will go). +feeds the other. The one flag that is *not* shared is `--seed`: in +`run_gdb_diag.py` it is the integer RNG seed for a random initial vector, and this +driver has none — `--subspace-from` selects the starting subspace and is unrelated. Seed — where the starting subspace comes from: diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py index bce60fc..3d771d2 100755 --- a/examples/gdb/run_gdb_heatbath.py +++ b/examples/gdb/run_gdb_heatbath.py @@ -35,10 +35,10 @@ Seeds, so the same driver produces every row of a seed comparison: - --seed files determinant files (default: upstream's four Fe4S4 files) - --seed hf the Hartree-Fock determinant alone, the no-input null - --seed from-alpha an alpha list interleaved with itself (a TPB-shaped space) - --seed strings an arbitrary bitstring file, e.g. sampled configurations + --subspace-from files determinant files (default: upstream's four Fe4S4 files) + --subspace-from hf the Hartree-Fock determinant alone, the no-input null + --subspace-from from-alpha an alpha list interleaved with itself (a TPB-shaped space) + --subspace-from strings an arbitrary bitstring file, e.g. sampled configurations Usage: # Fe4S4 from upstream's shipped subspace, one cutoff @@ -48,7 +48,7 @@ python run_gdb_heatbath.py --cutoffs 1e-3,1e-4,1e-5 --max_dim 2000000 # The null hypothesis: no input subspace at all, just Hartree-Fock - python run_gdb_heatbath.py --seed hf --cutoffs 1e-3,1e-4 + python run_gdb_heatbath.py --subspace-from hf --cutoffs 1e-3,1e-4 # Sharded across 4 ranks/GPUs. At --b_comm_size == ranks with t=1 the # expansion comes back already sharded for the next round, no gather needed. @@ -88,10 +88,7 @@ def parse_args(): help='FCIDUMP file defining the Hamiltonian') # --- the seed --------------------------------------------------------- - # --seed is the old spelling, kept working but deprecated: run_gdb_diag.py uses - # --seed for the integer RNG seed of its random initial vector, so the same flag - # meant two incompatible things across the two drivers. - parser.add_argument('--subspace-from', '--seed', default='files', + parser.add_argument('--subspace-from', default='files', dest='subspace_from', choices=['files', 'hf', 'from-alpha', 'strings'], help='Where the starting subspace comes from') @@ -158,17 +155,7 @@ def parse_args(): parser.add_argument('--log', default='', metavar='FILE', help='Write the per-round (dimension, energy) series as JSON') - args = parser.parse_args() - - # argparse's own `deprecated=` needs Python 3.13, and this supports 3.10, so - # detect the old spelling in argv instead. Warn rather than fail: --seed still - # works, it is only ambiguous across the two drivers. - if any(a == '--seed' or a.startswith('--seed=') for a in sys.argv[1:]): - print("NOTE: --seed is deprecated here and will be removed; use " - "--subspace-from. run_gdb_diag.py uses --seed for its integer RNG " - "seed, so the same flag meant two different things.", file=sys.stderr) - - return args + return parser.parse_args() def read_strings(paths): @@ -194,7 +181,7 @@ def interleave(alpha, beta): def hartree_fock_string(norb, nelec, ms2): """The Hartree-Fock determinant: the lowest orbitals doubly occupied. - Built from the FCIDUMP header alone, so ``--seed hf`` needs no input subspace at + Built from the FCIDUMP header alone, so ``--subspace-from hf`` needs no input subspace at all. That is the point of it: if expansion from HF reaches the same place as expansion from a sampled subspace, the sampling added nothing. """ @@ -229,7 +216,7 @@ def build_seed(args, sbd, norb, nelec, ms2, rank): source = f"Hartree-Fock determinant ({nelec} electrons, MS2={ms2})" elif args.subspace_from == 'from-alpha': if not args.alpha_file: - raise ValueError("--seed from-alpha requires --alpha-file") + raise ValueError("--subspace-from from-alpha requires --alpha-file") alpha = read_strings([args.alpha_file]) if args.alpha_limit: alpha = alpha[:args.alpha_limit] diff --git a/test/test_gdb_drivers.py b/test/test_gdb_drivers.py index 605b0c9..b7dedbd 100644 --- a/test/test_gdb_drivers.py +++ b/test/test_gdb_drivers.py @@ -249,28 +249,3 @@ def test_both_drivers_accept_either_alpha_flag_spelling(alpha_flag): ])) assert _dimension(stdout) == 576 assert _energy(stdout) == pytest.approx(H2O_576_ENERGY, abs=1e-8) - - -def test_heatbath_still_accepts_the_deprecated_seed_flag(tmp_path): - """``--seed`` keeps working, and says it is on the way out. - - It is the old spelling of ``--subspace-from``, renamed because - ``run_gdb_diag.py`` uses ``--seed`` for an integer RNG seed -- the same flag - taking a string here and an int there. Kept working so existing commands do - not break; warned about so they get updated. - """ - _require(H2O_DIR / "h2o-1em3-alpha.txt") - log = tmp_path / "ladder.json" - completed = _run("run_gdb_heatbath.py", [ - "--fcidump", str(H2O_DIR / "fcidump.txt"), - "--seed", "from-alpha", - "--alpha-file", str(H2O_DIR / "h2o-1em3-alpha.txt"), - "--alpha-limit", str(TINY_ALPHA_LIMIT), - "--cutoffs", "1e-3", "--max_rounds", "1", - "--log", str(log), - ]) - _assert_ok(completed) - assert "--seed is deprecated" in completed.stderr, ( - f"no deprecation notice:\n{completed.stderr[-2000:]}" - ) - assert json.loads(log.read_text())["rounds"][0]["dimension"] == TINY_ALPHA_LIMIT ** 2 From cb248bb2fdde47e76af0d7a1f9b50a03e7aa4b99 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 17:30:14 -0400 Subject: [PATCH 19/22] examples/gdb/README: how to feed a qiskit-addon-sqd counts file to --subspace-from --strings-file takes one determinant per line, not a counts JSON, and the conversion has a step that is easy to miss: qiskit-addon-sqd emits [beta | alpha] concatenated while GDB wants the bits interleaved. Document both steps with a snippet built on the driver's own interleave(), and note that a set/sorted() handles the duplicate rejection and ordering the shard contract expects. Skipping the interleave is worth spelling out because it does not fail cleanly: the concatenated strings diagonalize to a plausible number (-66.8043 where the correct conversion gives -76.0724) and only abort later inside the heatbath expansion with std::out_of_range. The electron-count check cannot catch it -- permuting bits preserves how many are set -- so the density still sums to 10. Snippet and command both run as written; the snippet reproduces the same file byte-for-byte as the conversion that was tested. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 51 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index a50fc58..c0d0041 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -197,6 +197,57 @@ different sizes from the same cutoff, so a cutoff-matched comparison mostly repo subspace size rather than seed quality. `--max_dim` and the `--log` series are there for that. +### Feeding QPU samples in: `--subspace-from strings` + +`--strings-file` takes a plain text file, one `2*norb`-bit determinant per line — not a +counts JSON. Converting a qiskit-addon-sqd counts file means two steps, and the second +one is not optional: + +1. **Postselect** on the target Hamming weight per spin sector. `recover_configurations` + and friends do this in the SQD loop; doing it here keeps determinants with the wrong + electron count out of the subspace. +2. **Interleave the bits.** qiskit-addon-sqd emits `[beta | alpha]` — the two halves + concatenated — while GDB wants them interleaved, bit `2*i` alpha orbital `i` and bit + `2*i + 1` beta orbital `i`. The drivers expose the conversion as `interleave()`. + +```python +import json, pathlib, sys +sys.path.insert(0, ".") # run from examples/gdb/ +from run_gdb_heatbath import interleave + +counts = json.loads(pathlib.Path("../tpb/count_dict_h2o.json").read_text()) +norb = len(next(iter(counts))) // 2 +n_alpha = n_beta = 5 # the sector you are solving + +dets = { + interleave(key[norb:], key[:norb]) # key is [beta | alpha] + for key in counts + if key[norb:].count("1") == n_alpha and key[:norb].count("1") == n_beta +} +pathlib.Path("dets.txt").write_text("\n".join(sorted(dets)) + "\n") +``` + +A `set` drops duplicate samples, which GDB rejects rather than silently dedupes, and +`sorted()` gives the ordering the shard contract wants. Then: + +```bash +python run_gdb_heatbath.py \ + --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \ + --subspace-from strings --strings-file dets.txt \ + --cutoffs 1e-3,1e-4 +``` + +On the bundled 275-bitstring h2o file that is a 275-determinant sparse subspace at +−76.0723973374, which the ladder then grows — as opposed to the 75,625-determinant +product space TPB would build from the same samples. + +**Skipping the interleave does not fail cleanly.** Feeding the concatenated strings +straight through diagonalizes to a plausible-looking number (−66.8043 for the file +above, against −76.0724 done right) and only aborts later, inside the heatbath +expansion, with `std::out_of_range`. The electron-count check does not catch it either: +permuting bits preserves how many are set, so the density still sums to 10. Check the +energy against a known reference before trusting a first conversion. + ### The loop needs no amplitudes A natural question, since `gdb_diag` does not return the wavefunction: the ladder never From 5c1f78aee615314709cab27392c386e08d6efc71 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 17:53:30 -0400 Subject: [PATCH 20/22] examples/gdb: validate alpha/beta electron counts before solving Both drivers now compare every determinant's per-spin electron count against the FCIDUMP's NELEC/MS2 and refuse a mismatch, naming how many determinants are wrong, the first offender, and the two likely causes: a concatenated [beta | alpha] bit order where GDB wants them interleaved, or samples that were never postselected. This is worth a guard because the failure is not self-announcing. Feeding qiskit-addon-sqd's concatenated strings straight in diagonalizes to -66.8043 where the correct conversion gives -76.0724, then aborts inside the heatbath expansion with std::out_of_range -- a crash whose message says nothing about bit order. The occupation density cannot catch it: permuting bits preserves how many are set, so it still sums to the right electron count. Cost is not a concern. The popcount runs on the already-packed words via a 256-entry uint8 table, 37 ms per million determinants measured in isolation, and end to end it disappears into run variance: a 1M-determinant run took 76.5/76.9/76.6 s without the check and 77.5/76.6/76.6 s with it. numpy.bitwise_count would be ~9x faster but needs numpy 2.0 and this package's floor is 1.19, so the table is the portable choice. --skip-weight-check opts out for a deliberately mixed-sector subspace. 28 passed serial, 12 with --run-slow, 13 under mpirun -n 2. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 18 ++++--- examples/gdb/run_gdb_diag.py | 76 ++++++++++++++++++++++++++++ examples/gdb/run_gdb_heatbath.py | 71 ++++++++++++++++++++++++++ test/test_gdb_drivers.py | 85 ++++++++++++++++++++++++++++++++ 4 files changed, 244 insertions(+), 6 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index c0d0041..a89ea38 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -241,12 +241,18 @@ On the bundled 275-bitstring h2o file that is a 275-determinant sparse subspace −76.0723973374, which the ladder then grows — as opposed to the 75,625-determinant product space TPB would build from the same samples. -**Skipping the interleave does not fail cleanly.** Feeding the concatenated strings -straight through diagonalizes to a plausible-looking number (−66.8043 for the file -above, against −76.0724 done right) and only aborts later, inside the heatbath -expansion, with `std::out_of_range`. The electron-count check does not catch it either: -permuting bits preserves how many are set, so the density still sums to 10. Check the -energy against a known reference before trusting a first conversion. +**Both drivers check this for you.** Every determinant's alpha and beta electron +counts are compared against the FCIDUMP's `NELEC`/`MS2` before anything is solved, and a +mismatch is refused with the count of offending determinants and the first one's index. +Skipping the interleave is worth guarding because it does not fail cleanly on its own: +the concatenated strings diagonalize to a plausible-looking number (−66.8043 for the +file above, against −76.0724 done right) and only abort later, inside the heatbath +expansion, with `std::out_of_range`. The occupation density cannot catch it either — +permuting bits preserves how many are set, so it still sums to 10. + +The check costs about 0.04 s per million determinants, which is not measurable against +the read and pack that precede it; `--skip-weight-check` turns it off if you are +deliberately mixing spin sectors. ### The loop needs no amplitudes diff --git a/examples/gdb/run_gdb_diag.py b/examples/gdb/run_gdb_diag.py index 89b6fe3..3e4a834 100755 --- a/examples/gdb/run_gdb_diag.py +++ b/examples/gdb/run_gdb_diag.py @@ -76,6 +76,68 @@ _DEFAULT_DETFILES = ",".join(f"{_GDB_APP}/det{i}.txt" for i in range(4)) +# One 256-entry popcount table, applied to the packed determinant words by viewing +# them as bytes. numpy.bitwise_count would be ~9x faster but needs numpy 2.0, and +# this package's floor is 1.19; at 0.018 s per million determinants the table is +# already a ~6% addition to the read-and-pack this replaces nothing of. +_POPCOUNT = np.array([bin(i).count("1") for i in range(256)], dtype=np.uint8) + +# Interleaved layout: alpha sits on the even bit positions, beta on the odd ones. +_ALPHA_MASK = np.uint64(0x5555555555555555) +_BETA_MASK = np.uint64(0xAAAAAAAAAAAAAAAA) + + +def _popcount_rows(words): + """Bits set per row of a (n, words) uint64 array.""" + return _POPCOUNT[words.view(np.uint8)].reshape(words.shape[0], -1).sum( + axis=1, dtype=np.int64) + + +def spin_weights(det): + """Per-determinant (n_alpha, n_beta) for a packed interleaved determinant array.""" + words = np.ascontiguousarray(det, dtype=np.uint64) + if words.ndim == 1: + words = words.reshape(1, -1) + return _popcount_rows(words & _ALPHA_MASK), _popcount_rows(words & _BETA_MASK) + + +def check_spin_weights(det, nelec, ms2, bit_length=64): + """Refuse a determinant list whose electron counts per spin are not uniform. + + A wrong bit order is the failure this catches, and it is worth catching here + because it does not fail cleanly downstream: concatenated ``[beta | alpha]`` + strings (what qiskit-addon-sqd emits) diagonalize to a plausible-looking + energy and only abort later inside the heatbath expansion. The occupation + density cannot catch it either, since permuting bits preserves how many are + set, so it still sums to the right electron count. + + Costs about 0.018 s per million determinants -- a few percent of the read and + pack that precede it. + """ + if bit_length != 64: + return # the masks above assume 64-bit words + want_a = (nelec + ms2) // 2 + want_b = nelec - want_a + got_a, got_b = spin_weights(det) + bad = (got_a != want_a) | (got_b != want_b) + n_bad = int(bad.sum()) + if not n_bad: + return + first = int(np.argmax(bad)) + raise ValueError( + f"{n_bad} of {len(bad)} determinants do not have {want_a} alpha and " + f"{want_b} beta electrons (first at index {first}: " + f"{int(got_a[first])} alpha, {int(got_b[first])} beta).\n" + " The usual cause is bit order: GDB expects alpha and beta INTERLEAVED " + "(bit 2*i alpha orbital i, bit 2*i+1 beta orbital i), while " + "qiskit-addon-sqd emits them CONCATENATED as [beta | alpha]. Interleave " + "before packing -- see examples/gdb/README.md.\n" + " The other cause is samples that were never postselected on the target " + "Hamming weight.\n" + " Pass --skip-weight-check to proceed anyway." + ) + + def parse_args(): """Parse command line arguments for all GDB_SBD parameters.""" parser = argparse.ArgumentParser( @@ -90,6 +152,13 @@ def parse_args(): "so GDB runs on the host there") # --- input ------------------------------------------------------------ + parser.add_argument('--skip-weight-check', action='store_true', + dest='skip_weight_check', + help='Do not verify that every determinant has the expected ' + 'alpha/beta electron count. The check costs about ' + '0.02 s per million determinants and catches a wrong ' + 'bit order, which otherwise returns a plausible wrong ' + 'energy with no error at all') parser.add_argument('--fcidump', default=_DEFAULT_FCIDUMP, help='FCIDUMP file defining the Hamiltonian') parser.add_argument('--detfiles', default=_DEFAULT_DETFILES, @@ -260,6 +329,7 @@ def main(): fcidump = sbd.LoadFCIDump(args.fcidump) norb = int(fcidump.header["NORB"]) nelec = int(fcidump.header["NELEC"]) + ms2 = int(fcidump.header.get("MS2", 0)) total_bits = 2 * norb b_size = args.b_comm_size @@ -331,6 +401,12 @@ def main(): if det.shape[0] == 0: print("ERROR: no determinants read", file=sys.stderr) return 1 + if not args.skip_weight_check: + try: + check_spin_weights(det, nelec, ms2, args.bit_length) + except ValueError as exc: + print(f"ERROR: {exc}", file=sys.stderr) + return 1 words = det.shape[1] config = sbd.GDB_SBD() diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py index 3d771d2..6483358 100755 --- a/examples/gdb/run_gdb_heatbath.py +++ b/examples/gdb/run_gdb_heatbath.py @@ -88,6 +88,13 @@ def parse_args(): help='FCIDUMP file defining the Hamiltonian') # --- the seed --------------------------------------------------------- + parser.add_argument('--skip-weight-check', action='store_true', + dest='skip_weight_check', + help='Do not verify that every seed determinant has the ' + 'expected alpha/beta electron count. The check costs ' + 'about 0.02 s per million determinants and catches a ' + 'wrong bit order, which otherwise diagonalizes to a ' + 'plausible wrong energy and aborts during expansion') parser.add_argument('--subspace-from', default='files', dest='subspace_from', choices=['files', 'hf', 'from-alpha', 'strings'], @@ -178,6 +185,68 @@ def interleave(alpha, beta): return ''.join(a_rev[i] + b_rev[i] for i in range(len(a_rev)))[::-1] +# One 256-entry popcount table, applied to the packed determinant words by viewing +# them as bytes. numpy.bitwise_count would be ~9x faster but needs numpy 2.0, and +# this package's floor is 1.19; at 0.018 s per million determinants the table is +# already a ~6% addition to the read-and-pack this replaces nothing of. +_POPCOUNT = np.array([bin(i).count("1") for i in range(256)], dtype=np.uint8) + +# Interleaved layout: alpha sits on the even bit positions, beta on the odd ones. +_ALPHA_MASK = np.uint64(0x5555555555555555) +_BETA_MASK = np.uint64(0xAAAAAAAAAAAAAAAA) + + +def _popcount_rows(words): + """Bits set per row of a (n, words) uint64 array.""" + return _POPCOUNT[words.view(np.uint8)].reshape(words.shape[0], -1).sum( + axis=1, dtype=np.int64) + + +def spin_weights(det): + """Per-determinant (n_alpha, n_beta) for a packed interleaved determinant array.""" + words = np.ascontiguousarray(det, dtype=np.uint64) + if words.ndim == 1: + words = words.reshape(1, -1) + return _popcount_rows(words & _ALPHA_MASK), _popcount_rows(words & _BETA_MASK) + + +def check_spin_weights(det, nelec, ms2, bit_length=64): + """Refuse a determinant list whose electron counts per spin are not uniform. + + A wrong bit order is the failure this catches, and it is worth catching here + because it does not fail cleanly downstream: concatenated ``[beta | alpha]`` + strings (what qiskit-addon-sqd emits) diagonalize to a plausible-looking + energy and only abort later inside the heatbath expansion. The occupation + density cannot catch it either, since permuting bits preserves how many are + set, so it still sums to the right electron count. + + Costs about 0.018 s per million determinants -- a few percent of the read and + pack that precede it. + """ + if bit_length != 64: + return # the masks above assume 64-bit words + want_a = (nelec + ms2) // 2 + want_b = nelec - want_a + got_a, got_b = spin_weights(det) + bad = (got_a != want_a) | (got_b != want_b) + n_bad = int(bad.sum()) + if not n_bad: + return + first = int(np.argmax(bad)) + raise ValueError( + f"{n_bad} of {len(bad)} determinants do not have {want_a} alpha and " + f"{want_b} beta electrons (first at index {first}: " + f"{int(got_a[first])} alpha, {int(got_b[first])} beta).\n" + " The usual cause is bit order: GDB expects alpha and beta INTERLEAVED " + "(bit 2*i alpha orbital i, bit 2*i+1 beta orbital i), while " + "qiskit-addon-sqd emits them CONCATENATED as [beta | alpha]. Interleave " + "before packing -- see examples/gdb/README.md.\n" + " The other cause is samples that were never postselected on the target " + "Hamming weight.\n" + " Pass --skip-weight-check to proceed anyway." + ) + + def hartree_fock_string(norb, nelec, ms2): """The Hartree-Fock determinant: the lowest orbitals doubly occupied. @@ -316,6 +385,8 @@ def main(): try: det, source = build_seed(args, sbd, norb, nelec, ms2, rank) + if not args.skip_weight_check: + check_spin_weights(det, nelec, ms2, args.bit_length) except ValueError as exc: if rank == 0: print(f"ERROR: {exc}", file=sys.stderr) diff --git a/test/test_gdb_drivers.py b/test/test_gdb_drivers.py index b7dedbd..9eec318 100644 --- a/test/test_gdb_drivers.py +++ b/test/test_gdb_drivers.py @@ -249,3 +249,88 @@ def test_both_drivers_accept_either_alpha_flag_spelling(alpha_flag): ])) assert _dimension(stdout) == 576 assert _energy(stdout) == pytest.approx(H2O_576_ENERGY, abs=1e-8) + + +# --- spin-weight validation --------------------------------------------------- + +COUNTS_FILE = REPO_ROOT / "examples" / "tpb" / "count_dict_h2o.json" + + +def _interleave(alpha: str, beta: str) -> str: + """Bit 2*i alpha orbital i, bit 2*i+1 beta orbital i, counting from the right. + + Spelled out here rather than imported from the driver so the test does not + agree with the code under test by construction. + """ + a, b = alpha[::-1], beta[::-1] + return "".join(a[i] + b[i] for i in range(len(a)))[::-1] + + +def _counts_as_files(tmp_path): + """The bundled counts file written both correctly and incorrectly. + + qiskit-addon-sqd emits ``[beta | alpha]`` concatenated; GDB wants the two + interleaved. The "wrong" file is that raw concatenation, which is the mistake + a first conversion actually makes. + """ + counts = json.loads(COUNTS_FILE.read_text()) + norb = len(next(iter(counts))) // 2 + right = sorted({_interleave(k[norb:], k[:norb]) for k in counts}) + wrong = sorted(set(counts)) + good, bad = tmp_path / "right.txt", tmp_path / "wrong.txt" + good.write_text("\n".join(right) + "\n") + bad.write_text("\n".join(wrong) + "\n") + return good, bad + + +def _heatbath_on(strings_file, tmp_path, extra=()): + return _run("run_gdb_heatbath.py", [ + "--fcidump", str(H2O_DIR / "fcidump.txt"), + "--subspace-from", "strings", "--strings-file", str(strings_file), + "--cutoffs", "1e-3", "--max_rounds", "1", + "--log", str(tmp_path / "ladder.json"), *extra, + ]) + + +def test_interleaved_counts_are_accepted(tmp_path): + """The documented conversion of a counts file must run.""" + _require(COUNTS_FILE) + _require(H2O_DIR / "fcidump.txt") + good, _ = _counts_as_files(tmp_path) + _assert_ok(_heatbath_on(good, tmp_path)) + rounds = json.loads((tmp_path / "ladder.json").read_text())["rounds"] + assert rounds[0]["dimension"] == 275 + + +def test_concatenated_counts_are_refused_with_a_diagnosis(tmp_path): + """Feeding [beta | alpha] straight through must be caught before it is solved. + + Without the check this diagonalizes to a plausible-looking energy and only + aborts later inside the heatbath expansion with std::out_of_range, so the + failure gave no hint of its cause. The occupation density cannot catch it + either: permuting bits preserves how many are set. + """ + _require(COUNTS_FILE) + _require(H2O_DIR / "fcidump.txt") + _, bad = _counts_as_files(tmp_path) + completed = _heatbath_on(bad, tmp_path) + assert completed.returncode != 0, "mis-ordered determinants were accepted" + combined = completed.stdout + completed.stderr + assert "do not have 5 alpha and 5 beta electrons" in combined, combined[-2000:] + assert "INTERLEAVED" in combined, "the message should name the likely cause" + + +def test_skip_weight_check_bypasses_the_validation(tmp_path): + """The escape hatch must skip the check, not merely survive it. + + Asserted by the absence of the refusal: the run still fails downstream, since + the determinants really are malformed, but it fails in SBD rather than here. + """ + _require(COUNTS_FILE) + _require(H2O_DIR / "fcidump.txt") + _, bad = _counts_as_files(tmp_path) + completed = _heatbath_on(bad, tmp_path, extra=("--skip-weight-check",)) + combined = completed.stdout + completed.stderr + assert "do not have 5 alpha and 5 beta electrons" not in combined, ( + "--skip-weight-check did not skip the check" + ) From 1fb7a75cc33aa7af3ae7c1926bce30d2fe8e948f Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Tue, 29 Sep 2026 18:43:43 -0400 Subject: [PATCH 21/22] examples/gdb/README: the OMP_NUM_THREADS example only worked on Open MPI `mpirun -np 4 -x OMP_NUM_THREADS=8 ...` is Open MPI syntax. On MPICH, Hydra rejects it outright: [mpiexec] match_arg (lib/utils/args.c:166): unrecognized argument x so that example could not run at all on an MPICH system -- including the h100 box we validate on. Set the variable in the shell instead, which both implementations honour for a single-node run, and say why rather than leaving the next person to rediscover it. Verified on both: Open MPI 5.0.10 and MPICH 5.0.0 give -76.0588897208 with 144 determinants per rank at --b_comm_size 4. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index a89ea38..7c4e231 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -37,7 +37,7 @@ python run_gdb_diag.py --fcidump $GDB/fcidump_Fe4S4.txt \ # Shard the basis: Fe4S4's four files, one per rank. For this data that is # already the balanced globally-sorted split, so no redistribution is needed. -mpirun -np 4 -x OMP_NUM_THREADS=8 python run_gdb_diag.py --b_comm_size 4 +OMP_NUM_THREADS=8 mpirun -np 4 python run_gdb_diag.py --b_comm_size 4 # Spend ranks on both named dimensions (t <= b, and t*b must divide the ranks) mpirun -np 8 python run_gdb_diag.py --b_comm_size 4 --t_comm_size 2 @@ -356,6 +356,11 @@ that quotient by integer division and never checks the remainder, which silently produces communicators of unequal size and a rank alone in its own basis ring, so `gdb_diag` refuses it. +These examples set `OMP_NUM_THREADS` in the shell rather than through `mpirun`, because +the flag for that is implementation-specific — `-x VAR=VAL` on Open MPI, `-env VAR VAL` +on MPICH, and each rejects the other's spelling. Setting it in the shell works on both +for a single-node run. + ### What each dimension buys | | parallelizes | shards memory? | constraint | From f36798e2fbcaa0c1b2e038f0b4cba8a30aa97e65 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Wed, 30 Sep 2026 13:23:10 -0400 Subject: [PATCH 22/22] examples/gdb: warn when Davidson never iterated, and when an hf seed cannot shard Three things a real run hit that the drivers said nothing about. SBD's Davidson can return without iterating. If the start vector has no coupling to the rest of the subspace the residual is zero immediately, so it reports one determinant's diagonal element as the ground state energy: a valid upper bound, but not the subspace's eigenvalue, and the run exits 0 with a plausible number. The occupancies give it away -- a one-determinant wavefunction has every occupancy 0 or 1, where a correlated one is fractional -- so both drivers check that and say what happened. Guarded on dimension > 1, since a single-determinant subspace legitimately looks that way; the heatbath driver only checks its seed round, because after an expansion the subspace contains its own excitation neighbours. --subspace-from hf with --b_comm_size > 1 is a trap: one determinant cannot be divided, so rank 0 owns it, the other ranks get empty shards, and a single parent leaves OpenMP nothing to split either -- the first rounds crawl on one core. The driver now warns, and the README example runs serially. Also note in the README that the last cutoff dominates a ladder's cost, and that --max_dim is what makes an unreachable rung a clean stop rather than an out-of-memory failure mid-round. Verified: detector true on [1,1,0,0], false on fractional occupancies; a healthy h2o run emits nothing. 34 passed serial, 13 under mpirun -n 2. Co-Authored-By: Claude Opus 5 (1M context) --- examples/gdb/README.md | 15 ++++++++++++--- examples/gdb/run_gdb_diag.py | 18 ++++++++++++++++++ examples/gdb/run_gdb_heatbath.py | 29 ++++++++++++++++++++++++++++- 3 files changed, 58 insertions(+), 4 deletions(-) diff --git a/examples/gdb/README.md b/examples/gdb/README.md index 7c4e231..9bd3458 100644 --- a/examples/gdb/README.md +++ b/examples/gdb/README.md @@ -81,11 +81,13 @@ new candidates, so one round's result **is** the next round's subspace. # Fe4S4 from upstream's shipped subspace (the default data above), one cutoff python run_gdb_heatbath.py --cutoffs 1e-3 -# A ladder, stopping before the subspace passes 2M determinants +# A ladder, stopping before the subspace passes 2M determinants. The last cutoff +# dominates the cost, so --max_dim is the brake. python run_gdb_heatbath.py --cutoffs 1e-3,1e-4,1e-5 --max_dim 2000000 -# The no-input null: start from the Hartree-Fock determinant alone -python run_gdb_heatbath.py --subspace-from hf --cutoffs 1e-3,1e-4 +# The no-input null: the Hartree-Fock determinant alone. Run it serially -- one +# determinant cannot be sharded. +OMP_NUM_THREADS=48 python run_gdb_heatbath.py --subspace-from hf --cutoffs 1e-3,1e-4 # Sharded. At --b_comm_size == ranks with t=1 the expansion comes back already # sharded for the next round, so nothing has to be gathered. @@ -171,6 +173,13 @@ which would mean the subspace shrank or a round failed to converge. For h2o and FCI limits in that basis (−76.24377680 and −109.04874199) are hard ceilings a correct run can never cross, which makes a sparse run self-checking even without a reference. +### What the rungs cost + +Each rung's expansion is much larger than the one before, so the last cutoff you name +dominates the run — a rung is reachable or not rather than merely slow, since past some +point the expansion outgrows what can be diagonalized. `--max_dim` is what makes that a +clean stop with a `stop_reason` instead of running out of memory mid-round. + ### Seeds `--subspace-from` selects where the starting subspace comes from, so one driver produces diff --git a/examples/gdb/run_gdb_diag.py b/examples/gdb/run_gdb_diag.py index 3e4a834..8683e6d 100755 --- a/examples/gdb/run_gdb_diag.py +++ b/examples/gdb/run_gdb_diag.py @@ -138,6 +138,15 @@ def check_spin_weights(det, nelec, ms2, bit_length=64): ) +def looks_like_a_single_determinant(density, tolerance=1e-6): + """True when every occupancy is 0 or 1: a one-determinant wavefunction. + + For a subspace of more than one determinant that means Davidson never iterated -- + the start vector had no coupling to the rest, so the residual was zero at once. + """ + return all(min(abs(x), abs(1.0 - x)) < tolerance for x in density) + + def parse_args(): """Parse command line arguments for all GDB_SBD parameters.""" parser = argparse.ArgumentParser( @@ -491,6 +500,15 @@ def main(): combined = [density[2 * i] + density[2 * i + 1] for i in range(len(density) // 2)] print(f"Density: {np.round(combined, 6).tolist()}") print(f" (sums to {sum(combined):.6f}; should equal the electron count {nelec})") + if results["global_dim"] > 1 and looks_like_a_single_determinant(density): + print(" WARNING: every occupancy is 0 or 1, so the wavefunction is a single " + "determinant") + print(" and Davidson did not iterate (look for 'tol=0' above). The energy is " + "a diagonal") + print(" element -- a valid upper bound, not the subspace's ground state. " + "Expanding the") + print(" subspace once (--carryover_type 2) gives Davidson something to " + "couple.") carryover = results['carryover_det'] if args.carryover_type: diff --git a/examples/gdb/run_gdb_heatbath.py b/examples/gdb/run_gdb_heatbath.py index 6483358..9ec03c2 100755 --- a/examples/gdb/run_gdb_heatbath.py +++ b/examples/gdb/run_gdb_heatbath.py @@ -162,7 +162,17 @@ def parse_args(): parser.add_argument('--log', default='', metavar='FILE', help='Write the per-round (dimension, energy) series as JSON') - return parser.parse_args() + args = parser.parse_args() + + # One determinant cannot be sharded: rank 0 owns it, the rest get empty shards, + # and a single parent gives OpenMP nothing to divide either. + if args.subspace_from == 'hf' and args.b_comm_size > 1: + print("WARNING: --subspace-from hf is a single determinant, so " + f"--b_comm_size {args.b_comm_size} leaves the", file=sys.stderr) + print(" other ranks with empty shards and runs the first rounds on one " + "core. Use -np 1.", file=sys.stderr) + + return args def read_strings(paths): @@ -247,6 +257,15 @@ def check_spin_weights(det, nelec, ms2, bit_length=64): ) +def looks_like_a_single_determinant(density, tolerance=1e-6): + """True when every occupancy is 0 or 1: a one-determinant wavefunction. + + For a subspace of more than one determinant that means Davidson never iterated -- + the start vector had no coupling to the rest, so the residual was zero at once. + """ + return all(min(abs(x), abs(1.0 - x)) < tolerance for x in density) + + def hartree_fock_string(norb, nelec, ms2): """The Hartree-Fock determinant: the lowest orbitals doubly occupied. @@ -492,6 +511,14 @@ def global_dim(local): shown = "" if delta is None else f"{delta:+.6f}" print(f" {rung:>4} {round_index:>5} {dim:>12} {energy:>18.10f} " f"{shown:>12} {elapsed:>7.1f}", flush=True) + # Only the seed round can be a no-op; after that the subspace contains + # its own excitation neighbours. + if (rank == 0 and rung == 0 and round_index == 0 and dim > 1 + and looks_like_a_single_determinant(result["density"])): + print(" NOTE: the seed diagonalization returned a single " + "determinant, so Davidson did") + print(" not iterate there -- the ladder effectively starts from one " + "determinant.", flush=True) history.append({ "rung": rung, "cutoff": cutoff, "round": round_index, "dimension": dim, "energy": energy,