Skip to content

README.md

Latest commit

 

History

History
515 lines (407 loc) · 29.4 KB

README.md

File metadata and controls

515 lines (407 loc) · 29.4 KB
Raw

TPB examples — SQD/SBD with a tensor-product basis

Examples demonstrating SBD's capabilities for quantum chemistry calculations. For how they fit with GDB and the variance check, see the workflows.

Overview

These examples all use TPB, where the subspace is the Cartesian product of an alpha and a beta determinant list.

The FCIDUMP may be restricted (RHF, ROHF) or unrestricted (UHF): sbd.tpb_diag, sbd.gdb_diag and the drivers read its layout and use the matching build of SBD, with no option to set. See Unrestricted (UHF) FCIDUMPs.

Examples

1. run_sbd_diag.py — Standalone SBD Diagonalization

Runs a single TPB diagonalization from an FCIDUMP file and alpha determinant file. No SQD loop, no Qiskit dependency.

# H2O (the default data) on 4 CPU ranks, a 2 x 2 alpha x beta grid
OMP_NUM_THREADS=8 mpirun -np 4 python -u run_sbd_diag.py \
    --adet_comm_size 2 --bdet_comm_size 2

# N2 on 8 GPUs, one rank each: the same command with --device gpu
mpirun -np 8 python -u run_sbd_diag.py --device gpu \
    --fcidump ../../vendor/sbd-upstream/data/n2/fcidump.txt \
    --adetfile ../../vendor/sbd-upstream/data/n2/1em3-alpha.txt \
    --adet_comm_size 4 --bdet_comm_size 2

# Also save the 1- and 2-particle RDMs
mpirun -np 4 python -u run_sbd_diag.py --adet_comm_size 2 --bdet_comm_size 2 \
    --rdm_output /tmp/h2o_rdms.npz

# Write the carryover SBD selects, then start the next run from it
mpirun -np 4 python -u run_sbd_diag.py --adet_comm_size 2 --bdet_comm_size 2 \
    --carryover_type 1 --carryover_ratio 0.5 \
    --carryover_adetfile /tmp/next-a.txt --carryover_bdetfile /tmp/next-b.txt
mpirun -np 4 python -u run_sbd_diag.py --adet_comm_size 2 --bdet_comm_size 2 \
    --symmetrize_spin 0 --adetfile /tmp/next-a.txt --bdetfile /tmp/next-b.txt

--rdm_output writes one .npz holding both rdm1 and rdm2 (np.load(path)["rdm1"]), where upstream's CLI writes two text files, and prints trace(rdm1) and the natural occupations. Without it, RDMs are not computed.

By default the beta list is the alpha one. --symmetrize_spin 0 reads --adetfile and --bdetfile as two separate lists; without it, --bdetfile is ignored with a warning.

Key options: --device, --fcidump, --adetfile, --bdetfile, --symmetrize_spin, --adet_comm_size, --bdet_comm_size, --task_comm_size, --method, --tolerance, --iteration, --initial_adet/--initial_bdet (the Davidson start; it must be in the subspace), --rdm_output. The SQD drivers prefix the solver options --sbd_*. python run_sbd_diag.py --help lists them all.

Requirements: sbd, mpi4py

2. run_sqd_sbd.py — SQD Loop with SBD Solver

Runs the self-consistent SQD workflow (qiskit-addon-sqd) using SBD as the eigensolver backend. Supports two bitstring input modes:

  • --counts FILE — load bitstrings from a count_dict.json
  • --samples N — generate N random bitstrings at the target Hamming weights (default). A plumbing check only: random determinants give a random subspace, so the energy is not meaningful. Use --counts for real results.
# H2O with the bundled counts file (275 bitstrings -> ~ -76.236 Ha)
mpirun -np 8 python -u run_sqd_sbd.py \
    --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \
    --counts count_dict_h2o.json \
    --device gpu \
    --adet_comm_size 4 --bdet_comm_size 2

# Random bitstrings for your own FCIDUMP: a plumbing check, not a result
mpirun -np 8 python -u run_sqd_sbd.py --fcidump /path/to/fci_dump.txt --device gpu \
    --samples_per_batch 800 --num_batches 3 --max_iterations 10 \
    --adet_comm_size 4 --bdet_comm_size 2

count_dict.json format: A JSON object mapping bitstrings to shot counts, as produced by a quantum device or simulator. Each bitstring has length 2 × NORB and is laid out as [beta | alpha]: the first NORB bits are beta (spin-down), the last NORB are alpha (spin-up), and within each half orbital 0 is the rightmost bit. qiskit-addon-sqd postselects the last NORB bits on num_elec_a and the first NORB on num_elec_b.

