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

python/__init__.py

@@ -727,6 +727,16 @@ def tpb_diag(fcidump, adet, bdet, sbd_data,
         them. ``sort_bitarray`` removes duplicates as well as sorting, so pass
         lists already in that form if you need to label the rows and columns
         yourself.
+
+        Both files use SBD's sign convention, which orders a determinant's
+        creation operators orbital by orbital (alpha 0, beta 0, alpha 1, ...).
+        PySCF, and so qiskit-addon-sqd's ``SCIState``, puts all alpha operators
+        before all beta ones, so some determinants differ in sign: multiply each
+        amplitude by ``(-1) ** n``, where ``n`` counts, over the occupied beta
+        orbitals ``p``, the occupied alpha orbitals above ``p``, before handing
+        it to PySCF. Magnitudes, and so selection by weight, are the same either
+        way. :func:`sbd.sbd_solver.solve_sci` converts the amplitudes in its
+        ``SCIState`` already.
     """
     _ensure_initialized()
     backend = get_backend(device)

python/sbd_solver.py

@@ -297,8 +297,10 @@ def _solve_sci_core(
                 "a different size means the dump and the subspace have diverged."
             )
 
+        amplitudes = flat.reshape(n_a, n_b)
+        _to_alpha_first_signs(amplitudes, solved_strings_a, solved_strings_b, norb)
         sci_state = SCIState(
-            amplitudes=flat.reshape(n_a, n_b),
+            amplitudes=amplitudes,
             ci_strs_a=solved_strings_a,
             ci_strs_b=solved_strings_b,
             norb=norb,
@@ -688,6 +690,41 @@ def _ci_strings_to_sbd_dets(
     return backend.sort_bitarray(dets)
 
 
+def _to_alpha_first_signs(amplitudes, strings_a, strings_b, norb, block_elements=1 << 24):
+    """Convert SBD's amplitudes, in place, to the sign convention SCIState assumes.
+
+    The two conventions order the creation operators of a determinant differently.
+    SBD interleaves the spins orbital by orbital (alpha 0, beta 0, alpha 1, ...);
+    PySCF, and so ``SCIState.rdm()`` and ``SCIState.spin_square()``, put every alpha
+    operator before every beta one. Moving each beta operator past the alpha
+    operators of higher orbitals gives the determinant the phase
+
+        (-1) ** sum over occupied beta orbitals p of (occupied alpha orbitals above p)
+
+    The magnitudes are the same, so energies, SBD's own RDMs, occupancies and any
+    selection by amplitude magnitude are unaffected, which is why this went unnoticed:
+    only what PySCF computes *from the amplitudes* was wrong -- e.g. ``spin_square()``
+    of 0.036 for an exact singlet, and an energy from ``rdm()`` 89 mHa off.
+
+    Args:
+        amplitudes: The ``(len(strings_a), len(strings_b))`` matrix, modified in place.
+        strings_a: Alpha CI strings labelling the rows; bit ``p`` is orbital ``p``.
+        strings_b: Beta CI strings labelling the columns.
+        norb: Number of spatial orbitals.
+        block_elements: Rows are converted in blocks of about this many elements, so
+            the temporary sign matrix stays bounded at large subspaces.
+    """
+    orbitals = np.arange(norb, dtype=np.int64)
+    occupied_b = ((np.asarray(strings_b, dtype=np.int64)[:, None] >> orbitals) & 1).astype(np.int32)
+    occupied_a = (np.asarray(strings_a, dtype=np.int64)[:, None] >> orbitals) & 1
+    # alpha electrons in orbitals strictly above p, kept only as a parity
+    above_a = ((np.cumsum(occupied_a[:, ::-1], axis=1)[:, ::-1] - occupied_a) & 1).astype(np.int32)
+    rows = max(1, block_elements // max(1, len(strings_b)))
+    for start in range(0, len(strings_a), rows):
+        parity = (above_a[start:start + rows] @ occupied_b.T) & 1
+        amplitudes[start:start + rows][parity == 1] *= -1.0
+
+
 def _sbd_dets_to_ci_strings(
     dets: list[list[int]], norb: int, backend,
     bit_length: int = SBD_DEFAULT_BIT_LENGTH,

releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml

@@ -0,0 +1,15 @@
+---
+fixes:
+  - |
+    ``SCIResult.sci_state`` returned by :func:`sbd.sbd_solver.solve_sci` and
+    :func:`sbd.sbd_solver.solve_sci_batch` now holds amplitudes in the sign
+    convention ``qiskit-addon-sqd``'s ``SCIState`` assumes. SBD orders a
+    determinant's creation operators orbital by orbital, PySCF puts all alpha
+    operators before all beta ones, so some determinants carried the opposite
+    sign. Everything ``SCIState`` computes through PySCF from the amplitudes
+    was wrong as a result: ``spin_square()`` returned 0.036 for an exact
+    singlet, and an energy rebuilt from ``rdm()`` was 89 mHa off. Energies,
+    ``SCIResult.rdm1``/``rdm2`` (SBD's own RDMs), orbital occupancies and
+    carryover selection, which uses only amplitude magnitudes, were not
+    affected. The raw wavefunction files SBD writes keep its own convention;
+    the :func:`sbd.tpb_diag` documentation gives the conversion.

test/test_sci_state.py

@@ -0,0 +1,126 @@
+# This code is a Qiskit project.
+#
+# (C) Copyright IBM 2026.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""``SCIResult.sci_state`` from ``solve_sci`` must be the state PySCF would return.
+
+``SCIState`` hands its amplitudes to PySCF for ``rdm()`` and ``spin_square()``, so
+they have to follow PySCF's sign convention, in which every alpha creation operator
+precedes every beta one. SBD orders them orbital by orbital instead, which gives
+some determinants the opposite sign. Energies, SBD's own RDMs (``SCIResult.rdm1`` /
+``rdm2``), occupancies and selection by amplitude magnitude never see the
+difference, which is why only these comparisons catch it.
+
+Each case solves one subspace twice -- with SBD through ``solve_sci``, and with
+qiskit-addon-sqd's PySCF solver -- and compares the states, not just the energies.
+The subspaces are full spaces and seeded random truncations of them: a small
+truncation can happen to need no sign change at all, so a lone small case can pass
+by accident.
+"""
+
+from __future__ import annotations
+
+import itertools
+
+import numpy as np
+import pytest
+
+pyscf = pytest.importorskip("pyscf", reason="pyscf builds the reference Hamiltonians")
+fermion = pytest.importorskip("qiskit_addon_sqd.fermion",
+                              reason="qiskit-addon-sqd provides SCIState and the reference solver")
+
+from pyscf import ao2mo, gto, scf  # noqa: E402
+
+from sbd.sbd_solver import solve_sci  # noqa: E402
+
+SYSTEMS = {
+    # (atom, basis, 2S): a doublet, a singlet, and one with more orbitals
+    "h3-doublet": ("H 0 0 0; H 0 0 0.8; H 0 0 2.2", "sto-3g", 1),
+    "h4-singlet": ("H 0 0 0; H 0 0 0.9; H 0 0 1.9; H 0 0 3.0", "sto-3g", 0),
+    "h6-singlet": ("; ".join(f"H 0 0 {0.95 * i:.2f}" for i in range(6)), "sto-3g", 0),
+}
+
+
+def _hamiltonian(name):
+    atom, basis, spin = SYSTEMS[name]
+    mol = gto.M(atom=atom, basis=basis, spin=spin, verbose=0)
+    mf = (scf.ROHF if spin else scf.RHF)(mol)
+    mf.conv_tol = 1e-12
+    mf.run()
+    c = mf.mo_coeff
+    norb = mol.nao
+    h1 = c.T @ scf.hf.get_hcore(mol) @ c
+    eri = ao2mo.restore(1, ao2mo.kernel(mol, c), norb)
+    return h1, eri, norb, mol.nelec
+
+
+def _strings(norb, n_electrons, keep=None, seed=0):
+    full = np.array(sorted(sum(1 << o for o in occ)
+                           for occ in itertools.combinations(range(norb), n_electrons)))
+    if keep is None or keep >= len(full):
+        return full
+    chosen = np.random.default_rng(seed).choice(len(full), size=keep, replace=False)
+    return np.sort(full[np.union1d(chosen, [0])])        # keep Hartree-Fock
+
+
+CASES = [
+    ("h3-doublet", None),
+    ("h4-singlet", None),
+    ("h6-singlet", None),
+    ("h6-singlet", 12),
+]
+
+
+def _energy_from_rdms(state, h1, eri):
+    return (np.einsum("pr,pr->", state.rdm(rank=1, spin_summed=True), h1)
+            + 0.5 * np.einsum("prqs,prqs->", state.rdm(rank=2, spin_summed=True), eri))
+
+
+def _check_sci_state_matches_pyscf(device_config, name, keep):
+    h1, eri, norb, nelec = _hamiltonian(name)
+    strings = (_strings(norb, nelec[0], keep, seed=1), _strings(norb, nelec[1], keep, seed=2))
+
+    # PySCF's default tolerance leaves the vector good to ~1e-6; converge both tightly
+    # so that the elementwise comparison measures the convention, not convergence.
+    expected = fermion.solve_sci(strings, h1, eri, norb=norb, nelec=nelec, tol=1e-14)
+    got = solve_sci(strings, h1, eri, norb=norb, nelec=nelec, device_config=device_config,
+                    sbd_config={"eps": 1e-12, "max_it": 300})
+    from mpi4py import MPI
+    if MPI.COMM_WORLD.Get_rank() != 0:
+        return                        # the amplitudes are gathered onto rank 0 only
+
+    assert got.energy == pytest.approx(expected.energy, abs=1e-8)
+    assert np.array_equal(got.sci_state.ci_strs_a, expected.sci_state.ci_strs_a)
+    assert np.array_equal(got.sci_state.ci_strs_b, expected.sci_state.ci_strs_b)
+
+    # The same vector up to one global sign -- the comparison that catches a
+    # per-determinant sign, where every magnitude-based one passes.
+    a, b = got.sci_state.amplitudes, expected.sci_state.amplitudes
+    overlap = np.vdot(b, a)
+    assert abs(overlap) == pytest.approx(1.0, abs=1e-7)
+    np.testing.assert_allclose(a * np.sign(overlap), b, atol=1e-6)
+
+    # And what users compute from it.
+    assert got.sci_state.spin_square() == pytest.approx(expected.sci_state.spin_square(), abs=1e-7)
+    assert _energy_from_rdms(got.sci_state, h1, eri) == pytest.approx(expected.energy, abs=1e-7)
+
+
+# One body, two wrappers: pytest-mpi skips unmarked tests under --only-mpi and marked
+# ones without it (see test_sqd_integration.py).
+@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES])
+def test_sci_state_matches_pyscf_standalone(device_config, name, keep):
+    _check_sci_state_matches_pyscf(device_config, name, keep)
+
+
+@pytest.mark.mpi
+@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES])
+def test_sci_state_matches_pyscf_mpi(device_config, name, keep):
+    _check_sci_state_matches_pyscf(device_config, name, keep)

