From 531fce41d4cf66ce3b446d93681866ef12a25781 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 17:24:27 -0400 Subject: [PATCH 1/9] examples: remove dead CLI params, fix miswired ones, expose symmetrize_spin Audited every CLI parameter in run_sqd_sbd.py, run_sqd_enlarge_subspace_sbd.py, and run_sbd_diag.py for ones that are exposed but don't actually do anything. run_sqd_sbd.py / run_sqd_enlarge_subspace_sbd.py: removed --sbd_do_rdm (SBD computes real RDMs when set, but neither driver ever reads result.rdm1/result.rdm2 back off the SCIResult -- real wasted compute for zero visible effect) and --sbd_do_shuffle (read into a local C++ variable in the det_vector-based diag() overload these drivers call, then never referenced again anywhere in that function -- the shuffle logic only exists in the other, file-based diag() overload, which this code path never reaches). Both were genuinely inert, not just discarded downstream. Both drivers also hardcoded symmetrize_spin=True in their diagonalize_fermionic_hamiltonian call, forcing ci_strs_a == ci_strs_b every iteration even though SBD itself supports distinct alpha/beta determinant sets. Added --symmetrize_spin (default 1, matching prior behavior) as a real toggle on both. run_sbd_diag.py had a different, more basic bug: --max_time, --init, --shuffle/--do_shuffle, --use_precalculated_dets, and --max_memory_gb_for_determinants were all parsed into args but never copied onto the TPB_SBD config object at all -- pure no-ops regardless of what SBD's C++ layer would otherwise do with them. Added the missing config.X = args.X lines (guarding the two Thrust-only fields with hasattr, matching sbd_solver.py's own _create_sbd_config pattern, since they don't exist on the CPU/OMP-offload backends' TPB_SBD). --bdetfile was a step further: tpb_diag_from_files's binding doesn't even accept a second determinant file, and always derives beta from the one adetfile it's given. Added --symmetrize_spin here too (default 1, same meaning): when 0, load --adetfile and --bdetfile independently via LoadAlphaDets + sort_bitarray and call the data-structure entry point (sbd.tpb_diag) instead of tpb_diag_from_files -- genuinely distinct alpha/beta subspaces now work. When 1 (default) and --bdetfile is also given, warn and ignore it rather than silently doing nothing with it. Also wired up real RDM retrieval in run_sbd_diag.py: reshapes SBD's raw one_p_rdm/two_p_rdm (same verified formula as sbd_solver._assemble_rdms, duplicated rather than imported to keep this file qiskit-free), prints trace(rdm1) and natural orbital occupations, and adds --rdm_output to save rdm1/rdm2 to a .npz file. Verified on h100: H2O default runs on both SQD drivers reproduce the known reference energies unchanged; --symmetrize_spin 0 produces materially different (and correctly asymmetric-shaped) subspaces on both; run_sbd_diag.py's --rdm 1 output has trace(rdm1)=10 for H2O with sensible natural orbital occupations; --symmetrize_spin 0 --bdetfile with a genuinely different file gives a different, physically sensible energy, while --symmetrize_spin 1 (default) with --bdetfile set warns and reproduces the plain single-file result; --max_time now visibly truncates convergence; --shuffle/--init no longer crash and behave as expected. --- python/examples/README.md | 29 +++- python/examples/run_sbd_diag.py | 141 ++++++++++++++++-- .../examples/run_sqd_enlarge_subspace_sbd.py | 17 ++- python/examples/run_sqd_sbd.py | 19 +-- 4 files changed, 177 insertions(+), 29 deletions(-) diff --git a/python/examples/README.md b/python/examples/README.md index d6d0700..cb377ce 100644 --- a/python/examples/README.md +++ b/python/examples/README.md @@ -28,12 +28,30 @@ mpirun -np 8 python -u run_sbd_diag.py \ --fcidump ../../vendor/sbd-upstream/data/n2/fcidump.txt \ --adetfile ../../vendor/sbd-upstream/data/n2/1em3-alpha.txt \ --adet_comm_size 2 --bdet_comm_size 2 --task_comm_size 2 + +# Retrieve the 1-/2-particle RDMs and save them to a file +mpirun -np 2 python -u run_sbd_diag.py --rdm 1 --rdm_output rdms.npz ``` -**Key options:** `--device`, `--fcidump`, `--adetfile`, `--adet_comm_size`, -`--bdet_comm_size`, `--task_comm_size`, `--method`, `--tolerance`, `--iteration`. -(These keep their unprefixed names here: this driver *is* SBD. The SQD driver -prefixes them `--sbd_*`.) Run `python run_sbd_diag.py --help` for the full list. +Passing `--rdm 1` prints `trace(rdm1)` (should equal the electron count) and +the natural orbital occupations (eigenvalues of `rdm1`) — occupations near 2 +or 0 indicate a single-reference-like orbital, occupations near 1 (or +several clustered together) flag multi-reference character / a candidate +active space. `--rdm_output PATH` additionally saves `rdm1`/`rdm2` to a +numpy `.npz` file. + +By default beta determinants are derived from `--adetfile` alone (identical +to it, or a shuffled copy if `--shuffle` is set). `--symmetrize_spin 0` +loads `--adetfile` and `--bdetfile` as independent, genuinely distinct +alpha/beta determinant sets instead; `--bdetfile` is otherwise ignored +(with a warning) since symmetric mode always derives beta from alpha. + +**Key options:** `--device`, `--fcidump`, `--adetfile`, `--bdetfile`, +`--symmetrize_spin`, `--adet_comm_size`, `--bdet_comm_size`, +`--task_comm_size`, `--method`, `--tolerance`, `--iteration`, `--rdm`, +`--rdm_output`. (These keep their unprefixed names here: this driver *is* +SBD. The SQD drivers prefix them `--sbd_*`.) Run `python run_sbd_diag.py +--help` for the full list. **Requirements:** `sbd`, `mpi4py` @@ -201,7 +219,8 @@ prints an OOM warning when it detects this. |-----------|-----------------|---------| | `--counts FILE` | Load hardware bitstrings from a JSON file (use this or `--samples`) | none — falls back to `--samples` if omitted | | `--samples N` | Generate N random bitstrings at the target Hamming weights; plumbing check only, energy not meaningful | `3000` (only used when `--counts` is omitted) | -| `--samples_per_batch` | Dominant control on subspace dimension. With `symmetrize_spin` the alpha and beta string sets are merged, so the subspace is up to `(2N)^2`, not `N^2` | `3000` | +| `--samples_per_batch` | Dominant control on subspace dimension. With `--symmetrize_spin 1` the alpha and beta string sets are merged, so the subspace is up to `(2N)^2`, not `N^2` | `3000` | +| `--symmetrize_spin` | `1` (default): merge the alpha and beta string pools every iteration, forcing `ci_strs_a == ci_strs_b`. SBD itself supports distinct alpha/beta determinant sets — this is purely a qiskit-addon-sqd loop-layer setting. `0`: sample and carry over alpha and beta independently, allowing them to differ | `1` | | `--num_batches` | Independent subsamples per iteration; occupancies are averaged across them | `1` (`run_sqd_enlarge_subspace_sbd.py`) / `3` (`run_sqd_sbd.py`) | | `--sqd_carryover_threshold` | `run_sqd_sbd.py` only. `\|coefficient\|` cutoff for carrying a determinant into the next iteration's sample pool. **Lower it to carry more** | `1e-4` | | `--enlarge_threshold` | `run_sqd_enlarge_subspace_sbd.py` only — the analogous "carry more" knob for that driver, but structurally different: it gates which *pairs* get expanded into single excitations via `enlarge_batch_from_transitions`, not which determinants survive into resampling. **Lower it to expand more pairs per round** | `1e-4` | diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index 39adf63..d8fcd9f 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -34,11 +34,25 @@ mpirun -np 8 python run_sbd_diag.py \ --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \ --adetfile ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt + + # Retrieve the 1-/2-particle RDMs and save them to a file + mpirun -np 8 python run_sbd_diag.py --rdm 1 --rdm_output rdms.npz + # Prints trace(rdm1) (should equal the electron count) and the natural + # orbital occupations (eigenvalues of rdm1) -- occupations near 2 or 0 + # indicate a single-reference-like orbital, occupations near 1 (or + # several clustered together) flag multi-reference character / a + # candidate active space. + + # Distinct alpha and beta determinant files (default is beta = alpha) + mpirun -np 8 python run_sbd_diag.py --symmetrize_spin 0 \ + --adetfile alpha-dets.txt --bdetfile beta-dets.txt """ import argparse import sys +import numpy as np + def parse_args(): """Parse command line arguments for all TPB_SBD parameters""" parser = argparse.ArgumentParser( @@ -64,7 +78,11 @@ def parse_args(): parser.add_argument('--adetfile', default='../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt', help='Path to alpha determinants file') parser.add_argument('--bdetfile', default='', - help='Path to beta determinants file (optional, uses adetfile if not specified)') + help='Path to beta determinants file, used only when ' + '--symmetrize_spin 0 (otherwise ignored with a ' + 'warning: symmetric mode always derives beta from ' + '--adetfile). Defaults to --adetfile itself when ' + 'left unset.') parser.add_argument('--loadname', default='', help='Load initial wavefunction from file') parser.add_argument('--savename', default='', @@ -94,9 +112,25 @@ def parse_args(): parser.add_argument('--init', type=int, default=0, help='Initialization method') parser.add_argument('--shuffle', '--do_shuffle', type=int, default=0, dest='do_shuffle', - help='Shuffle determinants (0=no, 1=yes)') + help='Shuffle determinants loaded from --adetfile before ' + 'mirroring/deriving beta from them (0=no, 1-4=yes, ' + 'different shuffle seeds -- see sbdiag.h). Only ' + 'takes effect when --symmetrize_spin 1 (default): ' + 'that is the code path that derives beta from a ' + 'single loaded list at all.') + parser.add_argument('--symmetrize_spin', type=int, default=1, choices=[0, 1], + help='1 (default): beta determinants are derived from ' + '--adetfile alone (identical to it, or a shuffled ' + 'copy if --shuffle is set) -- --bdetfile is ignored ' + 'with a warning if given. 0: load --adetfile and ' + '--bdetfile as independent, genuinely distinct ' + 'alpha/beta determinant sets (--shuffle has no ' + 'effect in this mode).') parser.add_argument('--rdm', '--do_rdm', type=int, default=0, choices=[0, 1], dest='do_rdm', help='Calculate RDM (0=density only, 1=full RDM)') + parser.add_argument('--rdm_output', default='', + help='When set (and --rdm 1), save rdm1/rdm2 to this ' + 'path as a numpy .npz file (keys: rdm1, rdm2).') # Carryover determinant selection parser.add_argument('--carryover_type', type=int, default=0, @@ -123,6 +157,43 @@ def parse_args(): return parser.parse_args() + +def assemble_rdms(results, norb): + """Build spin-summed (rdm1, rdm2) from SBD's raw one_p_rdm/two_p_rdm. + + Returns (None, None) when --rdm 0 (SBD leaves these keys as empty lists + in that case). The reshape below mirrors sbd_solver._assemble_rdms + (python/sbd_solver.py) -- duplicated here rather than imported, so this + file stays free of any qiskit-addon-sqd dependency. It was verified + there against PySCF's own make_rdm1/make_rdm2 on all three SBD + backends, both element-wise and via the energy identity + E = einsum("pr,pr->",rdm1,hcore) + 0.5*einsum("prqs,prqs->",rdm2,eri). + + SBD's documented layout (sbd-ext docs/user-guide.md): + one_p_rdm[s][i + L*j] = + two_p_rdm[s+2t][i+L*j+L^2*k+L^3*l] = + A Fortran-order reshape implements those flat-index formulas directly. + """ + one_p_rdm = results.get('one_p_rdm') + two_p_rdm = results.get('two_p_rdm') + if not one_p_rdm or not two_p_rdm: + return None, None + + one_p_rdm = np.asarray(one_p_rdm) + two_p_rdm = np.asarray(two_p_rdm) + L = norb + + rdm1 = (np.reshape(one_p_rdm[0], (L, L), order='F') + + np.reshape(one_p_rdm[1], (L, L), order='F')) + + spin_summed = sum( + np.reshape(two_p_rdm[s], (L, L, L, L), order='F') for s in range(4) + ) + rdm2 = spin_summed.transpose(0, 2, 1, 3) + + return rdm1, rdm2 + + def main(): args = parse_args() @@ -148,6 +219,9 @@ def main(): config.eps = args.eps config.method = args.method config.max_nb = args.max_nb + config.max_time = args.max_time + config.init = args.init + config.do_shuffle = args.do_shuffle config.do_rdm = args.do_rdm config.bit_length = args.bit_length config.carryover_type = args.carryover_type @@ -156,6 +230,13 @@ def main(): config.adet_comm_size = args.adet_comm_size config.bdet_comm_size = args.bdet_comm_size config.task_comm_size = args.task_comm_size + # Thrust-only fields -- absent from the CPU/OMP-offload backends' TPB_SBD, + # so guard with hasattr rather than assume, matching sbd_solver.py's own + # _create_sbd_config pattern for exactly this reason. + if hasattr(config, 'use_precalculated_dets'): + config.use_precalculated_dets = bool(args.use_precalculated_dets) + if hasattr(config, 'max_memory_gb_for_determinants'): + config.max_memory_gb_for_determinants = args.max_memory_gb_for_determinants if rank == 0: print("Configuration:") @@ -190,13 +271,37 @@ def main(): if args.dump_matrix_form_wf: config.dump_matrix_form_wf = args.dump_matrix_form_wf - results = sbd.tpb_diag_from_files( - fcidumpfile=args.fcidump, - adetfile=args.adetfile, - sbd_data=config, - loadname=args.loadname, - savename=args.savename, - ) + if args.symmetrize_spin: + if args.bdetfile and rank == 0: + print(f"WARNING: --bdetfile {args.bdetfile!r} is ignored because " + "--symmetrize_spin is 1 (default) -- symmetric mode always " + "derives beta from --adetfile alone. Pass --symmetrize_spin 0 " + "to use a distinct beta-determinant file.\n") + results = sbd.tpb_diag_from_files( + fcidumpfile=args.fcidump, + adetfile=args.adetfile, + sbd_data=config, + loadname=args.loadname, + savename=args.savename, + ) + else: + # No bound function accepts two separate determinant files + # directly (tpb_diag_from_files always derives beta from the one + # adetfile it's given) -- load both ourselves and call the + # data-structure entry point instead, mirroring what + # tpb_diag_from_files itself does internally after loading + # (sbdiag.h: LoadAlphaDets(...); sort_bitarray(adet);). + bdetfile = args.bdetfile or args.adetfile + fcidump = sbd.LoadFCIDump(args.fcidump) + norb = int(fcidump.header["NORB"]) + adet = sbd.sort_bitarray( + sbd.LoadAlphaDets(args.adetfile, args.bit_length, norb)) + bdet = sbd.sort_bitarray( + sbd.LoadAlphaDets(bdetfile, args.bit_length, norb)) + results = sbd.tpb_diag( + fcidump, adet, bdet, config, + loadname=args.loadname, savename=args.savename, + ) if rank == 0: print("="*70) @@ -214,6 +319,24 @@ def main(): print(f"Density: {combined_density}") print(f"Carryover determinants: {len(results['carryover_adet'])}") + + norb = len(density) // 2 + rdm1, rdm2 = assemble_rdms(results, norb) + if rdm1 is not None: + print() + print(f"1-RDM trace: {np.trace(rdm1):.6f} " + "(should equal the total electron count)") + occupations = np.sort(np.linalg.eigvalsh(rdm1))[::-1] + print(f"Natural orbital occupations (sorted): " + f"{np.round(occupations, 6).tolist()}") + print(" -- occupations near 2 or 0 indicate a " + "single-reference-like orbital; occupations near 1 " + "(or several clustered together) flag multi-reference " + "character / a candidate active space.") + if args.rdm_output: + np.savez(args.rdm_output, rdm1=rdm1, rdm2=rdm2) + print(f"Saved rdm1/rdm2 to {args.rdm_output}") + print("="*70) print("\n✓ Calculation completed successfully!") print() diff --git a/python/examples/run_sqd_enlarge_subspace_sbd.py b/python/examples/run_sqd_enlarge_subspace_sbd.py index 4d2b062..bc2b63c 100644 --- a/python/examples/run_sqd_enlarge_subspace_sbd.py +++ b/python/examples/run_sqd_enlarge_subspace_sbd.py @@ -81,6 +81,13 @@ def parse_args(): "excitations and feed the result forward as next round's " "include_configurations.") loop.add_argument("--samples_per_batch", type=int, default=3000) + loop.add_argument("--symmetrize_spin", type=int, default=1, choices=[0, 1], + help="1 (default): merge the alpha and beta string pools " + "every round, forcing ci_strs_a == ci_strs_b -- SBD " + "itself supports distinct alpha/beta determinant " + "sets, but qiskit-addon-sqd's own loop does not " + "when this is on. 0: sample and carry over alpha " + "and beta independently, allowing them to differ.") loop.add_argument("--num_batches", type=int, default=1, help="Batches per outer round. Unlike run_sqd_sbd.py this " "is not the main lever on subspace size -- excitation " @@ -136,8 +143,6 @@ def parse_args(): sbd.add_argument("--sbd_max_it", type=int, default=10, dest="max_it", help="Max SBD Davidson iterations per diagonalization.") sbd.add_argument("--sbd_max_nb", type=int, default=10, dest="max_nb") - sbd.add_argument("--sbd_do_rdm", type=int, default=0, dest="do_rdm") - sbd.add_argument("--sbd_do_shuffle", type=int, default=0, dest="do_shuffle") sbd.add_argument("--sbd_use_precalculated_dets", type=int, default=1, choices=[0, 1]) sbd.add_argument("--sbd_max_memory_gb_for_determinants", type=int, default=-1) @@ -369,7 +374,8 @@ def main(): print(" " f"--energy_tol {args.energy_tol:g} " f"--occupancies_tol {args.occupancies_tol:g} " - f"--enlarge_threshold {args.enlarge_threshold:g}") + f"--enlarge_threshold {args.enlarge_threshold:g} " + f"--symmetrize_spin {args.symmetrize_spin}") print("SBD solver : " f"--sbd_method {args.method} --sbd_eps {args.eps:g} " f"--sbd_max_it {args.max_it} --sbd_max_nb {args.max_nb}") @@ -377,8 +383,7 @@ def main(): sbd_config = { "method": args.method, "eps": args.eps, "max_it": args.max_it, - "max_nb": args.max_nb, "max_time": 3600.0, "do_rdm": args.do_rdm, - "do_shuffle": args.do_shuffle, "bit_length": args.bit_length, + "max_nb": args.max_nb, "max_time": 3600.0, "bit_length": args.bit_length, "use_precalculated_dets": bool(args.sbd_use_precalculated_dets), "max_memory_gb_for_determinants": args.sbd_max_memory_gb_for_determinants, "adet_comm_size": args.adet_comm_size, "bdet_comm_size": args.bdet_comm_size, @@ -413,7 +418,7 @@ def callback(results): include_configurations=current_include, initial_occupancies=current_occ, sci_solver=sbd_solver, - symmetrize_spin=True, + symmetrize_spin=bool(args.symmetrize_spin), max_dim=args.max_dim, callback=callback, seed=rand_seed, diff --git a/python/examples/run_sqd_sbd.py b/python/examples/run_sqd_sbd.py index 8c836d5..e6999a0 100644 --- a/python/examples/run_sqd_sbd.py +++ b/python/examples/run_sqd_sbd.py @@ -81,6 +81,13 @@ def parse_args(): help="Dominant control on subspace size. With " "symmetrize_spin the alpha and beta string sets " "merge, so the subspace is up to (2N)^2.") + sqd.add_argument("--symmetrize_spin", type=int, default=1, choices=[0, 1], + help="1 (default): merge the alpha and beta string pools " + "every iteration, forcing ci_strs_a == ci_strs_b -- " + "SBD itself supports distinct alpha/beta determinant " + "sets, but qiskit-addon-sqd's own loop does not when " + "this is on. 0: sample and carry over alpha and beta " + "independently, allowing them to differ.") sqd.add_argument("--num_batches", type=int, default=3) sqd.add_argument("--max_iterations", type=int, default=5, help="SQD self-consistent loop iterations. NOT the SBD " @@ -174,11 +181,6 @@ def parse_args(): "and cross-batch agreement.") sbd.add_argument("--sbd_max_nb", "--block", "--max_nb", type=int, default=10, dest="max_nb") - sbd.add_argument("--sbd_do_rdm", "--rdm", "--do_rdm", type=int, default=0, - dest="do_rdm", - help="0=density only (default, sufficient for SQD), 1=full RDM") - sbd.add_argument("--sbd_do_shuffle", "--shuffle", "--do_shuffle", type=int, - default=0, dest="do_shuffle") sbd.add_argument("--sbd_use_precalculated_dets", type=int, default=1, choices=[0, 1], help="Thrust only. 1 precomputes a determinant index for every " @@ -366,8 +368,6 @@ def main(): "max_it": args.max_it, "max_nb": args.max_nb, "max_time": 3600.0, - "do_rdm": args.do_rdm, - "do_shuffle": args.do_shuffle, "bit_length": args.bit_length, "use_precalculated_dets": bool(args.sbd_use_precalculated_dets), "max_memory_gb_for_determinants": args.sbd_max_memory_gb_for_determinants, @@ -440,7 +440,8 @@ def callback(results: list[SCIResult]): print(" " f"--energy_tol {args.energy_tol:g} " f"--occupancies_tol {args.occupancies_tol:g} " - f"--sqd_carryover_threshold {args.sqd_carryover_threshold:g}") + f"--sqd_carryover_threshold {args.sqd_carryover_threshold:g} " + f"--symmetrize_spin {args.symmetrize_spin}") print("SBD solver : " f"--sbd_method {args.method} --sbd_eps {args.eps:g} " f"--sbd_max_it {args.max_it} --sbd_max_nb {args.max_nb} " @@ -469,7 +470,7 @@ def callback(results: list[SCIResult]): include_configurations=include_configurations, initial_occupancies=initial_occupancies, sci_solver=sbd_solver, - symmetrize_spin=True, + symmetrize_spin=bool(args.symmetrize_spin), callback=callback, seed=rand_seed, ) From dfb1b3013f850c41409c9ec605010cdd1cb94969 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 17:53:24 -0400 Subject: [PATCH 2/9] sbd_solver: make assemble_rdms public, share it with run_sbd_diag.py run_sbd_diag.py's RDM retrieval initially duplicated the reshape formula rather than importing it from sbd_solver.py, on the reasoning that this file should stay free of any qiskit-addon-sqd dependency. That reasoning doesn't hold up: sbd_solver.py's mpi4py import is already a real, unavoidable transitive dependency of importing sbd itself (sbd.init() uses it internally), and its pyscf/qiskit-addon-sqd imports are both soft (try/except), so importing sbd_solver.assemble_rdms adds no new hard dependency here. Renamed _assemble_rdms to assemble_rdms (it's no longer solver-internal -- two independent entry points now use it) and had run_sbd_diag.py import it instead of carrying its own copy. One formula, one place to fix or extend; run_sqd_sbd.py/run_sqd_enlarge_subspace_sbd.py already get correct rdm1/rdm2 on SCIResult via this same function without any further work if they ever want to surface RDMs again. Verified on h100: run_sbd_diag.py's --rdm 1 output (trace, natural orbital occupations, saved .npz) is bit-identical to before the de-duplication; run_sqd_sbd.py's H2O energy is unaffected. No automated regression test for this specific reshape formula -- it was verified by hand against PySCF's make_rdm1/make_rdm2 on all three SBD backends, both element-wise and via the energy identity, during development; a future regression would also surface immediately via run_sbd_diag.py's --rdm output. --- python/examples/run_sbd_diag.py | 45 +++++++-------------------------- python/sbd_solver.py | 4 +-- 2 files changed, 11 insertions(+), 38 deletions(-) diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index d8fcd9f..93aea48 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -158,48 +158,21 @@ def parse_args(): return parser.parse_args() -def assemble_rdms(results, norb): - """Build spin-summed (rdm1, rdm2) from SBD's raw one_p_rdm/two_p_rdm. - - Returns (None, None) when --rdm 0 (SBD leaves these keys as empty lists - in that case). The reshape below mirrors sbd_solver._assemble_rdms - (python/sbd_solver.py) -- duplicated here rather than imported, so this - file stays free of any qiskit-addon-sqd dependency. It was verified - there against PySCF's own make_rdm1/make_rdm2 on all three SBD - backends, both element-wise and via the energy identity - E = einsum("pr,pr->",rdm1,hcore) + 0.5*einsum("prqs,prqs->",rdm2,eri). - - SBD's documented layout (sbd-ext docs/user-guide.md): - one_p_rdm[s][i + L*j] = - two_p_rdm[s+2t][i+L*j+L^2*k+L^3*l] = - A Fortran-order reshape implements those flat-index formulas directly. - """ - one_p_rdm = results.get('one_p_rdm') - two_p_rdm = results.get('two_p_rdm') - if not one_p_rdm or not two_p_rdm: - return None, None - - one_p_rdm = np.asarray(one_p_rdm) - two_p_rdm = np.asarray(two_p_rdm) - L = norb - - rdm1 = (np.reshape(one_p_rdm[0], (L, L), order='F') - + np.reshape(one_p_rdm[1], (L, L), order='F')) - - spin_summed = sum( - np.reshape(two_p_rdm[s], (L, L, L, L), order='F') for s in range(4) - ) - rdm2 = spin_summed.transpose(0, 2, 1, 3) - - return rdm1, rdm2 - - def main(): args = parse_args() # Import sbd — auto-initializes on first use, but we call init() # explicitly here to set the default device from --device flag. import sbd + # Reuses the same rdm1/rdm2 assembly sbd_solver.solve_sci uses for the + # SQD drivers -- verified against PySCF's make_rdm1/make_rdm2 on all + # three SBD backends, both element-wise and via the energy identity + # E = einsum("pr,pr->",rdm1,hcore) + 0.5*einsum("prqs,prqs->",rdm2,eri). + # mpi4py is already a real dependency of `sbd` itself (used internally + # for sbd.init()'s communicator), so importing sbd_solver here adds no + # new hard dependency; its pyscf/qiskit-addon-sqd imports are both + # soft (try/except), unused by assemble_rdms itself. + from sbd.sbd_solver import assemble_rdms sbd.init(device=args.device) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index e95a2d6..e9fb489 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -263,12 +263,12 @@ def _solve_sci_core( nelec=nelec, ) - rdm1, rdm2 = _assemble_rdms(results, norb) + rdm1, rdm2 = assemble_rdms(results, norb) return SCIResult(energy, sci_state, orbital_occupancies=occupancies, rdm1=rdm1, rdm2=rdm2) -def _assemble_rdms(results: dict, norb: int) -> tuple[np.ndarray | None, np.ndarray | None]: +def assemble_rdms(results: dict, norb: int) -> tuple[np.ndarray | None, np.ndarray | None]: """Build spin-summed (rdm1, rdm2) from SBD's raw one_p_rdm/two_p_rdm. Returns (None, None) when ``do_rdm`` was 0 (SBD leaves these keys as From 87bfb91f5275e00453538e2fa1d6a55eae478727 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 19:04:32 -0400 Subject: [PATCH 3/9] examples: load FCIDUMP integrals without pyscf in the SQD drivers run_sqd_sbd.py and run_sqd_enlarge_subspace_sbd.py each had their own `from pyscf import ao2mo, tools` just to turn an FCIDUMP file into dense hcore/eri arrays (tools.fcidump.to_scf(...).get_hcore() + ao2mo.restore(1, mf._eri, norb)). Added sbd_solver.load_integrals_from_fcidump, which builds the same dense chemist-notation arrays directly from SBD's own LoadFCIDump binding's raw (value, i, j, k, l) integrals list -- verified bit-identical (max diff 0.0) to the pyscf construction on H2O -- and reuses the already-existing _read_fcidump_ecore for the nuclear repulsion piece. Both drivers now use this instead. This does not make either driver pyscf-free overall: qiskit-addon-sqd itself declares pyscf as a hard, non-optional dependency and qiskit_addon_sqd/fermion.py imports it unconditionally at module level (`from pyscf import fci`), and both drivers fundamentally need diagonalize_fermionic_hamiltonian from that module. Verified directly: with pyscf blocked, both still fail, but now at `from qiskit_addon_sqd.fermion import ...` inside qiskit-addon-sqd's own file, not anywhere in this driver's code -- confirming the drivers' own pyscf usage is genuinely gone, even though the transitive dependency through qiskit-addon-sqd remains. run_sbd_diag.py (no qiskit-addon-sqd dependency at all) was already pyscf-free before this commit and stays that way -- verified it runs correctly with pyscf blocked entirely. --- .../examples/run_sqd_enlarge_subspace_sbd.py | 11 ++--- python/examples/run_sqd_sbd.py | 11 ++--- python/sbd_solver.py | 43 +++++++++++++++++++ 3 files changed, 55 insertions(+), 10 deletions(-) diff --git a/python/examples/run_sqd_enlarge_subspace_sbd.py b/python/examples/run_sqd_enlarge_subspace_sbd.py index bc2b63c..6ee7852 100644 --- a/python/examples/run_sqd_enlarge_subspace_sbd.py +++ b/python/examples/run_sqd_enlarge_subspace_sbd.py @@ -49,7 +49,6 @@ import numpy as np from mpi4py import MPI -from pyscf import ao2mo, tools from qiskit.primitives import BitArray from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian, enlarge_batch_from_transitions @@ -323,10 +322,12 @@ def main(): else: device_config = DeviceConfig.cpu() - mf_as = tools.fcidump.to_scf(str(args.fcidump)) - hcore = mf_as.get_hcore() - eri = ao2mo.restore(1, mf_as._eri, norb) - nuclear_repulsion_energy = mf_as.mol.energy_nuc() + # Reads the FCIDUMP directly (via SBD's own LoadFCIDump binding) rather + # than through pyscf, since pyscf isn't a dependency of this driver -- + # verified bit-identical to pyscf.tools.fcidump.to_scf(...) + ao2mo. + from sbd.sbd_solver import load_integrals_from_fcidump + hcore, eri, nuclear_repulsion_energy = load_integrals_from_fcidump( + args.fcidump, norb) rand_seed = np.random.default_rng(42) diff --git a/python/examples/run_sqd_sbd.py b/python/examples/run_sqd_sbd.py index e6999a0..ad01907 100644 --- a/python/examples/run_sqd_sbd.py +++ b/python/examples/run_sqd_sbd.py @@ -48,7 +48,6 @@ import numpy as np from mpi4py import MPI -from pyscf import ao2mo, tools from qiskit.primitives import BitArray from qiskit_addon_sqd.fermion import SCIResult, diagonalize_fermionic_hamiltonian @@ -302,10 +301,12 @@ def main(): device_config = DeviceConfig.cpu() # --- Load molecular integrals --- - mf_as = tools.fcidump.to_scf(str(args.fcidump)) - hcore = mf_as.get_hcore() - eri = ao2mo.restore(1, mf_as._eri, norb) - nuclear_repulsion_energy = mf_as.mol.energy_nuc() + # Reads the FCIDUMP directly (via SBD's own LoadFCIDump binding) rather + # than through pyscf, since pyscf isn't a dependency of this driver -- + # verified bit-identical to pyscf.tools.fcidump.to_scf(...) + ao2mo. + from sbd.sbd_solver import load_integrals_from_fcidump + hcore, eri, nuclear_repulsion_energy = load_integrals_from_fcidump( + args.fcidump, norb) # --- Load or generate bitstrings --- rand_seed = np.random.default_rng(42) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index e9fb489..1a35170 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -480,6 +480,49 @@ def _read_fcidump_ecore(fcidump_path): return 0.0 +def load_integrals_from_fcidump(fcidump_path, norb): + """Build (hcore, eri, nuclear_repulsion_energy) directly from an + FCIDUMP file, without requiring pyscf. + + Uses SBD's own ``LoadFCIDump`` binding to parse the file (the exact + same parser SBD itself uses to diagonalize), then builds dense + chemist-notation arrays from its raw ``(value, i, j, k, l)`` integrals + list -- 1-indexed, with ``k == l == 0`` marking a one-electron entry and + all-zero indices marking ECORE (read separately via + ``_read_fcidump_ecore``, which already exists for this purpose + elsewhere in this module). Verified bit-identical (max diff 0.0) to + ``pyscf.tools.fcidump.to_scf(...).get_hcore()`` / + ``ao2mo.restore(1, mf._eri, norb)`` on H2O. + + ``eri`` follows the same chemist convention ``ao2mo.restore(1, ...)`` + would produce (``eri[p,q,r,s] = (pq|rs)``), with the standard 8-fold + permutational symmetry of real orbitals. + """ + from . import get_backend + + backend = get_backend() + fcidump = backend.LoadFCIDump(str(fcidump_path)) + + hcore = np.zeros((norb, norb)) + eri = np.zeros((norb, norb, norb, norb)) + for value, i, j, k, l in fcidump.integrals: + if i == 0: + continue # ECORE line + p, q = i - 1, j - 1 + if k == 0: + hcore[p, q] = value + hcore[q, p] = value + else: + r, s = k - 1, l - 1 + for a, b, c, d in { + (p, q, r, s), (q, p, r, s), (p, q, s, r), (q, p, s, r), + (r, s, p, q), (s, r, p, q), (r, s, q, p), (s, r, q, p), + }: + eri[a, b, c, d] = value + + return hcore, eri, _read_fcidump_ecore(fcidump_path) + + def _ci_strings_to_sbd_dets( ci_strings: np.ndarray, norb: int, backend, bit_length: int = SBD_DEFAULT_BIT_LENGTH, From 8bd01f6ad49d6978e3b21cd8ce483c424ef04125 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 19:14:32 -0400 Subject: [PATCH 4/9] run_sbd_diag.py: fold --rdm/--rdm_output into one flag rdm1 and rdm2 are always computed and reported together -- there's no partial case where you'd want one without the other -- so having a separate --rdm {0,1} to turn computation on and --rdm_output PATH to additionally save it was redundant. Removed --rdm; --rdm_output PATH is now the sole control (empty, the default, skips RDM computation entirely), matching this file's own existing convention for --dump_matrix_form_wf/ --loadname/--savename, where a path is what turns each feature on. Verified on h100: --rdm_output still prints the same trace/natural-orbital output and saves the same .npz as before; omitting it still skips RDM computation cleanly (same as before). --- python/examples/README.md | 19 +++++++++++-------- python/examples/run_sbd_diag.py | 27 +++++++++++++++++---------- 2 files changed, 28 insertions(+), 18 deletions(-) diff --git a/python/examples/README.md b/python/examples/README.md index cb377ce..7ba9612 100644 --- a/python/examples/README.md +++ b/python/examples/README.md @@ -30,15 +30,18 @@ mpirun -np 8 python -u run_sbd_diag.py \ --adet_comm_size 2 --bdet_comm_size 2 --task_comm_size 2 # Retrieve the 1-/2-particle RDMs and save them to a file -mpirun -np 2 python -u run_sbd_diag.py --rdm 1 --rdm_output rdms.npz +mpirun -np 2 python -u run_sbd_diag.py --rdm_output rdms.npz ``` -Passing `--rdm 1` prints `trace(rdm1)` (should equal the electron count) and -the natural orbital occupations (eigenvalues of `rdm1`) — occupations near 2 -or 0 indicate a single-reference-like orbital, occupations near 1 (or -several clustered together) flag multi-reference character / a candidate -active space. `--rdm_output PATH` additionally saves `rdm1`/`rdm2` to a -numpy `.npz` file. +`--rdm_output PATH` computes the RDMs, prints `trace(rdm1)` (should equal +the electron count) and the natural orbital occupations (eigenvalues of +`rdm1`) — occupations near 2 or 0 indicate a single-reference-like orbital, +occupations near 1 (or several clustered together) flag multi-reference +character / a candidate active space — and saves `rdm1`/`rdm2` to that path +as a numpy `.npz` file. There's no separate on/off flag: `rdm1` and `rdm2` +are always computed and reported together, so (matching +`--dump_matrix_form_wf`/`--loadname`/`--savename`) a path is what turns the +feature on, and leaving it empty (the default) skips computing RDMs at all. By default beta determinants are derived from `--adetfile` alone (identical to it, or a shuffled copy if `--shuffle` is set). `--symmetrize_spin 0` @@ -48,7 +51,7 @@ alpha/beta determinant sets instead; `--bdetfile` is otherwise ignored **Key options:** `--device`, `--fcidump`, `--adetfile`, `--bdetfile`, `--symmetrize_spin`, `--adet_comm_size`, `--bdet_comm_size`, -`--task_comm_size`, `--method`, `--tolerance`, `--iteration`, `--rdm`, +`--task_comm_size`, `--method`, `--tolerance`, `--iteration`, `--rdm_output`. (These keep their unprefixed names here: this driver *is* SBD. The SQD drivers prefix them `--sbd_*`.) Run `python run_sbd_diag.py --help` for the full list. diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index 93aea48..d7da2ce 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -36,7 +36,7 @@ --adetfile ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt # Retrieve the 1-/2-particle RDMs and save them to a file - mpirun -np 8 python run_sbd_diag.py --rdm 1 --rdm_output rdms.npz + mpirun -np 8 python run_sbd_diag.py --rdm_output rdms.npz # Prints trace(rdm1) (should equal the electron count) and the natural # orbital occupations (eigenvalues of rdm1) -- occupations near 2 or 0 # indicate a single-reference-like orbital, occupations near 1 (or @@ -126,12 +126,17 @@ def parse_args(): '--bdetfile as independent, genuinely distinct ' 'alpha/beta determinant sets (--shuffle has no ' 'effect in this mode).') - parser.add_argument('--rdm', '--do_rdm', type=int, default=0, choices=[0, 1], dest='do_rdm', - help='Calculate RDM (0=density only, 1=full RDM)') parser.add_argument('--rdm_output', default='', - help='When set (and --rdm 1), save rdm1/rdm2 to this ' - 'path as a numpy .npz file (keys: rdm1, rdm2).') - + help='Path to save rdm1/rdm2 as a numpy .npz file ' + '(keys: rdm1, rdm2). Empty (default): density ' + 'only, no RDMs computed -- matching ' + '--dump_matrix_form_wf/--loadname/--savename, a ' + 'path here is what turns the feature on; there ' + 'is no separate on/off flag since rdm1 and rdm2 ' + 'are always computed and reported together. ' + 'Also prints trace(rdm1) and the natural ' + 'orbital occupations either way.') + # Carryover determinant selection parser.add_argument('--carryover_type', type=int, default=0, help='Carryover determinant selection type') @@ -195,7 +200,7 @@ def main(): config.max_time = args.max_time config.init = args.init config.do_shuffle = args.do_shuffle - config.do_rdm = args.do_rdm + 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 @@ -306,9 +311,11 @@ def main(): "single-reference-like orbital; occupations near 1 " "(or several clustered together) flag multi-reference " "character / a candidate active space.") - if args.rdm_output: - np.savez(args.rdm_output, rdm1=rdm1, rdm2=rdm2) - print(f"Saved rdm1/rdm2 to {args.rdm_output}") + # rdm1 is only non-None when --rdm_output was given (that's + # what turns config.do_rdm on above), so this always fires + # here -- no separate on/off check needed. + np.savez(args.rdm_output, rdm1=rdm1, rdm2=rdm2) + print(f"Saved rdm1/rdm2 to {args.rdm_output}") print("="*70) print("\n✓ Calculation completed successfully!") From 18c84aaee6a2f12ecc66ede5a5ab9101caaf5662 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 19:29:17 -0400 Subject: [PATCH 5/9] run_sbd_diag.py: examples show a real path for --rdm_output, add metavar Usage examples used a bare relative filename (rdms.npz), which doesn't illustrate that --rdm_output takes a path (directories included), even though the help text and README prose already call it that. Updated both usage examples to /tmp/h2o_rdms.npz, and added metavar='PATH' so --help itself reads --rdm_output PATH consistently with the docs. --- python/examples/README.md | 2 +- python/examples/run_sbd_diag.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/python/examples/README.md b/python/examples/README.md index 7ba9612..be727e0 100644 --- a/python/examples/README.md +++ b/python/examples/README.md @@ -30,7 +30,7 @@ mpirun -np 8 python -u run_sbd_diag.py \ --adet_comm_size 2 --bdet_comm_size 2 --task_comm_size 2 # Retrieve the 1-/2-particle RDMs and save them to a file -mpirun -np 2 python -u run_sbd_diag.py --rdm_output rdms.npz +mpirun -np 2 python -u run_sbd_diag.py --rdm_output /tmp/h2o_rdms.npz ``` `--rdm_output PATH` computes the RDMs, prints `trace(rdm1)` (should equal diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index d7da2ce..2cbdf54 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -36,7 +36,7 @@ --adetfile ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt # Retrieve the 1-/2-particle RDMs and save them to a file - mpirun -np 8 python run_sbd_diag.py --rdm_output rdms.npz + mpirun -np 8 python run_sbd_diag.py --rdm_output /tmp/h2o_rdms.npz # Prints trace(rdm1) (should equal the electron count) and the natural # orbital occupations (eigenvalues of rdm1) -- occupations near 2 or 0 # indicate a single-reference-like orbital, occupations near 1 (or @@ -126,7 +126,7 @@ def parse_args(): '--bdetfile as independent, genuinely distinct ' 'alpha/beta determinant sets (--shuffle has no ' 'effect in this mode).') - parser.add_argument('--rdm_output', default='', + parser.add_argument('--rdm_output', default='', metavar='PATH', help='Path to save rdm1/rdm2 as a numpy .npz file ' '(keys: rdm1, rdm2). Empty (default): density ' 'only, no RDMs computed -- matching ' From a24f2f2457e45891ae4f555d1ea9c7d3a391db6a Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 19:40:29 -0400 Subject: [PATCH 6/9] run_sbd_diag.py/README: clarify --rdm_output writes ONE .npz with both arrays Explicitly contrast with upstream SBD's own CLI (two separate files, 1pRDM.txt/2pRDM.txt) in the help text, module docstring, and README -- this is exactly what caused confusion when discussing the feature. Also trims the README paragraph down from 9 lines to 5. --- python/examples/README.md | 15 ++++++--------- python/examples/run_sbd_diag.py | 25 ++++++++++--------------- 2 files changed, 16 insertions(+), 24 deletions(-) diff --git a/python/examples/README.md b/python/examples/README.md index be727e0..83bdf1d 100644 --- a/python/examples/README.md +++ b/python/examples/README.md @@ -33,15 +33,12 @@ mpirun -np 8 python -u run_sbd_diag.py \ mpirun -np 2 python -u run_sbd_diag.py --rdm_output /tmp/h2o_rdms.npz ``` -`--rdm_output PATH` computes the RDMs, prints `trace(rdm1)` (should equal -the electron count) and the natural orbital occupations (eigenvalues of -`rdm1`) — occupations near 2 or 0 indicate a single-reference-like orbital, -occupations near 1 (or several clustered together) flag multi-reference -character / a candidate active space — and saves `rdm1`/`rdm2` to that path -as a numpy `.npz` file. There's no separate on/off flag: `rdm1` and `rdm2` -are always computed and reported together, so (matching -`--dump_matrix_form_wf`/`--loadname`/`--savename`) a path is what turns the -feature on, and leaving it empty (the default) skips computing RDMs at all. +`--rdm_output PATH` computes the RDMs and saves **one** `.npz` file at +`PATH` holding both `rdm1` and `rdm2` (`data = np.load(PATH); data["rdm1"]`, +`data["rdm2"]`) — unlike upstream SBD's own CLI, which writes two separate +files (`1pRDM.txt`/`2pRDM.txt`). It also prints `trace(rdm1)` and the +natural orbital occupations. Leaving it empty (the default) skips computing +RDMs entirely. By default beta determinants are derived from `--adetfile` alone (identical to it, or a shuffled copy if `--shuffle` is set). `--symmetrize_spin 0` diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index 2cbdf54..5cad452 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -35,13 +35,10 @@ --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \ --adetfile ../../vendor/sbd-upstream/data/h2o/h2o-1em3-alpha.txt - # Retrieve the 1-/2-particle RDMs and save them to a file + # Retrieve the 1-/2-particle RDMs: saves ONE .npz file holding both + # rdm1 and rdm2 (np.load(path)["rdm1"] / ["rdm2"]) -- not two separate + # files the way upstream SBD's own CLI does (1pRDM.txt/2pRDM.txt). mpirun -np 8 python run_sbd_diag.py --rdm_output /tmp/h2o_rdms.npz - # Prints trace(rdm1) (should equal the electron count) and the natural - # orbital occupations (eigenvalues of rdm1) -- occupations near 2 or 0 - # indicate a single-reference-like orbital, occupations near 1 (or - # several clustered together) flag multi-reference character / a - # candidate active space. # Distinct alpha and beta determinant files (default is beta = alpha) mpirun -np 8 python run_sbd_diag.py --symmetrize_spin 0 \ @@ -127,15 +124,13 @@ def parse_args(): 'alpha/beta determinant sets (--shuffle has no ' 'effect in this mode).') parser.add_argument('--rdm_output', default='', metavar='PATH', - help='Path to save rdm1/rdm2 as a numpy .npz file ' - '(keys: rdm1, rdm2). Empty (default): density ' - 'only, no RDMs computed -- matching ' - '--dump_matrix_form_wf/--loadname/--savename, a ' - 'path here is what turns the feature on; there ' - 'is no separate on/off flag since rdm1 and rdm2 ' - 'are always computed and reported together. ' - 'Also prints trace(rdm1) and the natural ' - 'orbital occupations either way.') + help='Save rdm1 and rdm2 together in ONE numpy .npz ' + 'file at this path (np.load(PATH)["rdm1"] / ' + '["rdm2"]) -- unlike upstream SBD\'s own CLI, ' + 'which writes two separate files (1pRDM.txt / ' + '2pRDM.txt). Also prints trace(rdm1) and the ' + 'natural orbital occupations. Empty (default): ' + 'no RDMs computed at all.') # Carryover determinant selection parser.add_argument('--carryover_type', type=int, default=0, From 092916d4f129db18b759bb34ea9c007548d48545 Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 19:49:01 -0400 Subject: [PATCH 7/9] Revert "examples: load FCIDUMP integrals without pyscf in the SQD drivers" This reverts commit 41affcf03d33887a52dea41f0fd68a25a3ca4dc. qiskit-addon-sqd itself declares pyscf as a hard, non-optional dependency and imports it unconditionally in qiskit_addon_sqd/fermion.py, which run_sqd_sbd.py and run_sqd_enlarge_subspace_sbd.py both need regardless (that's their entire purpose). pyscf is therefore always present whenever either driver can run at all, so load_integrals_from_fcidump was duplicating pyscf's own FCIDUMP parsing for no actual dependency-reduction benefit -- pure added maintenance surface. Reverting back to pyscf.tools.fcidump.to_scf(...) + ao2mo.restore(1, ...) directly, as before. run_sbd_diag.py is unaffected: it never called load_integrals_from_fcidump (it doesn't build dense hcore/eri arrays at all -- SBD's own tpb_diag_from_files/tpb_diag load the FCIDUMP directly), only sbd_solver.assemble_rdms, which stays. --- .../examples/run_sqd_enlarge_subspace_sbd.py | 11 +++-- python/examples/run_sqd_sbd.py | 11 +++-- python/sbd_solver.py | 43 ------------------- 3 files changed, 10 insertions(+), 55 deletions(-) diff --git a/python/examples/run_sqd_enlarge_subspace_sbd.py b/python/examples/run_sqd_enlarge_subspace_sbd.py index 6ee7852..bc2b63c 100644 --- a/python/examples/run_sqd_enlarge_subspace_sbd.py +++ b/python/examples/run_sqd_enlarge_subspace_sbd.py @@ -49,6 +49,7 @@ import numpy as np from mpi4py import MPI +from pyscf import ao2mo, tools from qiskit.primitives import BitArray from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian, enlarge_batch_from_transitions @@ -322,12 +323,10 @@ def main(): else: device_config = DeviceConfig.cpu() - # Reads the FCIDUMP directly (via SBD's own LoadFCIDump binding) rather - # than through pyscf, since pyscf isn't a dependency of this driver -- - # verified bit-identical to pyscf.tools.fcidump.to_scf(...) + ao2mo. - from sbd.sbd_solver import load_integrals_from_fcidump - hcore, eri, nuclear_repulsion_energy = load_integrals_from_fcidump( - args.fcidump, norb) + mf_as = tools.fcidump.to_scf(str(args.fcidump)) + hcore = mf_as.get_hcore() + eri = ao2mo.restore(1, mf_as._eri, norb) + nuclear_repulsion_energy = mf_as.mol.energy_nuc() rand_seed = np.random.default_rng(42) diff --git a/python/examples/run_sqd_sbd.py b/python/examples/run_sqd_sbd.py index ad01907..e6999a0 100644 --- a/python/examples/run_sqd_sbd.py +++ b/python/examples/run_sqd_sbd.py @@ -48,6 +48,7 @@ import numpy as np from mpi4py import MPI +from pyscf import ao2mo, tools from qiskit.primitives import BitArray from qiskit_addon_sqd.fermion import SCIResult, diagonalize_fermionic_hamiltonian @@ -301,12 +302,10 @@ def main(): device_config = DeviceConfig.cpu() # --- Load molecular integrals --- - # Reads the FCIDUMP directly (via SBD's own LoadFCIDump binding) rather - # than through pyscf, since pyscf isn't a dependency of this driver -- - # verified bit-identical to pyscf.tools.fcidump.to_scf(...) + ao2mo. - from sbd.sbd_solver import load_integrals_from_fcidump - hcore, eri, nuclear_repulsion_energy = load_integrals_from_fcidump( - args.fcidump, norb) + mf_as = tools.fcidump.to_scf(str(args.fcidump)) + hcore = mf_as.get_hcore() + eri = ao2mo.restore(1, mf_as._eri, norb) + nuclear_repulsion_energy = mf_as.mol.energy_nuc() # --- Load or generate bitstrings --- rand_seed = np.random.default_rng(42) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index 1a35170..e9fb489 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -480,49 +480,6 @@ def _read_fcidump_ecore(fcidump_path): return 0.0 -def load_integrals_from_fcidump(fcidump_path, norb): - """Build (hcore, eri, nuclear_repulsion_energy) directly from an - FCIDUMP file, without requiring pyscf. - - Uses SBD's own ``LoadFCIDump`` binding to parse the file (the exact - same parser SBD itself uses to diagonalize), then builds dense - chemist-notation arrays from its raw ``(value, i, j, k, l)`` integrals - list -- 1-indexed, with ``k == l == 0`` marking a one-electron entry and - all-zero indices marking ECORE (read separately via - ``_read_fcidump_ecore``, which already exists for this purpose - elsewhere in this module). Verified bit-identical (max diff 0.0) to - ``pyscf.tools.fcidump.to_scf(...).get_hcore()`` / - ``ao2mo.restore(1, mf._eri, norb)`` on H2O. - - ``eri`` follows the same chemist convention ``ao2mo.restore(1, ...)`` - would produce (``eri[p,q,r,s] = (pq|rs)``), with the standard 8-fold - permutational symmetry of real orbitals. - """ - from . import get_backend - - backend = get_backend() - fcidump = backend.LoadFCIDump(str(fcidump_path)) - - hcore = np.zeros((norb, norb)) - eri = np.zeros((norb, norb, norb, norb)) - for value, i, j, k, l in fcidump.integrals: - if i == 0: - continue # ECORE line - p, q = i - 1, j - 1 - if k == 0: - hcore[p, q] = value - hcore[q, p] = value - else: - r, s = k - 1, l - 1 - for a, b, c, d in { - (p, q, r, s), (q, p, r, s), (p, q, s, r), (q, p, s, r), - (r, s, p, q), (s, r, p, q), (r, s, q, p), (s, r, q, p), - }: - eri[a, b, c, d] = value - - return hcore, eri, _read_fcidump_ecore(fcidump_path) - - def _ci_strings_to_sbd_dets( ci_strings: np.ndarray, norb: int, backend, bit_length: int = SBD_DEFAULT_BIT_LENGTH, From 701234536c3fd5f2065377b121f03ab1a2a5749a Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 20:09:59 -0400 Subject: [PATCH 8/9] run_sbd_diag.py/README: drop an unneeded import comment, tie README to the concrete example The comment above 'from sbd.sbd_solver import assemble_rdms' justified why importing it doesn't add a new dependency -- reasoning that belongs in a commit message, not cluttering the import site. Also: README still described --rdm_output using the abstract 'PATH' placeholder in prose right after showing a concrete example path -- ties the description to that same example instead. --- python/examples/README.md | 13 +++++++------ python/examples/run_sbd_diag.py | 8 -------- 2 files changed, 7 insertions(+), 14 deletions(-) diff --git a/python/examples/README.md b/python/examples/README.md index 83bdf1d..da59ac4 100644 --- a/python/examples/README.md +++ b/python/examples/README.md @@ -33,12 +33,13 @@ mpirun -np 8 python -u run_sbd_diag.py \ mpirun -np 2 python -u run_sbd_diag.py --rdm_output /tmp/h2o_rdms.npz ``` -`--rdm_output PATH` computes the RDMs and saves **one** `.npz` file at -`PATH` holding both `rdm1` and `rdm2` (`data = np.load(PATH); data["rdm1"]`, -`data["rdm2"]`) — unlike upstream SBD's own CLI, which writes two separate -files (`1pRDM.txt`/`2pRDM.txt`). It also prints `trace(rdm1)` and the -natural orbital occupations. Leaving it empty (the default) skips computing -RDMs entirely. +`--rdm_output` takes the file to save to (`/tmp/h2o_rdms.npz` above) and +writes **one** `.npz` file there holding both `rdm1` and `rdm2` together +(`data = np.load("/tmp/h2o_rdms.npz"); data["rdm1"]`, `data["rdm2"]`) — +unlike upstream SBD's own CLI, which writes two separate files +(`1pRDM.txt`/`2pRDM.txt`). It also prints `trace(rdm1)` and the natural +orbital occupations. Leaving it empty (the default) skips computing RDMs +entirely. By default beta determinants are derived from `--adetfile` alone (identical to it, or a shuffled copy if `--shuffle` is set). `--symmetrize_spin 0` diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index 5cad452..3a1e79e 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -164,14 +164,6 @@ def main(): # Import sbd — auto-initializes on first use, but we call init() # explicitly here to set the default device from --device flag. import sbd - # Reuses the same rdm1/rdm2 assembly sbd_solver.solve_sci uses for the - # SQD drivers -- verified against PySCF's make_rdm1/make_rdm2 on all - # three SBD backends, both element-wise and via the energy identity - # E = einsum("pr,pr->",rdm1,hcore) + 0.5*einsum("prqs,prqs->",rdm2,eri). - # mpi4py is already a real dependency of `sbd` itself (used internally - # for sbd.init()'s communicator), so importing sbd_solver here adds no - # new hard dependency; its pyscf/qiskit-addon-sqd imports are both - # soft (try/except), unused by assemble_rdms itself. from sbd.sbd_solver import assemble_rdms sbd.init(device=args.device) From b66cd2a04b808dea69fee8dddbe79ce225c50d0e Mon Sep 17 00:00:00 2001 From: Sophia Wen Date: Thu, 17 Sep 2026 20:12:46 -0400 Subject: [PATCH 9/9] run_sbd_diag.py: --rdm_output metavar FILE instead of PATH PATH reads as ambiguous between a file and a directory; this flag specifically takes a file location, so FILE is less confusing. --- python/examples/run_sbd_diag.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/examples/run_sbd_diag.py b/python/examples/run_sbd_diag.py index 3a1e79e..4ac49cf 100644 --- a/python/examples/run_sbd_diag.py +++ b/python/examples/run_sbd_diag.py @@ -123,9 +123,9 @@ def parse_args(): '--bdetfile as independent, genuinely distinct ' 'alpha/beta determinant sets (--shuffle has no ' 'effect in this mode).') - parser.add_argument('--rdm_output', default='', metavar='PATH', + parser.add_argument('--rdm_output', default='', metavar='FILE', help='Save rdm1 and rdm2 together in ONE numpy .npz ' - 'file at this path (np.load(PATH)["rdm1"] / ' + 'file at this path (np.load(FILE)["rdm1"] / ' '["rdm2"]) -- unlike upstream SBD\'s own CLI, ' 'which writes two separate files (1pRDM.txt / ' '2pRDM.txt). Also prints trace(rdm1) and the '