For H2O (NORB=24, 5α+5β) the Hartree–Fock configuration — the five lowest orbitals doubly occupied — is therefore "0"*19 + "1"*5 in both halves:

{
  "000000000000000000011111000000000000000000011111": 16,
  "010000000010001010000001010000000001000010100100": 12,
  "000010001110000000000010001001000110000000000100": 8
}

Bitstrings whose halves do not hold exactly num_elec_a / num_elec_b ones are dropped by postselection, so a file of uniform-random strings yields nothing usable — for H2O only C(24,5)² / 4²⁴ ≈ 6e-6 of them qualify.

count_dict_h2o.json in this directory is a ready-made H2O example: 275 bitstrings taken from the vendored h2o-1em3-alpha.txt determinant list, giving a 275 × 275 = 75,625-determinant subspace at ≈ -76.236 Ha.

Key options: --fcidump (required), --counts, --samples, --samples_per_batch, --num_batches, --max_iterations, --device, MPI decomposition flags, and the SQD tolerances --energy_tol / --occupancies_tol / --sqd_carryover_threshold. Inner-solver flags are prefixed (--sbd_eps, --sbd_max_it, --sbd_method, ...) and have sensible defaults; the old unprefixed spellings still work. --sbd_carryover_type lets SBD choose what carries into the next iteration instead of SQD; see SBD-selected carryover. Run python run_sqd_sbd.py --help for the full list, which is grouped by layer.

Requirements: see SBD in an SQD loop in the top-level README (pyscf, qiskit, qiskit-addon-sqd). --sbd_carryover_type needs qiskit-addon-sqd >= 0.14.0; everything else works from 0.13.1.

See SQD Parameters below for the full reference, grouped by SQD loop / SBD solver / MPI grid / checkpointing.

3. run_sqd_enlarge_subspace_sbd.py — SQD that grows its own subspace

The SQD loop of run_sqd_sbd.py, plus one step between rounds: it expands the dominant determinants of the state just solved with qiskit-addon-sqd's enlarge_batch_from_transitions (same-spin single excitations) and passes them to the next round as include_configurations. It runs one SQD iteration per round in its own loop, and stops when the expansion adds nothing new or the energy and occupancies stop moving (--energy_tol/--occupancies_tol); --max_iterations is only a cap.

mpirun -np 8 python -u run_sqd_enlarge_subspace_sbd.py \
    --fcidump ../../vendor/sbd-upstream/data/h2o/fcidump.txt \
    --counts count_dict_h2o.json \
    --device gpu \
    --adet_comm_size 4 --bdet_comm_size 2 --enlarge_threshold 1e-4

JAX on GPUs is handled for you. The expansion step is JAX, which by default reserves 75% of a GPU on first use, so several ranks would collide on one card. On a GPU backend the driver puts each rank's JAX on the card SBD chose for it (sbd.get_device_id()) and sets XLA_PYTHON_CLIENT_PREALLOCATE=false, so the command above needs no extra environment. A value you set yourself is kept, and JAX_PLATFORMS=cpu moves the expansion to the CPU.

Same bundled 275-bitstring H2O pool as run_sqd_sbd.py's own example above: plain SQD reaches ≈ -76.236 Ha and stops there; this driver keeps going past that fixed pool on its own and converges to -76.2421767512 Ha.

See SQD Parameters below for the flags it shares with run_sqd_sbd.py and the ones that differ (--enlarge_threshold in place of --sqd_carryover_threshold, and --max_dim's risk profile is sharper here).

4. run_sqd_sbd.ipynb — Jupyter walkthrough (serial)

Interactive single-rank companion to run_sqd_sbd.py. Same SQD self-consistent loop on h2o, but inside a Jupyter kernel (MPI.COMM_WORLD size 1). Uses the bundled count_dict_h2o.json (275 bitstrings → 75,625 determinants) and reaches ≈ −76.236 Ha in a few seconds on CPU.

pytest --nbmake run_sqd_sbd.ipynb      # what CI runs; needs the nbtest extra
# or open it in JupyterLab and step through the cells

SQD Parameters

Reference for every flag run_sqd_sbd.py accepts, grouped the way --help groups them: SQD loop, SBD solver, MPI grid, checkpointing.

How each iteration builds its subspace. SQD samples bitstrings from a quantum device, repairs the noisy ones against an orbital-occupancy estimate (configuration recovery), subsamples them into batches, and diagonalizes each batch. What makes it a loop is that two results feed back into the next iteration. Three sources, concatenated in this priority order inside qiskit-addon-sqd's diagonalize_fermionic_hamiltonian (fermion.py):

strs_a = include_a  ++  carryover_strings_a  ++  samples_a    then dedupe, truncate to max_dim, sort
  1. include_a — configurations passed as include_configurations. Static: fixed before the loop, present every iteration, never updated.
  2. carryover_strings_a — from the previous iteration's result, in descending order of marginal weight (a string's |c|^2 summed over the other spin sector). By default SQD selects them: every determinant whose |coefficient| is at least --sqd_carryover_threshold survives. With --sbd_carryover_type set, SBD selects them instead and SQD takes that selection as given; see SBD-selected carryover.
  3. samples_a — drawn fresh this iteration, sorted by how often each was sampled.