tox.ini

@@ -29,6 +29,25 @@ passenv =
     OMP_NUM_THREADS
     PATH
     LD_LIBRARY_PATH
+# Work around Qiskit/qiskit#17213: with qiskit 2.6.0 installed, importing our
+# extension fails on Linux with "cannot allocate memory in static TLS block".
+#
+# glibc sizes the per-thread TLS block at startup and keeps only a small surplus for
+# libraries dlopen'ed afterwards that need initial-exec TLS. Every Python extension
+# arrives that way, so they all draw on that surplus. qiskit 2.6.0's _accelerate grew
+# from 736 to 1544 bytes of it (matrixmultiply 0.3.11 enabling AVX-512 kernels, which
+# quadruple one thread-local buffer), which leaves too little for the libgomp our
+# extension pulls in. The error names libgomp, but it is whoever asks last that fails.
+#
+# Raising the surplus restores the headroom. It is reserved per thread and only commits
+# the pages actually used, so a generous value costs little. Drop this once 17213 is
+# fixed upstream.
+#
+# Linux-only: this is a glibc mechanism, and macOS links libomp with no equivalent
+# limit. Note glibc silently ignores an unrecognized tunable name, so a typo here looks
+# like a working config.
+setenv =
+    linux: GLIBC_TUNABLES = glibc.rtld.optional_static_tls=65536
 extras =
     test
 commands =
@@ -52,6 +71,9 @@ allowlist_externals =
     mpirun
 setenv =
     SBD_TEST_NPROCS = {env:SBD_TEST_NPROCS:2}
+    # setenv is replaced, not merged, in a child environment, so repeat the
+    # static-TLS headroom [testenv] sets; see the comment there.
+    linux: GLIBC_TUNABLES = glibc.rtld.optional_static_tls=65536
 commands =
     mpirun -n {env:SBD_TEST_NPROCS} pytest --only-mpi {posargs}