diff --git a/python/__init__.py b/python/__init__.py index d0a9b78..e0f1c7e 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -727,6 +727,16 @@ def tpb_diag(fcidump, adet, bdet, sbd_data, them. ``sort_bitarray`` removes duplicates as well as sorting, so pass lists already in that form if you need to label the rows and columns yourself. + + Both files use SBD's sign convention, which orders a determinant's + creation operators orbital by orbital (alpha 0, beta 0, alpha 1, ...). + PySCF, and so qiskit-addon-sqd's ``SCIState``, puts all alpha operators + before all beta ones, so some determinants differ in sign: multiply each + amplitude by ``(-1) ** n``, where ``n`` counts, over the occupied beta + orbitals ``p``, the occupied alpha orbitals above ``p``, before handing + it to PySCF. Magnitudes, and so selection by weight, are the same either + way. :func:`sbd.sbd_solver.solve_sci` converts the amplitudes in its + ``SCIState`` already. """ _ensure_initialized() backend = get_backend(device) diff --git a/python/sbd_solver.py b/python/sbd_solver.py index b346158..59f975f 100644 --- a/python/sbd_solver.py +++ b/python/sbd_solver.py @@ -297,8 +297,10 @@ def _solve_sci_core( "a different size means the dump and the subspace have diverged." ) + amplitudes = flat.reshape(n_a, n_b) + _to_alpha_first_signs(amplitudes, solved_strings_a, solved_strings_b, norb) sci_state = SCIState( - amplitudes=flat.reshape(n_a, n_b), + amplitudes=amplitudes, ci_strs_a=solved_strings_a, ci_strs_b=solved_strings_b, norb=norb, @@ -688,6 +690,41 @@ def _ci_strings_to_sbd_dets( return backend.sort_bitarray(dets) +def _to_alpha_first_signs(amplitudes, strings_a, strings_b, norb, block_elements=1 << 24): + """Convert SBD's amplitudes, in place, to the sign convention SCIState assumes. + + The two conventions order the creation operators of a determinant differently. + SBD interleaves the spins orbital by orbital (alpha 0, beta 0, alpha 1, ...); + PySCF, and so ``SCIState.rdm()`` and ``SCIState.spin_square()``, put every alpha + operator before every beta one. Moving each beta operator past the alpha + operators of higher orbitals gives the determinant the phase + + (-1) ** sum over occupied beta orbitals p of (occupied alpha orbitals above p) + + The magnitudes are the same, so energies, SBD's own RDMs, occupancies and any + selection by amplitude magnitude are unaffected, which is why this went unnoticed: + only what PySCF computes *from the amplitudes* was wrong -- e.g. ``spin_square()`` + of 0.036 for an exact singlet, and an energy from ``rdm()`` 89 mHa off. + + Args: + amplitudes: The ``(len(strings_a), len(strings_b))`` matrix, modified in place. + strings_a: Alpha CI strings labelling the rows; bit ``p`` is orbital ``p``. + strings_b: Beta CI strings labelling the columns. + norb: Number of spatial orbitals. + block_elements: Rows are converted in blocks of about this many elements, so + the temporary sign matrix stays bounded at large subspaces. + """ + orbitals = np.arange(norb, dtype=np.int64) + occupied_b = ((np.asarray(strings_b, dtype=np.int64)[:, None] >> orbitals) & 1).astype(np.int32) + occupied_a = (np.asarray(strings_a, dtype=np.int64)[:, None] >> orbitals) & 1 + # alpha electrons in orbitals strictly above p, kept only as a parity + above_a = ((np.cumsum(occupied_a[:, ::-1], axis=1)[:, ::-1] - occupied_a) & 1).astype(np.int32) + rows = max(1, block_elements // max(1, len(strings_b))) + for start in range(0, len(strings_a), rows): + parity = (above_a[start:start + rows] @ occupied_b.T) & 1 + amplitudes[start:start + rows][parity == 1] *= -1.0 + + def _sbd_dets_to_ci_strings( dets: list[list[int]], norb: int, backend, bit_length: int = SBD_DEFAULT_BIT_LENGTH, diff --git a/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml new file mode 100644 index 0000000..c32ec5c --- /dev/null +++ b/releasenotes/notes/sci-state-sign-convention-58b57b05535e1ca2.yaml @@ -0,0 +1,15 @@ +--- +fixes: + - | + ``SCIResult.sci_state`` returned by :func:`sbd.sbd_solver.solve_sci` and + :func:`sbd.sbd_solver.solve_sci_batch` now holds amplitudes in the sign + convention ``qiskit-addon-sqd``'s ``SCIState`` assumes. SBD orders a + determinant's creation operators orbital by orbital, PySCF puts all alpha + operators before all beta ones, so some determinants carried the opposite + sign. Everything ``SCIState`` computes through PySCF from the amplitudes + was wrong as a result: ``spin_square()`` returned 0.036 for an exact + singlet, and an energy rebuilt from ``rdm()`` was 89 mHa off. Energies, + ``SCIResult.rdm1``/``rdm2`` (SBD's own RDMs), orbital occupancies and + carryover selection, which uses only amplitude magnitudes, were not + affected. The raw wavefunction files SBD writes keep its own convention; + the :func:`sbd.tpb_diag` documentation gives the conversion. diff --git a/test/test_sci_state.py b/test/test_sci_state.py new file mode 100644 index 0000000..1f39aac --- /dev/null +++ b/test/test_sci_state.py @@ -0,0 +1,126 @@ +# This code is a Qiskit project. +# +# (C) Copyright IBM 2026. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +"""``SCIResult.sci_state`` from ``solve_sci`` must be the state PySCF would return. + +``SCIState`` hands its amplitudes to PySCF for ``rdm()`` and ``spin_square()``, so +they have to follow PySCF's sign convention, in which every alpha creation operator +precedes every beta one. SBD orders them orbital by orbital instead, which gives +some determinants the opposite sign. Energies, SBD's own RDMs (``SCIResult.rdm1`` / +``rdm2``), occupancies and selection by amplitude magnitude never see the +difference, which is why only these comparisons catch it. + +Each case solves one subspace twice -- with SBD through ``solve_sci``, and with +qiskit-addon-sqd's PySCF solver -- and compares the states, not just the energies. +The subspaces are full spaces and seeded random truncations of them: a small +truncation can happen to need no sign change at all, so a lone small case can pass +by accident. +""" + +from __future__ import annotations + +import itertools + +import numpy as np +import pytest + +pyscf = pytest.importorskip("pyscf", reason="pyscf builds the reference Hamiltonians") +fermion = pytest.importorskip("qiskit_addon_sqd.fermion", + reason="qiskit-addon-sqd provides SCIState and the reference solver") + +from pyscf import ao2mo, gto, scf # noqa: E402 + +from sbd.sbd_solver import solve_sci # noqa: E402 + +SYSTEMS = { + # (atom, basis, 2S): a doublet, a singlet, and one with more orbitals + "h3-doublet": ("H 0 0 0; H 0 0 0.8; H 0 0 2.2", "sto-3g", 1), + "h4-singlet": ("H 0 0 0; H 0 0 0.9; H 0 0 1.9; H 0 0 3.0", "sto-3g", 0), + "h6-singlet": ("; ".join(f"H 0 0 {0.95 * i:.2f}" for i in range(6)), "sto-3g", 0), +} + + +def _hamiltonian(name): + atom, basis, spin = SYSTEMS[name] + mol = gto.M(atom=atom, basis=basis, spin=spin, verbose=0) + mf = (scf.ROHF if spin else scf.RHF)(mol) + mf.conv_tol = 1e-12 + mf.run() + c = mf.mo_coeff + norb = mol.nao + h1 = c.T @ scf.hf.get_hcore(mol) @ c + eri = ao2mo.restore(1, ao2mo.kernel(mol, c), norb) + return h1, eri, norb, mol.nelec + + +def _strings(norb, n_electrons, keep=None, seed=0): + full = np.array(sorted(sum(1 << o for o in occ) + for occ in itertools.combinations(range(norb), n_electrons))) + if keep is None or keep >= len(full): + return full + chosen = np.random.default_rng(seed).choice(len(full), size=keep, replace=False) + return np.sort(full[np.union1d(chosen, [0])]) # keep Hartree-Fock + + +CASES = [ + ("h3-doublet", None), + ("h4-singlet", None), + ("h6-singlet", None), + ("h6-singlet", 12), +] + + +def _energy_from_rdms(state, h1, eri): + return (np.einsum("pr,pr->", state.rdm(rank=1, spin_summed=True), h1) + + 0.5 * np.einsum("prqs,prqs->", state.rdm(rank=2, spin_summed=True), eri)) + + +def _check_sci_state_matches_pyscf(device_config, name, keep): + h1, eri, norb, nelec = _hamiltonian(name) + strings = (_strings(norb, nelec[0], keep, seed=1), _strings(norb, nelec[1], keep, seed=2)) + + # PySCF's default tolerance leaves the vector good to ~1e-6; converge both tightly + # so that the elementwise comparison measures the convention, not convergence. + expected = fermion.solve_sci(strings, h1, eri, norb=norb, nelec=nelec, tol=1e-14) + got = solve_sci(strings, h1, eri, norb=norb, nelec=nelec, device_config=device_config, + sbd_config={"eps": 1e-12, "max_it": 300}) + from mpi4py import MPI + if MPI.COMM_WORLD.Get_rank() != 0: + return # the amplitudes are gathered onto rank 0 only + + assert got.energy == pytest.approx(expected.energy, abs=1e-8) + assert np.array_equal(got.sci_state.ci_strs_a, expected.sci_state.ci_strs_a) + assert np.array_equal(got.sci_state.ci_strs_b, expected.sci_state.ci_strs_b) + + # The same vector up to one global sign -- the comparison that catches a + # per-determinant sign, where every magnitude-based one passes. + a, b = got.sci_state.amplitudes, expected.sci_state.amplitudes + overlap = np.vdot(b, a) + assert abs(overlap) == pytest.approx(1.0, abs=1e-7) + np.testing.assert_allclose(a * np.sign(overlap), b, atol=1e-6) + + # And what users compute from it. + assert got.sci_state.spin_square() == pytest.approx(expected.sci_state.spin_square(), abs=1e-7) + assert _energy_from_rdms(got.sci_state, h1, eri) == pytest.approx(expected.energy, abs=1e-7) + + +# One body, two wrappers: pytest-mpi skips unmarked tests under --only-mpi and marked +# ones without it (see test_sqd_integration.py). +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_standalone(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep) + + +@pytest.mark.mpi +@pytest.mark.parametrize("name, keep", CASES, ids=[f"{n}-{k or 'full'}" for n, k in CASES]) +def test_sci_state_matches_pyscf_mpi(device_config, name, keep): + _check_sci_state_matches_pyscf(device_config, name, keep)