The order matters when max_dim is set: include and carryover are kept ahead of fresh samples, so if those two already fill the cap, this iteration's new samples are truncated away entirely.

The samples are not re-used raw counts. Each iteration re-runs configuration recovery from the original bitstrings using the occupancies from the previous iteration's best batch (_prepare_ci_strings in fermion.py), then subsamples. Recovery is not cumulative — it always re-derives from the raw samples, just with a better occupancy estimate each time. On iteration 1 there are no occupancies yet, so the raw samples are only filtered by electron count (Hamming-weight postselection).

So exactly two things flow from iteration N to N+1, and neither is a tolerance: the average orbital occupancies (into recovery, source 3) and the carryover (source 2). By default the carryover is derived from the wavefunction amplitudes; with --sbd_carryover_type 1 SBD hands it over directly and the amplitudes never leave SBD.

SQD loop parameters

Shared by both run_sqd_sbd.py and run_sqd_enlarge_subspace_sbd.py except where noted. --max_dim is the one that most needs attention: it has no universally safe default (see below), and in run_sqd_enlarge_subspace_sbd.py leaving it unset is riskier still, since each round's subspace can grow from the previous one rather than being resampled at a fixed size — the driver prints an OOM warning when it detects this.

Shapes the subspace — changes the numbers you compute:

Parameter What it controls Default
--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 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. Has no effect when --sbd_carryover_type is non-zero (the driver warns) 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
--max_dim Critical. Cap on strings per spin sector, so the subspace cannot exceed max_dim^2. The main brake on runaway cost — no fixed value is safe for every system, since the right cap depends on available memory and orbital count. Start from a value known to work at a similar orbital count (e.g. 15000 was used for a 45-orbital system) and adjust down if you see an OOM unset (no cap)
--include_hf Force the single Slater determinant with the lowest num_elec_a/num_elec_b orbital indices occupied into include_configurations, every iteration. Cheap correctness check: that determinant's own diagonal energy is an exact lower bound on what a subspace containing it can do — if forcing it in moves the result, the sampled pool was missing it (and probably its low-excitation neighbors too) off

Decides when to stop — changes nothing about the subspace:

Parameter What it controls Default
--max_iterations Hard cap on loop iterations (not the inner --sbd_max_it). In run_sqd_enlarge_subspace_sbd.py this is a safety cap only — the loop normally stops earlier, once a round adds no new determinants or both tolerances below are met 30 (run_sqd_enlarge_subspace_sbd.py) / 5 (run_sqd_sbd.py)
--energy_tol Iteration-to-iteration change in energy 1e-8
--occupancies_tol Largest change in any single orbital occupancy — an infinity norm, not an average 1e-5

Both stopping criteria must hold in the same iteration. fermion.py's convergence check combines the energy-change test and the occupancy-change test with a logical and, so the loop only stops once both are satisfied at once — not whichever one happens first. A run that reaches --max_iterations may be converged in energy while one stubborn orbital's occupancy is still moving, and loosening only one tolerance will not stop it. Watch the per-batch energies: while they still disagree, the loop has not converged regardless of what the total says.

SBD solver parameters

Per diagonalization, not per loop:

