Fermionic pools¶

Enumerate the fermionic excitation pools that ADAPT-GCiM selects from, apply their anti-Hermitian generators to state vectors without forming matrices, build exact circuits of their exponentials, and read XACC fermionic Hamiltonians into Qiskit Pauli operators. Import from nwqlib.subroutines.fermionic_pool, nwqlib.subroutines.fermionic_circuits and nwqlib.subroutines.xacc. The GCiM guide describes how the ADAPT Method uses the pools.

import numpy as np
from qiskit.quantum_info import Statevector
from nwqlib.subroutines.fermionic_circuits import build_generator_circuit
from nwqlib.subroutines.fermionic_pool import (
    apply_generator_exponential,
    enumerate_spin_adapted_gsd_pool,
)

pool = enumerate_spin_adapted_gsd_pool(2)  # 2 spatial orbitals, 4 qubits
for generator in pool:
    print(generator.family, generator.spatial_indices)
state = np.zeros(16, dtype=complex)
state[0b0011] = 1  # spatial orbital 0 doubly occupied
circuit = build_generator_circuit(pool[0], 0.3)
action = apply_generator_exponential(pool[0], state, 0.3)
print(np.allclose(Statevector(state).evolve(circuit).data, action))

single (0, 1)
double_singlet (0, 0, 0, 1)
double_singlet (0, 0, 1, 1)
double_singlet (0, 1, 1, 1)
True

The circuit of exp(0.3 A) for the first pool member agrees with its matrix-free Taylor action.

Pools and generator actions¶

Fermionic excitation pools for ADAPT-GCiM and their matrix-free actions.

Every pool member is an anti-Hermitian generator A stored as its Pauli sum and, for the fermionic families, also as normal-ordered fermionic terms. ADAPT-GCiM prepares exp(theta A) from it.

Conventions. Spin orbital (g, up) is mode 2*g and (g, down) is mode 2*g + 1 for zero-based spatial orbital g, as in Zheng et al. (2024), arXiv:2312.07691v3, Appendix E.1 (text after Eq. (E3), p. 13). Mode j is qubit j, and the Jordan-Wigner images are those of Eq. (E8), a_j = (X_j + i Y_j)/2 prod_{k<j} Z_k. In a Qiskit Pauli label qubit 0 is the rightmost character, so qubit j is label position num_qubits - 1 - j. Equation labels E1-E8 below refer to that appendix.

Default pool (enumerate_spin_adapted_gsd_pool). The generalized spin-adapted singles and doubles of Appendix E.1, Eqs. (E1)-(E3). Singles come first, over spatial pairs p < q. Doubles follow over spatial pairs (p, q) and (r, s) with p <= q, r <= s and rank(p, q) <= rank(r, s) in lexicographic order. For each index quadruple the triplet of Eq. (E3) precedes the singlet of Eq. (E2), and an operator that vanishes identically is skipped. Pool indices follow this order. The singlet and triplet couplings are those of Magoulas and Evangelista, arXiv:2511.13485v2, Sec. III, Eqs. (6) and (7).

Orientation and sign. The builders write each generator as B - B^dagger, normal-order it with the canonical anticommutation relations and scale it to unit Euclidean norm over the normal-ordered coefficients. A single (p, q) is a_{p up}^dagger a_{q up} + a_{p dn}^dagger a_{q dn} - h.c., Eq. (E1). The double base terms are interleaved products a_r^dagger a_p a_s^dagger a_q with spin labels. Since a_p a_s^dagger = delta_ps - a_s^dagger a_p, the base operator is B = -X plus one-body terms proportional to delta_ps, where X = a_r^dagger a_s^dagger a_p a_q (spin labels implied). The generator is then B - B^dagger = X^dagger - X with X^dagger = a_q^dagger a_p^dagger a_s a_r = a_p^dagger a_q^dagger a_r a_s, which creates in (p, q) and annihilates in (r, s). So the library double (p, q, r, s) is Eq. (E2) or (E3) with the same index order and a positive normalization factor. Reading only the first creation index of a base term gives the opposite sign. The one-body terms never enter the pool, because p == s with the ordering above forces p == q == r == s, where B is Hermitian and the operator is skipped. This agreement holds for every enumerated double, including repeated indices. The test test_generalized_pool_has_the_signed_operators_of_zheng_table_5 checks all operators of Table V (Table 5 in the npj version, doi:10.1038/s41534-024-00916-8). The q < p condition of Eq. (E4) belongs to the elementary circuit decomposition and does not restrict this enumeration.