Parameter What it controls Default
--sbd_eps Davidson stop: norm of the residual vector, not an energy. Error in the energy goes roughly as |R|^2/gap, so this already implies far better energy accuracy than --energy_tol asks for. Tighten it for a near-degenerate system 1e-5
--sbd_max_it Cap on Davidson iterations. Reaching it before --sbd_eps returns a partially converged vector with no warning — watch the tol= values SBD prints, and cross-batch agreement 10
--sbd_max_nb Davidson basis vectors (block size). Peak memory during Davidson grows with how many sub-iterations it actually needs, not just --max_dim — a run that succeeds for several iterations at a fixed dim can still later need more basis vectors and run out of memory even though the subspace itself did not grow. Lowering this trades some convergence robustness for a lower memory ceiling 10
--sbd_method 0=Davidson, 1=Davidson+Ham, 2=Lanczos, 3=Lanczos+Ham 0
--sbd_use_precalculated_dets Thrust only. 1 precomputes a determinant index for every (α,β) pair — the whole subspace, on the GPU. 0 uses per-thread storage: slower per matvec, far less memory 1
--sbd_max_memory_gb_for_determinants Thrust only, and only consulted when --sbd_use_precalculated_dets 0 (mult_thrust.h:257-273). Caps the per-thread buffer in GB -1 (uncapped)
--sbd_thrust_collapse_loops Thrust only, and only consulted when --sbd_use_precalculated_dets 1. 1 runs each task's loops as one flattened GPU loop; 0 as nested loops, which need less GPU memory. Upstream switches to 0 by itself if the flattened size would overflow. Energies are unchanged 1
--sbd_cpu_subspace Thrust only. 1 keeps the Davidson subspace vectors in CPU memory instead of GPU memory; H is still applied on the GPU. For runs whose subspace vectors do not fit in GPU memory. Needs OMP_NUM_THREADS of 2 or more; ignored on the other backends 0

For reference, upstream's own TPB_SBD struct defaults are looser still (max_it=1, eps=1e-4), and run_sbd_diag.py uses eps=1e-3. On the h2o counts case, eps=1e-5 and eps=1e-8 give the same energy to ten decimal places.

SBD-selected carryover

By default SQD picks the next iteration's carryover by thresholding the amplitudes SBD returns. These flags let SBD pick it instead, using its own carryover rules, and report the result through SCIResult.carryover (qiskit-addon-sqd >= 0.14.0). SQD then uses SBD's selection as given: --sqd_carryover_threshold stops applying, and the order SBD's carryover comes in is the order a --max_dim truncation keeps.

Parameter What it controls Default
--sbd_carryover_type 0: off, SQD selects. 1: rank half-determinants by marginal weight. 2: type 1 plus all single excitations of what it selected. 3: rank whole determinants by amplitude, then single-excite 0
--sbd_carryover_ratio Types 1 and 2: fraction of half-determinants to keep, used only when non-zero. Ignored by type 3 0.1 (upstream SBD's own default is 0.0)
--sbd_carryover_threshold Types 1 and 2: cumulative-weight cutoff, used only when the ratio is 0. Type 3: a direct cutoff on determinant weight 1e-4

SBD reads exactly one of ratio and threshold for each type and silently ignores the other, so with the defaults, --sbd_carryover_threshold does nothing for types 1 and 2. The driver warns when you pass a setting that cannot take effect.

What each type costs. Type 1 is the one that saves work. Its ranking is exactly the marginal weight SQD orders carryover by, so nothing else needs the eigenvector, and the driver skips the |adet| × |bdet| wavefunction dump entirely. At large subspace size that dump is a multi-GB write-then-read every iteration. Types 2 and 3 still write the dump. They add determinants by single excitation, and SBD re-sorts the result into its canonical order (which it needs for binary search), so the driver re-ranks their carryover from the amplitudes: strings that were in the subspace come first, by weight, and the generated ones, which have no weight yet, come after.

A checkpoint written with --sbd_carryover_type set holds SBD's carryover strings, not the whole solved subspace, so --resume_from seeds from that smaller set.

Using a different subspace policy (Trim SQD)

qiskit-addon-sqd >= 0.14.0 lets the loop run a different per-iteration schedule, passed as a policy. The driver always uses the default, StandardPolicy, and does not expose others as flags. To try Trim SQD (qiskit_addon_sqd.trim.TrimPolicy), pass it from Python, with SBD as the solver as usual:

from functools import partial

from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
from qiskit_addon_sqd.trim import TrimPolicy
from sbd.sbd_solver import solve_sci_batch

num_batches = 4
policy = TrimPolicy(
    trim_ratio=1 / num_batches,   # see the note on trim_ratio below
    carryover_threshold=1e-4,     # the threshold goes INSIDE the policy
)

result = diagonalize_fermionic_hamiltonian(
    hcore, eri, bit_array,
    samples_per_batch=3000, num_batches=num_batches,
    norb=norb, nelec=nelec, symmetrize_spin=True,
    sci_solver=partial(solve_sci_batch, sbd_config={"eps": 1e-5, "max_it": 10}),
    policy=policy,
)
# Pass the threshold to the policy only. Also passing carryover_threshold= here
# raises a ValueError.

What changes compared to the default:

  • The batches are disjoint. One pool of samples_per_batch × num_batches bitstrings is divided among them, instead of each batch being subsampled independently.
  • Each iteration costs num_batches + 1 SBD calls, not num_batches. Every batch is diagonalized, its highest-weight strings are kept, the survivors of all batches are merged, and the merged subspace is diagonalized once more. That last diagonalization's energy is what the iteration reports.
  • trim_ratio sets how big the merged subspace is. It holds roughly trim_ratio × num_batches of one batch's strings per spin sector (less where batches share strings). TrimPolicy's default of 0.1 suits about 10 batches; with fewer, the merged subspace comes out much smaller than a batch. 1 / num_batches keeps it about one batch in size, which is the regime Trim SQD is designed for. That is our rule of thumb, not a tuned value. The addon's defaults are not tuned either.
  • With carryover_type != 0 in sbd_config, the policy's carryover settings do nothing. The trim schedule still runs, but each batch is screened by SBD's own selection rather than by trim_ratio, so tune ratio / threshold in sbd_config instead.

Trim SQD pays off when the samples far outnumber the subspace you can afford to diagonalize.

MPI grid parameters

Parameter What it controls Default
--adet_comm_size Ranks spanning the alpha-determinant dimension 1
--bdet_comm_size Ranks spanning the beta-determinant dimension 1
--task_comm_size Ranks spanning task-level parallelism 1

All ranks diagonalize each batch together, then move to the next batch sequentially. See MPI Decomposition below for the full 4D grid (a fourth, derived dimension — helper — is not set directly).

Checkpointing parameters

Parameter What it controls Default
--checkpoint_path Write ci_strs_a/ci_strs_b/orbital_occupancies/energy to this path as JSON text (rank 0 only), every --checkpoint_frequency iterations. Must be visible under the same path from every rank — --resume_from has no rank-0-reads-then-broadcasts step, every rank opens this path itself unset
--checkpoint_frequency Write every this many iterations, plus always on the last one regardless of alignment 1 (every iteration)
--resume_from Seed a new run's include_configurations/initial_occupancies from a previous --checkpoint_path's last recorded iteration unset

What's preserved is which determinants (ci_strs_a/ci_strs_b, as plain integers — independent of MPI decomposition, since they carry no rank/grid information) and the derived per-orbital occupancy averages. Not the wavefunction amplitudes matrix itself (dim_a × dim_b floats — at large --max_dim this alone would dwarf everything else; neither consuming parameter needs it).

Not a bit-identical continuation: RNG state is fresh in the new process, and every string from the resumed iteration becomes a permanent include for every iteration of the new run — unlike a true single-process continuation, where --sqd_carryover_threshold would keep pruning low-weight determinants each iteration. A resumed run is seeded richer than an actual continuation would have been at that point, not identical to one.

MPI Decomposition

Total MPI ranks must be a multiple of task_comm_size × adet_comm_size × bdet_comm_size — not equal to it. SBD splits the ranks you asked for across those three dimensions and puts whatever remains into a fourth, "helper" dimension, computed as ranks / (task_comm_size × adet_comm_size × bdet_comm_size).

So 8 ranks with --adet_comm_size 2 --bdet_comm_size 2 is valid: the grid is 1 × 2 × 2 and the helper dimension absorbs the remaining factor of 2.

When using more than one rank, specify at least --adet_comm_size. Examples:

Ranks Decomposition Helper
1 default (all = 1) 1
2 --adet_comm_size 2 1
4 --adet_comm_size 2 --bdet_comm_size 2 1
8 --adet_comm_size 2 --bdet_comm_size 2 2
8 --adet_comm_size 2 --bdet_comm_size 2 --task_comm_size 2 1

Expected Results

  • H2O: ground state energy ≈ -76.236 Hartree
  • N2: ground state energy ≈ -109.042 Hartree (with 1e-3 dets)

Backend Selection

Every driver takes --device: cpu (the default), gpu (Thrust, NVIDIA), gpu-omp (OpenMP offload, NVIDIA or AMD), or auto, which takes gpu, then gpu-omp, when a GPU is present, and otherwise cpu. Only the backend asked for is loaded. From Python, pass device= to each call; no sbd.init() is needed:

result_cpu = sbd.tpb_diag(fcidump, adet, bdet, config, device="cpu")
result_gpu = sbd.tpb_diag(fcidump, adet, bdet, config, device="gpu")   # same process: fine

Do not use cpu and gpu-omp in one process. They share an OpenMP runtime, and whichever loads first decides it: if the CPU backend loads first, gpu-omp then runs on the host, with correct energies and the GPU idle. sbd.has_backend_conflict() reports it. cpu and gpu can share a process.

Unrestricted (UHF) FCIDUMPs

TPB and GDB take a UHF FCIDUMP through the same calls and drivers as a restricted one. Each call reads the file's layout from its largest integral index (below) and, for a UHF file, uses the unrestricted build of the device asked for: upstream SBD built with -D_UHF, which every install builds next to the cpu and gpu backends. There is no option to set.

The layout is upstream SBD's own:

  • NORB stays the number of spatial orbitals.
  • Integrals are indexed by interleaved, 1-based spin orbitals: alpha of spatial orbital p (0-based) is 2p+1, beta is 2p+2, so indices run to 2*NORB.
  • One-body lines are h i j 0 0; two-body lines are chemist-notation (ij|kl) with i, j of one spin and k, l of one spin, each unique integral written once.

A restricted file indexes spatial orbitals, up to NORB, so the largest index tells the two apart; sbd.fcidump_format(path) reports 'uhf' or 'restricted'. A file that fits neither is refused: an index above 2*NORB, indices above NORB whose pair mixes spins, or a UHF=/IUHF= header flag, which marks another program's layout.

Making one. ../tools/make_uhf_fcidump.py writes a UHF FCIDUMP with PySCF (which you install): its write_uhf_fcidump takes the spin blocks as pyscf.ao2mo.general gives them. As an example, the NH2 radical in 6-31G, with every alpha and beta string (920,205 determinants) and its exact energy:

python ../tools/make_uhf_fcidump.py nh2 --outdir nh2
python run_sbd_diag.py --fcidump nh2/nh2_uhf.fcidump --symmetrize_spin 0 \
    --adetfile nh2/nh2_alpha.txt --bdetfile nh2/nh2_beta.txt --tolerance 1e-8
# Ground state energy: the e_fci in nh2/reference.json

On NVIDIA GPUs the same command takes --device gpu, under mpirun as for a restricted file.

--symmetrize_spin 0 is needed for an open-shell system, restricted or not: the default reuses the alpha list as the beta one, which holds the wrong number of electrons. The driver says so if it is left out.

Limits:

  • Backends: cpu and gpu (Thrust). gpu-omp does not support UHF files: it has no unrestricted build, and passing one raises an error. sbd.has_uhf(device) reports whether a device has one; every install builds them for cpu and gpu.
  • SQD: the SQD loop (sbd.sbd_solver.solve_sci) is restricted only for now.
  • RDMs: --rdm_output saves rdm1/rdm2 summed over spin, as upstream SBD writes them. The alpha and beta orbitals of a UHF state differ, so the per-spin blocks are the useful ones: sbd.tpb_diag and sbd.gdb_diag return them as one_p_rdm/two_p_rdm.
  • Variance: ../analysis/variance.py detects a UHF file and uses the stochastic method; the exact method is restricted only.
  • Memory: the unrestricted two-body store is 4*NORB^4 doubles, about 32 times the restricted one, and like it is replicated on every rank.

Available Test Data

Paths are relative to this folder, as the drivers' defaults are.

H2O (../../vendor/sbd-upstream/data/h2o/): h2o-1em3 through h2o-1em8 alpha determinant files. N2 (../../vendor/sbd-upstream/data/n2/): 1em3 through 1em7 and 3em4 through 3em7 alpha determinant files.

Smaller thresholds = more determinants = higher accuracy.

Performance Tips

  • CPU: set OMP_NUM_THREADS to the cores per rank, e.g. 8 ranks × 4 threads on 32 cores.
  • GPU: one rank per GPU; each rank takes gpu_id = rank % num_gpus by itself. Use method 0 (matrix-free Davidson).
  • Keep threads on GPU runs too. Building the helper lists runs on the host, so divide the cores among the ranks as for a CPU run (96 cores, 8 ranks: OMP_NUM_THREADS=12); one thread per rank makes that part serial.

See Also