Reversing an orientation replaces A by -A. At a fixed ADAPT angle this changes the product states and the adaptive trajectory. Flipping the sign of a molecular orbital negates every generator whose ladder operators act on that orbital an odd number of times, which is why the chemistry builder fixes the sign of every orbital.

UCCSD-SD pool (enumerate_uccsd_sd_pool). A smaller, reference-dependent comparison pool. Occupied and virtual spin orbitals are read from a reference occupation string, indexed by spin orbital as above. Singles a_a^dagger a_i - h.c. run over occupied i (outer loop, ascending) and virtual a (inner loop, ascending) with equal spin. Doubles a_a^dagger a_b^dagger a_i a_j - h.c. run over occupied pairs i < j (outer) and virtual pairs a < b (inner) with equal total S_z. These are the excitation operators of Romero et al., arXiv:1701.02691v2, Eqs. (8)-(9), in the unitary form of Eq. (17). When occupied orbitals precede virtual ones they are also Zheng et al., arXiv:2312.07691v3, Eqs. (E4)-(E5). Members are normalized like the default pool. spin_orbital_indices holds (i, a) or (i, j, a, b). The occupied-to-virtual restriction with a fixed angle need not span the target sector, so an adaptive run on this pool can stop above the sector minimum.

QEB-SD and OVP-CEO pools. The QEB-SD pool uses the same occupied/virtual enumeration with qubit excitations, which are Jordan-Wigner excitations without parity strings (Yordanov et al., arXiv:2011.10540v2, Eqs. (17)-(18)). The OVP-CEO pool holds sums and differences of qubit excitations on one spin-orbital set, each with a single rotation angle (Ramôa et al., arXiv:2407.08696v3, Sec. II.B.3, Eqs. (23)-(24)). ADAPT-GCiM assigns one angle to each selected operator, so the multi-parameter MVP-CEOs of that paper are not represented.

FermionicGenerator ¶

FermionicGenerator(*, family: PoolFamily, spatial_indices: tuple[int, ...], pool_index: int, num_qubits: int, fermion_terms: tuple[tuple[complex, NormalTerm], ...], pauli_terms: tuple[tuple[str, complex], ...], spin_orbital_indices: tuple[int, ...] | None = None)

One anti-Hermitian ADAPT-GCiM pool generator A, stored as its Pauli sum and, for the fermionic families, its normal-ordered terms.

The enumerate_*_pool functions return tuples of them in pool-index order. The fields below are read-only.

Attributes:

family (PoolFamily) –

Generator family, which fixes its excitation pattern and orientation (see the module text above).
spatial_indices (tuple[int, ...]) –

Ordered spatial-orbital indices defining the generator, (p, q) for singles and (p, q, r, s) for doubles.
pool_index (int) –

Position in the enumerated pool.
num_qubits (int) –

Spin-orbital register width of the Pauli representation.
fermion_terms (tuple[tuple[complex, NormalTerm], ...]) –

(coefficient, operators) pairs of the normal-ordered generator. Each operator is (mode, action) with action 1 for creation and 0 for annihilation, and the rightmost operator acts first. Empty for the qubit-excitation families.
pauli_terms (tuple[tuple[str, complex], ...]) –

(label, coefficient) pairs of the Pauli sum of A. Pauli words are Hermitian and A is anti-Hermitian, so every coefficient is purely imaginary.
spin_orbital_indices (tuple[int, ...] | None) –

Spin orbitals of a reference-dependent member. UCCSD and QEB members and CEO singles store the occupied, then the virtual indices, ascending within each group. Opposite-spin CEO doubles store (occupied alpha, occupied beta, virtual alpha, virtual beta), and same-spin CEO doubles store their four modes sorted. None for the spin-adapted families, whose spatial indices identify them.

spin_orbital ¶

spin_orbital(spatial_orbital: int, spin: Spin) -> int

Return the spin-orbital (mode) index 2*g for spin up and 2*g + 1 for spin down of spatial orbital g.

This is the convention of Zheng et al., arXiv:2312.07691v3, Appendix E.1.

Parameters:

spatial_orbital (int) –

Nonnegative spatial-orbital index g.
spin (str) –

"up" or "down".

Returns:

mode ( int ) –

The spin-orbital index, which is also the qubit index.

enumerate_spin_adapted_gsd_pool ¶

enumerate_spin_adapted_gsd_pool(n_spatial_orbitals: int) -> tuple[FermionicGenerator, ...]

Enumerate the generalized spin-adapted singles and doubles of Zheng et al., arXiv:2312.07691v3, Eqs. (E1)-(E3).

The order and orientation are described in the module docstring. With P = n (n + 1) / 2 spatial pairs for n orbitals, the pool has at most P singles and P (P + 1) doubles, because each pair of pairs contributes at most one triplet and one singlet. Operators that vanish identically are skipped, so the actual count is smaller.

Parameters:

n_spatial_orbitals (int) –

Number of spatial orbitals n. The generators act on 2 n interleaved spin-orbital qubits.

Returns:

tuple[FermionicGenerator, ...] –

The generators in pool-index order.

enumerate_uccsd_sd_pool ¶

enumerate_uccsd_sd_pool(reference_occupations) -> tuple[FermionicGenerator, ...]

Enumerate the occupied-to-virtual UCCSD singles and doubles of a reference.

This smaller, reference-dependent pool is a comparison baseline beside enumerate_spin_adapted_gsd_pool. The module text gives the ordering, index meaning and spin selection, and the source equations (Romero et al., arXiv:1701.02691v2, Eqs. (8)-(9) and (17)).

Parameters:

reference_occupations (str | Sequence[int]) –

Nonempty 0/1 string or sequence whose entry j is the occupation of spin orbital j, in the interleaved order of spin_orbital. It fixes the occupied and virtual spin orbitals.

Returns:

pool ( tuple[FermionicGenerator, ...] ) –

The generators in pool-index order.

enumerate_qeb_sd_pool ¶

enumerate_qeb_sd_pool(reference_occupations) -> tuple[FermionicGenerator, ...]

Enumerate the reference-dependent qubit-excitation singles and doubles (QEB-SD).

The occupied/virtual enumeration mirrors enumerate_uccsd_sd_pool, but the operators are the qubit excitations of Yordanov et al., arXiv:2011.10540v2, Eqs. (17)-(18), p. 6: Q_a^dagger Q_i - Q_i^dagger Q_a and Q_a^dagger Q_b^dagger Q_i Q_j - Q_i^dagger Q_j^dagger Q_a Q_b, with Q_j = (X_j + i Y_j)/2 and no Jordan-Wigner parity string. Eqs. (21)-(22) of the same paper give their exponentials as Pauli rotations.

enumerate_ceo_ovp_pool ¶

enumerate_ceo_ovp_pool(reference_occupations) -> tuple[FermionicGenerator, ...]

Enumerate the one-variational-parameter coupled exchange operators (OVP-CEO) of a reference.

MVP-CEOs are deliberately not represented: they need per-component parameters inside one selected generator, while ADAPT-GCiM records one rotation parameter per selected operator.

The pool holds the single qubit excitations, then for each opposite-spin quadruple the plus and minus combinations of the direct and crossed double qubit excitations (Ramôa et al., arXiv:2407.08696v3, Sec. II.B.3, Eqs. (23)-(24), p. 8). For each same-spin quadruple it holds the sums and differences of every pair of its three qubit excitations, the six OVP-CEOs described on p. 9.

apply_generator ¶

apply_generator(generator: FermionicGenerator, state: ndarray) -> ndarray

Return A @ state for a pool generator A, from its Pauli sum, without forming a matrix.

The cost of this matrix-free action is given on the fermionic pools page.

Parameters:

generator (FermionicGenerator) –

The generator A.
state (ndarray) –

Complex vector of length 2**generator.num_qubits.

Returns:

result ( ndarray ) –

A new vector A @ state.

apply_generator_exponential ¶

apply_generator_exponential(generator: FermionicGenerator, state: ndarray, theta: float) -> ndarray

Return exp(theta * A) @ state for a pool generator A, by scaled Taylor steps without forming a matrix.

The state is not normalized, and a new vector of the same length is returned. The step count is max(1, ceil(2 |theta| sum_k |c_k|)) over the Pauli coefficients c_k of A (0 when that product is zero), so each step h has ||h A|| <= 1/2. Each step is a Taylor polynomial of degree 18, whose remainder is below exp(1/2) (1/2)**19 / 19! < 2.6e-23 per step. The work is 18 matrix-free Pauli-sum actions per step. ADAPT's classical work counts use the same step count.

Parameters:

generator (FermionicGenerator) –

Anti-Hermitian generator A.
state (ndarray) –

Complex vector of length 2**generator.num_qubits.
theta (float) –

Rotation angle in radians.

Returns:

ndarray –

exp(theta * A) @ state up to the Taylor remainder above and
ndarray –

rounding errors.

apply_spin_squared ¶

apply_spin_squared(state: ndarray, *, num_qubits: int) -> ndarray

Return S^2 @ state for the total spin, without forming a matrix, in the same Jordan-Wigner convention.

Uses S^2 = S_- S_+ + S_z (S_z + 1) with S_+ = sum_p a_{p up}^dagger a_{p down} and S_- = S_+^dagger. The ladder product is a two-factor matrix-free action, and S_z is diagonal in the interleaved occupation basis.

Parameters:

state (array_like) –

The 2**num_qubits amplitudes in the interleaved occupation basis, qubit 2p for spatial orbital p with spin up and qubit 2p + 1 with spin down. It is converted to complex.
num_qubits (int) –

Number of spin orbitals, even.

Returns:

vector ( ndarray ) –

S^2 @ state, complex, of length 2**num_qubits.

Raises:

ValueError –

If num_qubits is odd, if the state length is not 2**num_qubits, or if the action needs more than 10 GB (decimal, 10_000_000_000 bytes) of known arrays or more than 1,000,000,000 scalar products. These two limits are fixed and checked before the action.

Cost of the matrix-free Pauli action¶

The Pauli action of apply_generator groups terms by their X/Y flip mask, and NumPy population counts give the Jordan–Wigner parity signs. A group sums its coefficient phases once per coordinate and applies the resulting diagonal to every supplied state column, and the output has the same coordinates as the input. For L Pauli terms, G distinct flip masks and N=2**n coordinates:

The number of state gathers depends on the number of distinct flip masks, while phase formation still visits every Pauli term at every coordinate, so the arithmetic costs O(LN + GN) for one state.
The action keeps one output vector of N complex values and tiles of at most 1,024 coordinates. A required complex input conversion is an additional state-sized allocation. It constructs no dense operator or stored L * 2**n table.

XACC Hamiltonian input¶

parse_xacc and read_xacc accept two-body XACC input (at most four ladder operators per line) and return a SparsePauliOp with its complete identity contribution. The parse_xacc entry below maps 2 a_0^dagger a_0 - 3 I to -2 I - Z, and read_xacc(path, num_modes=1) reads the same text from a local file.

Input format:

Each line's coefficient multiplies the literal ordered operator product. Creation is ^, and annihilation has no suffix. Blank lines, whitespace, complex coefficients, scientific notation and an optional trailing + are supported.
Both entry points strip at most one leading UTF-8 byte-order mark (BOM). Two leading BOMs leave the second in the input and raise ValueError on line 1.
No extra one-half or one-quarter prefactors, symmetry completion or Hermitian projection are applied.

Coefficients and the cutoff:

The default coefficient_cutoff=0.0 keeps every nonzero accumulated coefficient. Normal ordering coalesces equivalent terms before the Jordan–Wigner expansion. Both cross-input normal-order sums and cross-monomial Pauli sums use real and imaginary math.fsum, to keep small residues through large cancellations. Coefficients and individual bounded products use float64, and the result is not exact real arithmetic.
Overflow or a nonfinite accumulated component raises ValueError. fsum also rejects intermediate overflow even if a later term could cancel it.
To approximate the mapped operator deliberately, supply an absolute coefficient_cutoff to either entry point, for example parse_xacc(text, num_modes=24, coefficient_cutoff=1e-12). It removes final Pauli coefficients with magnitude at most the cutoff, after all contributions have been combined, and does not discard individual input terms or normal-order terms.
The removed coefficient magnitudes sum to an upper bound on the operator-norm change, because each Pauli has norm one. The parser does not compute or certify a downstream energy-error bound. For an LCU, the kept label count L sets the SELECT index width ceil(log2(L)). Deliberate pruning can reduce this width and circuit work, and small coefficients are not automatically scientifically irrelevant.

Modes and width:

File mode j maps directly to qubit j. Unlike the spin-adapted pool's interleaved convention, XACC has no implied spin ordering, and blocked and interleaved indices remain as supplied. Electron count, spin ordering and reference occupations must be supplied separately by the code that uses the operator.
The required positive num_modes is both the output width and the caller's bound on the input: every written mode, including those of zero-coefficient terms, must be smaller. The width is never inferred from input indices, and invalid indices fail before the Jordan–Wigner expansion. Valid wide models used only for resource counts are not subject to local simulation limits.
Omitting the required num_modes keyword raises Python's ordinary TypeError. An invalid supplied bound or malformed XACC syntax raises ValueError.
The mapped identity includes contributions from number operators and quartic terms, so the file's scalar alone is not the qubit identity shift.

Cost:

Parsing streams lines and keeps coefficient buckets for stable summation: O(R) coefficient storage for R raw input terms. Packed real and imaginary parts use 16 payload bytes per contribution, plus array capacity and per-key overhead, while keeping fsum accuracy.
These buckets are released before mapping, which costs O(T n) time and storage for T coalesced terms on n modes, with at most 16 branches per term. Pauli coefficient buckets keep at most 16 T contributions during mapping and are released before Qiskit allocates the returned Pauli arrays. No dense matrix, statevector or backend is used.

Bounded XACC fermionic input, mapped literally to Jordan–Wigner Paulis.

parse_xacc ¶

parse_xacc(text: str, *, num_modes: int, coefficient_cutoff: float = 0.0) -> SparsePauliOp

Map XACC text to a qubit Hamiltonian, keeping every nonzero term by default.

Each nonblank line is (real,imag) p^ q ... +; the trailing + is optional. Operators multiply in their literal written order, ^ denotes creation, and a line without operators is a scalar. Decimal and scientific notation are supported. At most four ladder operators per term are supported, keeping normal ordering and JW expansion bounded per term.

Coefficients already include all prefactors. No symmetry completion, Hermitian projection, spin reordering, or extra factors are applied. Cross-input sums use real/imaginary math.fsum at both normal-order and Pauli coalescing stages. Only exactly zero accumulated coefficients are removed by default. An explicit nonnegative coefficient_cutoff removes final Pauli coefficients of magnitude at most that absolute cutoff, after stable coalescing; it never prunes individual input or normal-order terms. Individual coefficients and bounded products use float64 arithmetic. The identity includes all mapped contributions, not only the file scalar.

Mode j maps to qubit j (rightmost Pauli character for mode zero). The required positive num_modes bounds every written index before JW expansion and supplies padding. It is a caller-owned input bound, not a local-simulation limit. UTF-8 text may have an initial byte-order mark. Electron count, spin ordering, and reference occupations are separate caller metadata: none can be inferred from this syntax.

Parameters:

text (str) –

XACC text.
num_modes (int) –

Positive number of modes, which is the output width. Every written mode index must be smaller.
coefficient_cutoff (float, default: 0.0 ) –

Default 0.0. Nonnegative absolute cutoff on the final Pauli coefficients.

Returns:

hamiltonian ( SparsePauliOp ) –

The Jordan-Wigner image on num_modes qubits, including its complete identity term.

Raises:

TypeError –

If num_modes is omitted or text is not a string.
ValueError –

With a line number, for malformed or nonfinite input, an invalid num_modes or a mode index outside it.

Examples:

2 a_0^dagger a_0 - 3 I maps to -2 I - Z, because a_0^dagger a_0 = (I - Z)/2:

>>> from nwqlib.subroutines.xacc import parse_xacc
>>> parse_xacc("(2,0)0^ 0 +\n(-3,0)", num_modes=1).to_list()
[('I', (-2+0j)), ('Z', (-1+0j))]

read_xacc ¶

read_xacc(path: str | PathLike[str], *, num_modes: int, coefficient_cutoff: float = 0.0) -> SparsePauliOp

Read a UTF-8 XACC file and map it to a qubit Hamiltonian, as parse_xacc does for text.

Lines stream into coefficient buckets for stable normal-order summation, then the coalesced polynomial is JW mapped with stable Pauli summation. Buckets keep O(R) coefficients for R input terms and O(T) coefficients for T coalesced terms (bounded degree four); no statevector or dense operator is constructed.

Compact generator circuits¶

build_generator_circuit(generator, theta, ...) constructs exp(theta*A) for a supported anti-Hermitian generator as a QuantumCircuit. An outer control occupies qubit zero and shifts the system modes by one. For an inverse, the builder evaluates the complete parameter map at -theta, including the even angle functions, to produce exp(-theta*A) directly. This avoids synthesizing the controls again through a generic circuit inverse.

The builder supports commuting Pauli generators and the default spin-adapted GSD index classes:

Pair/split and four-distinct doubles use the closed-form excitation factors of Magoulas and Evangelista, arXiv:2511.13485v2. Pair/split doubles follow Sec. V, Eqs. (15), (21), (22) and (25)–(27). Four-distinct doubles follow Sec. VI, Tables I and II, with the generator mapping of Supplemental Tables SI and SII. The library angle is divided by sqrt(2) to give the paper angle, because the paper operators have coefficient norm sqrt(2) in the library's normalization.
Shared-index doubles use NWQLib's own exact block synthesis instead of those closed forms. They use at most a 64 by 64 real active-mode generator and matrix exponentials of dimension at most five.
Fermionic routing includes the signs on occupied spectator modes and cancels when the outer control is zero.
Noncommuting custom records must match a supported default family. The builder does not synthesize them through amplitudes or a generic dense unitary.

These are exact factorization formulas evaluated in floating-point arithmetic. Shared-index two-level synthesis and occupation-conditioned four-distinct rotations can require substantially more gates than small amplitude preparations. Construction is bounded in active support, and no claim of optimal gate count is made. ADAPT composes the generator circuits it chose with its initial preparation and ordered parameter values.

build_generator_circuit ¶

build_generator_circuit(generator: FermionicGenerator, theta: float, *, controlled: bool = False, inverse: bool = False, _plan: _GeneratorCircuitPlan | None = None, max_bytes: int = DEFAULT_INPUT_BYTES) -> QuantumCircuit

Build an exact circuit of exp(theta*A) for one pool generator A, optionally controlled or inverted.

Inversion evaluates the complete parameter map at -theta, including the even angle functions such as alpha_5 of Magoulas and Evangelista, arXiv:2511.13485v2, Eq. (27), because exp(theta*A)^dagger = exp(-theta*A) for anti-Hermitian A. The multi-controlled gates are therefore never passed through a generic circuit inverse, which test_compact_inverse_evaluates_full_map_without_reinverting_controls checks. Custom commuting Pauli sums are supported. A noncommuting custom operator has no exact compact decomposition here and is rejected. Construction never expands a problem state.

For the commuting route A = sum_k i c_k P_k with real c_k, so exp(theta*A) = prod_k exp(i theta c_k P_k). apply_pauli_rotation(P, phi) implements exp(-i phi P / 2), hence phi = -2 theta c_k. An identity term gives the phase exp(i theta c_k). It is a global phase without control and a phase gate on the control qubit under control, where the phase is physical.

Parameters:

generator (FermionicGenerator) –

Anti-Hermitian pool generator of a supported family, or a custom commuting Pauli sum.
theta (float) –

Finite rotation angle in radians. On the shared_index route, theta times each occupation block must have a 1-norm of at most 2**37, the limit of the matrix exponential.
controlled (bool, default: False ) –

Default False. Add one control qubit. The control is qubit 0 and system qubit j becomes qubit j + 1.
inverse (bool, default: False ) –

Default False. Build exp(-theta*A).
_plan (_GeneratorCircuitPlan | None, default: None ) –

Internal reuse of a precomputed plan across angles. Leave it unset.
max_bytes (int, default: DEFAULT_INPUT_BYTES ) –

Default 10 GB (decimal, 10_000_000_000 bytes). Byte limit for choosing the construction.

Returns:

circuit ( QuantumCircuit ) –

Circuit on generator.num_qubits + int(controlled) qubits.