From 18952cf87adb156381dc175eee445d5fc4edc0d9 Mon Sep 17 00:00:00 2001 From: Max Rossmannek Date: Fri, 4 Sep 2026 11:31:15 +0200 Subject: [PATCH] refactor(circuit)!: build UCC from an ffsim UCCSD operator Mirrors the UCJ change: the gate reimplemented the amplitude conventions and the parameter-vector packing alongside the part that is actually this package's concern, namely turning the amplitudes into a mapper-agnostic fermionic circuit. The first half duplicated ffsim, so this narrows the gate to the second half and takes the ffsim operator directly. Verified rather than assumed: the removed `num_parameters`/`from_parameters` agreed with ffsim's `n_params`/`from_parameters` (14 parameters at norb=4, nocc=2, with identical t1/t2), and a restricted operator applied through this gate reproduces ffsim's own state vector exactly on coupled-cluster amplitudes. Unlike UCJ, this drops capability ffsim cannot currently replace, which is why it is a separate commit: `from_t_amplitudes`, the `spinless` variant and the opt-in `antisymmetric` parameterization all go. Users needing those can build an Evolution over their own cluster operator, which also hands them control over the Trotter ordering of its terms. Note that ffsim ships no Qiskit gate for UCCSD at all, so this gate remains the only route from one of its UCCSD operators to a circuit -- under Jordan-Wigner or any other encoding. What remains is the part with no ffsim equivalent: the cluster generator and its conjugate-paired grouping, which keeps every factor of the product formula Hermitian and hence unitary. A term-by-term split would not even preserve the norm, so the grouping tests are kept verbatim; the two that previously reached this machinery only through the spinless variant are ported to the restricted one, which exercises the same paths on a spinful register. One pre-existing convention difference is now documented rather than changed: the cluster operator only ever sees the part of a same-spin t2 that is symmetric under the simultaneous exchange t2[i,j,a,b] = t2[j,i,b,a], while ffsim reads the raw tensor. Coupled-cluster amplitudes always carry that symmetry (verified against PySCF), so the two agree exactly on any physical input; only a hand-built asymmetric tensor can tell them apart. As with UCJ, the class docstring's example now needs ffsim, so it is gated with `.. skip: start if(not HAS_FFSIM)` (the SkipParser this relies on is registered in python/conftest.py by the preceding commit). Without the guard the example fails collection on Windows, where ffsim cannot be installed. Refs #318 Co-Authored-By: Claude Opus 5 --- python/qiskit_fermions/circuit/library/ucc.py | 859 +++--------------- .../ucc-from-ffsim-op-88e7c005f5390543.yaml | 35 + tests/python/circuit/library/test_ucc.py | 499 ++-------- .../circuit/library/test_ucc_apply_unitary.py | 136 +-- 4 files changed, 273 insertions(+), 1256 deletions(-) create mode 100644 releasenotes/notes/ucc-from-ffsim-op-88e7c005f5390543.yaml diff --git a/python/qiskit_fermions/circuit/library/ucc.py b/python/qiskit_fermions/circuit/library/ucc.py index 26de818d7..5e9182fd7 100644 --- a/python/qiskit_fermions/circuit/library/ucc.py +++ b/python/qiskit_fermions/circuit/library/ucc.py @@ -10,453 +10,169 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Unitary coupled cluster (UCC) ansatz gate.""" +"""Unitary coupled-cluster (UCC) ansatz gate.""" from __future__ import annotations import itertools -import sys -from enum import Enum -from typing import TYPE_CHECKING, cast +from collections.abc import Sequence +from typing import TYPE_CHECKING, Any, cast import numpy as np from qiskit_fermions._lib.operators.fermion_operator import FermionOperator from qiskit_fermions.operators.fermion_action import ann, cre +from qiskit_fermions.utils.optionals import HAS_FFSIM from .. import FermionicGate from .evolution import Evolution - -if sys.version_info >= (3, 11): - from typing import Self -else: - from typing_extensions import Self +from .orbital_rotation import OrbitalRotation if TYPE_CHECKING: - from collections.abc import Sequence - from qiskit_fermions.circuit import FermionicCircuit class UCC(FermionicGate): - r"""Implements the unitary coupled cluster (UCC) ansatz. + r"""Implements the unitary coupled-cluster singles and doubles (UCCSD) ansatz. - A unitary coupled cluster operator has the form + A unitary coupled-cluster operator has the form .. math:: - e^{T - T^\dagger} - - where :math:`T = T_1 + T_2` is the cluster operator built from the single and double fermionic - excitations, parameterized by the :math:`t_1` and :math:`t_2` amplitudes. Since - :math:`T - T^\dagger` is anti-Hermitian, its exponential is unitary. - - This gate supports three spin variants (see :class:`.UCC.Variant`), selected explicitly by the - ``variant`` argument and validated against the shapes of the supplied amplitudes (mirroring - ffsim's :external:class:`~ffsim.UCCSDOpRestrictedReal` and - :external:class:`~ffsim.UCCSDOpUnrestrictedReal`): - - - **restricted** -- a single spin-summed amplitude pair. The cluster operator is - - .. math:: - - \begin{align} - T_1 &= \sum_{ia} t_{ia}\left( - a^\dagger_{a\alpha} a_{i\alpha} + a^\dagger_{a\beta} a_{i\beta}\right), \\ - T_2 &= \sum_{ijab} t_{ijab}\left[ - \frac12\left( - a^\dagger_{a\alpha} a^\dagger_{b\alpha} a_{j\alpha} a_{i\alpha} - + a^\dagger_{a\beta} a^\dagger_{b\beta} a_{j\beta} a_{i\beta}\right) - + a^\dagger_{a\alpha} a^\dagger_{b\beta} a_{j\beta} a_{i\alpha}\right], - \end{align} - - with ``t1`` of shape ``(nocc, nvrt)`` and ``t2`` of shape ``(nocc, nocc, nvrt, nvrt)``. Acts - on ``2 * norb`` block-spin modes. - - **unrestricted** -- independent per-spin amplitudes. The cluster operator is - - .. math:: - - \begin{align} - T_1 &= \sum_{ia} t^{(\alpha)}_{ia} a^\dagger_{a\alpha} a_{i\alpha} - + \sum_{IA} t^{(\beta)}_{IA} a^\dagger_{A\beta} a_{I\beta}, \\ - T_2 &= \frac14 \sum_{ijab} t^{(\alpha\alpha)}_{ijab} - a^\dagger_{a\alpha} a^\dagger_{b\alpha} a_{j\alpha} a_{i\alpha} - + \frac14 \sum_{IJAB} t^{(\beta\beta)}_{IJAB} - a^\dagger_{A\beta} a^\dagger_{B\beta} a_{J\beta} a_{I\beta} - + \sum_{iJaB} t^{(\alpha\beta)}_{iJaB} - a^\dagger_{a\alpha} a^\dagger_{B\beta} a_{J\beta} a_{i\alpha}, - \end{align} - - with ``t1`` a pair ``(t1a, t1b)`` and ``t2`` a triple ``(t2aa, t2ab, t2bb)``. Acts on - ``2 * norb`` block-spin modes. Note that the occupied/virtual split is resolved **per spin - sector**, so the two sectors may have different numbers of occupied orbitals. - - **spinless** -- a single register of ``norb`` spinless modes, - - .. math:: - - T_1 = \sum_{ia} t_{ia}\, a^\dagger_a a_i, \qquad - T_2 = \frac14 \sum_{ijab} t_{ijab}\, - a^\dagger_a a^\dagger_b a_j a_i, - - with the same amplitude shapes as the ``"restricted"`` variant. Acts on ``norb`` modes. - - In every variant the occupied orbitals are ordered before the virtual ones, so orbital - :math:`i < n_\text{occ}` is occupied and orbital :math:`n_\text{occ} + a` is virtual. + e^{T - T^\dagger}, \qquad + T = \sum_{ia} t^a_i\, a^\dagger_a a_i + + \sum_{ijab} t^{ab}_{ij}\, a^\dagger_a a^\dagger_b a_j a_i, + + with :math:`i, j` occupied and :math:`a, b` virtual orbitals, and the amplitudes + :math:`t_1 = t^a_i` and :math:`t_2 = t^{ab}_{ij}` supplied by the operator this gate wraps. + + The operator itself is built by `ffsim `__, and this + gate turns it into a :class:`.FermionicCircuit`. That division of labor is deliberate: ffsim owns + the ansatz math (the amplitude conventions and the parameter-vector packing that a variational + optimizer drives), while this gate expresses the result as fermionic modes so that the transpiler + can lower it through *any* fermion-to-qubit encoding. ffsim ships no Qiskit gate for UCCSD at all, + so this is the only route from one of its UCCSD operators to a circuit; see the + :ref:`ffsim guide `. + + Accepts any of ffsim's four UCCSD operators, whose type fixes the spin variant and the number of + modes this gate acts on. All four act on ``2 * norb`` block-spin modes (mode ``p`` is alpha + orbital ``p``, mode ``norb + p`` is beta orbital ``p``), with the occupied orbitals ordered before + the virtual ones: + + - :external:class:`~ffsim.UCCSDOpRestrictedReal` and + :external:class:`~ffsim.UCCSDOpRestricted` share one spatial parametrization between both spin + sectors: ``t1`` has shape ``(nocc, nvrt)`` and ``t2`` has shape ``(nocc, nocc, nvrt, nvrt)``. + - :external:class:`~ffsim.UCCSDOpUnrestrictedReal` and + :external:class:`~ffsim.UCCSDOpUnrestricted` parameterize the spin sectors independently: + ``t1`` is a pair ``(t1a, t1b)`` and ``t2`` a triple ``(t2aa, t2ab, t2bb)``. .. note:: - Unlike :class:`.UCJ`, this ansatz carries no final orbital rotation: its :math:`t_1` - amplitudes already provide the single excitations, so a trailing rotation would be redundant - freedom. Append an :class:`.OrbitalRotation` explicitly if you want one. + The cluster operator only ever sees the part of a same-spin :math:`t_2` block that is + symmetric under the *simultaneous* exchange :math:`t_2[i,j,a,b] = t_2[j,i,b,a]`, because the + underlying excitation :math:`a^\dagger_a a^\dagger_b a_j a_i` is invariant under relabeling + the pairs :math:`(i,a) \leftrightarrow (j,b)`. Coupled-cluster amplitudes (from PySCF, or + ffsim's own :external:func:`~ffsim.uccsd_generator_restricted`) always carry that symmetry, so + this gate and ffsim agree exactly on them. A hand-built ``t2`` without it describes the same + ansatz here as its symmetrized counterpart, whereas ffsim reads the raw tensor. .. note:: - By default only the symmetry the cluster operator actually enforces is imposed on the - same-spin :math:`t_2` blocks. The stricter antisymmetry of the standard coupled-cluster - convention is available opt-in via the ``antisymmetric`` flag (see :attr:`.antisymmetric`), - which both validates supplied amplitudes and shrinks the parameter vector accordingly. It is - not supported for the ``"restricted"`` variant, whose single :math:`t_2` also carries the - cross-spin amplitudes. + Unlike :class:`.UCJ`, this ansatz carries no final orbital rotation of its own: its + :math:`t_1` amplitudes already provide the single excitations, so a trailing rotation would be + redundant freedom. ffsim's operators do expose a ``final_orbital_rotation``; when one is set, + this gate appends it as a closing :class:`.OrbitalRotation`. .. note:: Because the individual excitation terms of :math:`T - T^\dagger` do **not** commute, the circuit :meth:`~qiskit.circuit.Gate.definition` this gate produces is a *first-order product - formula* (Trotter) approximation of the exponential, not an exact decomposition -- the usual - situation for UCC ansatz circuits. The state-vector simulation path + formula* (Trotter) approximation of the exponential, not an exact decomposition (the usual + situation for UCC ansatz circuits). The state-vector simulation path (:meth:`_apply_unitary_placed_`), by contrast, applies the exponential *exactly* via ``scipy``'s ``expm_multiply``. Consequently the simulated gate and its synthesized circuit agree only up to the Trotter error; use a higher-order product formula during transpilation to tighten it. + .. note:: + ffsim does not support Windows (through its unconditional PySCF dependency), so this gate + requires the ``ffsim`` extra (``pip install "qiskit-fermions[ffsim]"``) and is unavailable + there. Use `WSL `__ on Windows. + .. caution:: This is an early development prototype. Beware of changes to its interface without warning during the pre-release development of this package. - """ - class Variant(Enum): - """The spin variant of a :class:`.UCC` ansatz. + .. invisible-code-block: python - Mirrors ffsim's UCCSD operator flavors - (:external:class:`~ffsim.UCCSDOpRestrictedReal` and - :external:class:`~ffsim.UCCSDOpUnrestrictedReal`), plus a spinless variant. Passed explicitly - to :meth:`.UCC.__init__` (or its string value) and validated against the supplied amplitude - shapes. - """ + >>> from qiskit_fermions.utils.optionals import HAS_FFSIM - SPINLESS = "spinless" - RESTRICTED = "restricted" - UNRESTRICTED = "unrestricted" + .. skip: start if(not HAS_FFSIM) - def __init__( - self, - variant: UCC.Variant | str, - t1: np.ndarray | tuple[np.ndarray, np.ndarray], - t2: np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray], - *, - antisymmetric: bool = False, - atol: float = 1e-8, - ) -> None: - r"""Initializing an instance of this gate can be done with the arguments listed below. + .. doctest:: - Args: - variant: the spin variant, a :class:`.UCC.Variant` (or its string value - ``"restricted"``, ``"unrestricted"``, or ``"spinless"``). Determines the number of - modes this gate acts on (see the class docstring) and the expected amplitude shapes. - t1: the :math:`t_1` (singles) amplitudes. For the ``"restricted"`` and ``"spinless"`` - variants, a single array of shape ``(nocc, nvrt)``. For the ``"unrestricted"`` - variant, a pair ``(t1a, t1b)``. - t2: the :math:`t_2` (doubles) amplitudes. For the ``"restricted"`` and ``"spinless"`` - variants, a single array of shape ``(nocc, nocc, nvrt, nvrt)``. For the - ``"unrestricted"`` variant, a triple ``(t2aa, t2ab, t2bb)``. - antisymmetric: whether the same-spin :math:`t_2` blocks obey the *separate* occupied and - virtual antisymmetry (the standard coupled-cluster convention, see - :attr:`.antisymmetric`). When ``True`` the supplied blocks are validated against it - and the parameter vector is restricted to the corresponding subspace. Not supported - for the ``"restricted"`` variant. - atol: the absolute tolerance for the ``antisymmetric`` validation. + >>> import ffsim + >>> import numpy as np + >>> from qiskit_fermions.circuit.library import UCC + >>> uccsd_op = ffsim.UCCSDOpRestrictedReal( + ... t1=np.zeros((1, 1)), t2=np.zeros((1, 1, 1, 1)) + ... ) + >>> gate = UCC(uccsd_op) + >>> gate.norb, gate.num_modes + (2, 4) - Raises: - ValueError: if ``variant`` is not recognized, if the amplitude shapes are inconsistent - with each other or with ``variant``, if ``antisymmetric`` is requested for the - ``"restricted"`` variant, or if ``antisymmetric`` is requested but a same-spin - :math:`t_2` block violates that antisymmetry. - """ - variant = self._normalize_variant(variant) - self._variant = variant - self._validate_antisymmetric_supported(variant, antisymmetric) - - self.antisymmetric = antisymmetric - r"""Whether the same-spin :math:`t_2` blocks obey the separate occupied/virtual antisymmetry. - - The cluster operator only ever sees the part of a same-spin :math:`t_2` block that is - symmetric under the *simultaneous* exchange - :math:`t_2[i,j,a,b] = t_2[j,i,b,a]`, because the underlying excitation - :math:`a^\dagger_a a^\dagger_b a_j a_i` is invariant under relabeling the pairs - :math:`(i,a) \leftrightarrow (j,b)`. That weaker symmetry is therefore always imposed. - - The standard coupled-cluster convention additionally makes the block antisymmetric in each - index pair *separately*, :math:`t_2[i,j,a,b] = -t_2[j,i,a,b] = -t_2[i,j,b,a]`, which is a - strict subspace of the above. Setting this flag opts into that convention: the supplied - amplitudes are validated against it, and :meth:`.num_parameters` / - :meth:`.from_parameters` / :meth:`.to_parameters` switch to the smaller parameter basis that - spans exactly this subspace. - """ - self._atol = atol - - if variant is UCC.Variant.UNRESTRICTED: - # the variant selects which argument shapes are valid; mypy cannot narrow the unions - self.t1: np.ndarray | tuple[np.ndarray, ...] = tuple( - np.asarray(t, dtype=complex) for t in cast("tuple", t1) - ) - """The :math:`t_1` (singles) amplitudes.""" - self.t2: np.ndarray | tuple[np.ndarray, ...] = tuple( - np.asarray(t, dtype=complex) for t in cast("tuple", t2) - ) - """The :math:`t_2` (doubles) amplitudes.""" - else: - self.t1 = np.asarray(cast(np.ndarray, t1), dtype=complex) - self.t2 = np.asarray(cast(np.ndarray, t2), dtype=complex) - - norb = self._validate_shapes() - self.norb = norb - """The number of spatial orbitals (or spinless modes, for the spinless variant).""" - - if antisymmetric: - # validated only after the shapes are known to be consistent, so the transpose-based - # checks below cannot fail on a malformed tensor instead of a genuinely asymmetric one - self._validate_antisymmetric_amplitudes() - - num_modes = norb if variant is UCC.Variant.SPINLESS else 2 * norb - - super().__init__("UCC", num_modes, []) - - @classmethod - def from_t_amplitudes( - cls, - t2: np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray], - *, - t1: np.ndarray | tuple[np.ndarray, np.ndarray] | None = None, - variant: UCC.Variant | str = "restricted", - antisymmetric: bool = False, - atol: float = 1e-8, - ) -> Self: - r"""Constructs a UCC ansatz from coupled-cluster :math:`t_2` (and optional :math:`t_1`) amplitudes. - - This is a convenience constructor mirroring :meth:`.UCJ.from_t_amplitudes`. Unlike the (L)UCJ - ansatz -- which *factorizes* the amplitudes into diagonal Coulomb layers -- the UCC ansatz - uses the amplitudes directly as its parameters, so this simply defaults an omitted ``t1`` to - zeros of the shape implied by ``t2``, giving a doubles-only (UCCD) ansatz. - - Args: - t2: the :math:`t_2` amplitudes. For the ``"restricted"`` and ``"spinless"`` variants, a - single array of shape ``(nocc, nocc, nvrt, nvrt)``. For the ``"unrestricted"`` - variant, a triple ``(t2aa, t2ab, t2bb)``. - t1: the optional :math:`t_1` amplitudes. For ``"unrestricted"``, a pair ``(t1a, t1b)``; - otherwise a single array of shape ``(nocc, nvrt)``. Defaults to zeros. - variant: the spin variant to build, a :class:`.UCC.Variant` (or its string value - ``"restricted"``, ``"unrestricted"``, or ``"spinless"``). - antisymmetric: whether to assert the standard coupled-cluster antisymmetry of the - same-spin :math:`t_2` blocks (see :attr:`.antisymmetric`). Amplitudes from a genuine - coupled-cluster calculation satisfy it, so this is a cheap way to confirm they - survived whatever preprocessing produced them. Not supported for the ``"restricted"`` - variant. - atol: the absolute tolerance for the ``antisymmetric`` validation. - - Returns: - The constructed :class:`.UCC` gate. - - Raises: - ValueError: if ``variant`` is not recognized, if the amplitude shapes are inconsistent - with each other or with ``variant``, if ``antisymmetric`` is requested for the - ``"restricted"`` variant, or if ``antisymmetric`` is requested but a same-spin - :math:`t_2` block violates that antisymmetry. - """ - variant = cls._normalize_variant(variant) - - if t1 is None: - t1 = cls._zero_t1_like(t2, variant) - - return cls(variant, t1, t2, antisymmetric=antisymmetric, atol=atol) - - @classmethod - def num_parameters( - cls, - norb: int, - nocc: int | tuple[int, int], - variant: UCC.Variant | str, - *, - antisymmetric: bool = False, - ) -> int: - r"""Returns the number of parameters of a UCC ansatz with the given settings. - - Args: - norb: the number of spatial orbitals (or spinless modes, for the spinless variant). - nocc: the number of occupied orbitals. For the ``"unrestricted"`` variant a pair - ``(nocc_a, nocc_b)`` giving the per-spin occupations; otherwise a single integer. - variant: the spin variant, a :class:`.UCC.Variant` (or its string value - ``"restricted"``, ``"unrestricted"``, or ``"spinless"``). - antisymmetric: whether the same-spin :math:`t_2` blocks are restricted to the standard - coupled-cluster antisymmetric subspace (see :attr:`.antisymmetric`), which needs - strictly fewer parameters. Not supported for the ``"restricted"`` variant. - - Returns: - The number of parameters. - - Raises: - ValueError: if ``variant`` is not recognized, if ``nocc`` is a pair for a variant other - than :attr:`.UCC.Variant.UNRESTRICTED` (or an integer for that variant), or if - ``antisymmetric`` is requested for the ``"restricted"`` variant. - """ - variant = cls._normalize_variant(variant) - nocc = cls._normalize_nocc(nocc, variant) - cls._validate_antisymmetric_supported(variant, antisymmetric) - - if variant is UCC.Variant.UNRESTRICTED: - nocc_a, nocc_b = cast("tuple[int, int]", nocc) - nvrt_a, nvrt_b = norb - nocc_a, norb - nocc_b - return ( - # t1a and t1b - nocc_a * nvrt_a - + nocc_b * nvrt_b - # the same-spin t2aa and t2bb blocks (the cross-spin t2ab is unconstrained) - + cls._same_spin_block_num_parameters(nocc_a, nvrt_a, antisymmetric) - + nocc_a * nocc_b * nvrt_a * nvrt_b - + cls._same_spin_block_num_parameters(nocc_b, nvrt_b, antisymmetric) - ) + .. skip: end + """ - # restricted and spinless both carry a single t1/t2 pair - nocc_int = cast(int, nocc) - nvrt = norb - nocc_int - return nocc_int * nvrt + cls._same_spin_block_num_parameters(nocc_int, nvrt, antisymmetric) - - @classmethod - def from_parameters( - cls, - params: np.ndarray, - norb: int, - nocc: int | tuple[int, int], - variant: UCC.Variant | str, - *, - antisymmetric: bool = False, - ) -> Self: - r"""Constructs a UCC ansatz from a real-valued parameter vector. - - With ``antisymmetric=False`` (the default) the parameter ordering matches ffsim's - :external:class:`~ffsim.UCCSDOpRestrictedReal` / - :external:class:`~ffsim.UCCSDOpUnrestrictedReal` convention, so a vector produced by ffsim's - own ``to_parameters`` round-trips through this method. - - With ``antisymmetric=True`` the same-spin :math:`t_2` blocks are instead built from the - smaller basis spanning the standard coupled-cluster antisymmetric subspace (see - :attr:`.antisymmetric`), so the expected vector length differs and ffsim's vectors no longer - apply. + def __init__(self, uccsd_op: Any) -> None: + """Initializing an instance of this gate can be done with the argument listed below. Args: - params: the real-valued parameter vector. - norb: the number of spatial orbitals (or spinless modes, for the spinless variant). - nocc: the number of occupied orbitals. For the ``"unrestricted"`` variant a pair - ``(nocc_a, nocc_b)``; otherwise a single integer. - variant: the spin variant, a :class:`.UCC.Variant` (or its string value - ``"restricted"``, ``"unrestricted"``, or ``"spinless"``). - antisymmetric: whether to build the same-spin :math:`t_2` blocks in the antisymmetric - subspace (see :attr:`.antisymmetric`). The resulting amplitudes then satisfy that - antisymmetry by construction. Not supported for the ``"restricted"`` variant. - - Returns: - The constructed :class:`.UCC` gate. + uccsd_op: the ffsim UCCSD operator to build the circuit from, one of + :external:class:`~ffsim.UCCSDOpRestrictedReal`, + :external:class:`~ffsim.UCCSDOpRestricted`, + :external:class:`~ffsim.UCCSDOpUnrestrictedReal` or + :external:class:`~ffsim.UCCSDOpUnrestricted`. Its type determines the spin variant + (see the class docstring). Raises: - ValueError: if ``variant`` is not recognized, if ``antisymmetric`` is requested for the - ``"restricted"`` variant, or if ``len(params)`` does not match - :meth:`.num_parameters` for the given settings. + MissingOptionalLibraryError: if ``ffsim`` is not installed. + TypeError: if ``uccsd_op`` is not one of ffsim's four UCCSD operator types. """ - variant = cls._normalize_variant(variant) - nocc = cls._normalize_nocc(nocc, variant) - cls._validate_antisymmetric_supported(variant, antisymmetric) - expected = cls.num_parameters(norb, nocc, variant, antisymmetric=antisymmetric) - if len(params) != expected: - raise ValueError( - "The number of parameters passed did not match the number expected based on the " - f"given settings. Expected {expected} but got {len(params)}." + HAS_FFSIM.require_now("UCC") + import ffsim + + restricted = (ffsim.UCCSDOpRestricted, ffsim.UCCSDOpRestrictedReal) + unrestricted = (ffsim.UCCSDOpUnrestricted, ffsim.UCCSDOpUnrestrictedReal) + if not isinstance(uccsd_op, restricted + unrestricted): + raise TypeError( + "UCC requires one of ffsim's UCCSD operators (UCCSDOpRestricted, " + "UCCSDOpRestrictedReal, UCCSDOpUnrestricted or UCCSDOpUnrestrictedReal), but got " + f"{type(uccsd_op).__name__}." ) - index = 0 - t1: np.ndarray | tuple[np.ndarray, np.ndarray] - t2: np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray] - - if variant is UCC.Variant.UNRESTRICTED: - nocc_a, nocc_b = cast("tuple[int, int]", nocc) - nvrt_a, nvrt_b = norb - nocc_a, norb - nocc_b - # ffsim's ordering is t1a, t1b, t2aa, t2ab, t2bb -- note that the cross-spin block sits - # *between* the two same-spin blocks, not after them - t1a, index = cls._t1_from_parameters(params, index, nocc_a, nvrt_a) - t1b, index = cls._t1_from_parameters(params, index, nocc_b, nvrt_b) - t2aa, index = cls._same_spin_t2_from_parameters( - params, index, nocc_a, nvrt_a, antisymmetric - ) - shape_ab = (nocc_a, nocc_b, nvrt_a, nvrt_b) - n_ab = int(np.prod(shape_ab)) - t2ab = np.asarray(params[index : index + n_ab], dtype=float).reshape(shape_ab) - index += n_ab - t2bb, index = cls._same_spin_t2_from_parameters( - params, index, nocc_b, nvrt_b, antisymmetric - ) - return cls(variant, (t1a, t1b), (t2aa, t2ab, t2bb), antisymmetric=antisymmetric) - - nocc_int = cast(int, nocc) - nvrt = norb - nocc_int - t1, index = cls._t1_from_parameters(params, index, nocc_int, nvrt) - t2, index = cls._same_spin_t2_from_parameters(params, index, nocc_int, nvrt, antisymmetric) - - return cls(variant, t1, t2, antisymmetric=antisymmetric) + self.uccsd_op = uccsd_op + """The ffsim UCCSD operator this gate builds its circuit from.""" - def to_parameters(self) -> np.ndarray: - r"""Converts this UCC ansatz to a real-valued parameter vector. + self._unrestricted = isinstance(uccsd_op, unrestricted) - The inverse of :meth:`.from_parameters`, using the same ordering and the same basis this - gate's :attr:`.antisymmetric` flag selects -- so ``from_parameters(gate.to_parameters(), ...)`` - round-trips as long as the flag is passed consistently. + super().__init__("UCC", 2 * uccsd_op.norb, []) - .. note:: - Only the independent amplitude entries implied by the variant's symmetries (and by - :attr:`.antisymmetric`) are written out; see :meth:`.num_parameters`. Amplitudes violating - those symmetries -- or carrying a non-negligible imaginary part -- are therefore not - recoverable from the parameter vector. + @property + def norb(self) -> int: + """The number of spatial orbitals.""" + return int(self.uccsd_op.norb) - .. note:: - The round-trip is two-sided and holds at *any* parameter scale: the amplitudes are this - ansatz's parameters directly, so both directions are a plain re-indexing. - - Returns: - The real-valued parameter vector. - """ - cls = type(self) - variant = self._variant - antisymmetric = self.antisymmetric - params = np.empty( - cls.num_parameters(self.norb, self._nocc(), variant, antisymmetric=antisymmetric) - ) - - index = 0 - if variant is UCC.Variant.UNRESTRICTED: - t1a, t1b = cast("tuple[np.ndarray, np.ndarray]", self.t1) - t2aa, t2ab, t2bb = cast("tuple[np.ndarray, np.ndarray, np.ndarray]", self.t2) - index = cls._t1_to_parameters(params, index, t1a) - index = cls._t1_to_parameters(params, index, t1b) - index = cls._same_spin_t2_to_parameters(params, index, t2aa, antisymmetric) - params[index : index + t2ab.size] = t2ab.real.ravel() - index += t2ab.size - cls._same_spin_t2_to_parameters(params, index, t2bb, antisymmetric) - return params - - index = cls._t1_to_parameters(params, index, cast(np.ndarray, self.t1)) - cls._same_spin_t2_to_parameters(params, index, cast(np.ndarray, self.t2), antisymmetric) - return params + def _nocc(self) -> int | tuple[int, int]: + """Returns the number of occupied orbitals, per spin sector when unrestricted.""" + if self._unrestricted: + t1a, t1b = cast("tuple[np.ndarray, np.ndarray]", self.uccsd_op.t1) + return int(t1a.shape[0]), int(t1b.shape[0]) + return int(cast(np.ndarray, self.uccsd_op.t1).shape[0]) def cluster_operator(self) -> FermionOperator: r"""Returns the anti-Hermitian cluster generator :math:`T - T^\dagger`. The generator is expressed in the block-spin mode convention (mode ``p`` is alpha orbital - ``p``, mode ``norb + p`` is beta orbital ``p``) for the spinful variants, and directly on the - ``norb`` modes for the spinless variant. Occupied orbitals are ordered before virtual ones. + ``p``, mode ``norb + p`` is beta orbital ``p``). Occupied orbitals are ordered before virtual + ones. Being anti-Hermitian, this generator relates to the ansatz unitary by :math:`e^{T - T^\dagger} = e^{-i H}` with the Hermitian :math:`H = i (T - T^\dagger)`. That @@ -513,319 +229,6 @@ def support(actions: Sequence[tuple[bool, int]]) -> tuple[int, ...]: return FermionOperator.from_terms_with_groups(terms_with_groups) - @staticmethod - def _normalize_variant(variant: UCC.Variant | str) -> UCC.Variant: - """Normalizes a variant argument, raising the same ``ValueError`` as ``__init__``.""" - try: - return UCC.Variant(variant) - except ValueError: - raise ValueError( - f"Unknown UCC variant {variant!r}; expected 'restricted', 'unrestricted', or " - "'spinless'." - ) from None - - @staticmethod - def _normalize_nocc(nocc: int | tuple[int, int], variant: UCC.Variant) -> int | tuple[int, int]: - """Validates ``nocc`` against ``variant``, rejecting the wrong arity. - - The unrestricted variant resolves its occupied/virtual split per spin sector, so it requires - a ``(nocc_a, nocc_b)`` pair; the other variants take a single integer. A mismatch would - otherwise be silently misinterpreted (a pair truncated to its first element, or an integer - broadcast to both sectors), so reject it explicitly. - """ - is_pair = isinstance(nocc, tuple) - if variant is UCC.Variant.UNRESTRICTED and not is_pair: - raise ValueError( - f"The 'unrestricted' variant requires a (nocc_a, nocc_b) pair; got nocc={nocc!r}." - ) - if variant is not UCC.Variant.UNRESTRICTED and is_pair: - raise ValueError( - f"A tuple nocc={nocc!r} is only valid for the 'unrestricted' variant; pass a single " - f"integer for the '{variant.value}' variant." - ) - return nocc - - @staticmethod - def _validate_antisymmetric_supported(variant: UCC.Variant, antisymmetric: bool) -> None: - r"""Rejects ``antisymmetric=True`` for the ``restricted`` variant. - - The restricted ``t2`` does double duty: the *same* tensor supplies both the same-spin - amplitudes (where the separate occupied/virtual antisymmetry is the standard convention) and - the alpha-beta amplitudes (where it is not -- there is no antisymmetry to impose between an - alpha and a beta excitation, since exchanging them is not an exchange of identical operators). - Imposing the antisymmetry on the restricted tensor would therefore silently constrain the - cross-spin channel too, yielding an over-constrained ansatz. Refuse explicitly rather than - offer a flag that means something different here than everywhere else. - """ - if antisymmetric and variant is UCC.Variant.RESTRICTED: - raise ValueError( - "antisymmetric=True is not supported for the 'restricted' variant: its single t2 " - "tensor supplies both the same-spin amplitudes (where the antisymmetry is the " - "standard convention) and the alpha-beta amplitudes (where it does not apply), so " - "imposing it would also constrain the cross-spin channel. Use the 'unrestricted' " - "variant to constrain only the same-spin blocks, or pass antisymmetric=False." - ) - - def _validate_antisymmetric_amplitudes(self) -> None: - r"""Validates that the same-spin ``t2`` blocks obey the separate occupied/virtual antisymmetry. - - Only the same-spin blocks are checked; the unrestricted cross-spin ``t2ab`` carries no such - symmetry. The ``restricted`` variant never reaches here (see - :meth:`._validate_antisymmetric_supported`). - """ - if self._variant is UCC.Variant.UNRESTRICTED: - t2aa, _, t2bb = cast("tuple[np.ndarray, np.ndarray, np.ndarray]", self.t2) - blocks = [("t2aa", t2aa), ("t2bb", t2bb)] - else: # SPINLESS: the single register is exactly one same-spin sector - blocks = [("t2", cast(np.ndarray, self.t2))] - - for name, t2 in blocks: - self._check_antisymmetric_block(t2, name, self._atol) - - @staticmethod - def _check_antisymmetric_block(t2: np.ndarray, name: str, atol: float) -> None: - """Raises if ``t2`` is not antisymmetric under occupied and virtual exchange separately. - - ``t2 + t2.transpose(...)`` is exactly the symmetric (i.e. antisymmetry-violating) part of the - tensor under each exchange, so the check needs no index loops. - """ - for label, violation in ( - ("occupied (i <-> j)", t2 + t2.transpose(1, 0, 2, 3)), - ("virtual (a <-> b)", t2 + t2.transpose(0, 1, 3, 2)), - ): - worst = float(np.abs(violation).max(initial=0.0)) - if worst > atol: - raise ValueError( - f"antisymmetric=True was requested but the {name} amplitudes are not " - f"antisymmetric under {label} exchange: largest violation {worst:.3e} exceeds " - f"atol={atol:.3e}. Pass antisymmetric=False to allow the full " - "exchange-symmetric family." - ) - - @classmethod - def _same_spin_block_num_parameters(cls, nocc: int, nvrt: int, antisymmetric: bool) -> int: - """Returns the parameter count of one same-spin ``t2`` block, per symmetry convention.""" - if antisymmetric: - # the independent entries are the strictly-upper-triangular occupied pairs combined with - # the strictly-upper-triangular virtual pairs; everything else follows by antisymmetry - # (and the i == j / a == b entries must vanish) - return len(cls._antisymmetric_double_indices(nocc, nvrt)) - # the upper triangle of the combined occupied-virtual pair index - n_pairs = nocc * nvrt - return n_pairs * (n_pairs + 1) // 2 - - @staticmethod - def _antisymmetric_double_indices(nocc: int, nvrt: int) -> list[tuple[int, int, int, int]]: - """Returns the independent ``(i, j, a, b)`` entries of an antisymmetric same-spin ``t2``.""" - return [ - (i, j, a, b) - for i, j in itertools.combinations(range(nocc), 2) - for a, b in itertools.combinations(range(nvrt), 2) - ] - - @staticmethod - def _zero_t1_like( - t2: np.ndarray | tuple[np.ndarray, np.ndarray, np.ndarray], variant: UCC.Variant - ) -> np.ndarray | tuple[np.ndarray, np.ndarray]: - """Builds zero ``t1`` amplitudes of the shape implied by ``t2`` (a doubles-only ansatz).""" - if variant is UCC.Variant.UNRESTRICTED: - t2ab = np.asarray(cast("tuple", t2)[1]) - if t2ab.ndim != 4: - raise ValueError( - "The unrestricted t2ab amplitudes must be 4-dimensional " - f"(nocc_a, nocc_b, nvrt_a, nvrt_b); got shape {t2ab.shape}." - ) - nocc_a, nocc_b, nvrt_a, nvrt_b = t2ab.shape - return np.zeros((nocc_a, nvrt_a)), np.zeros((nocc_b, nvrt_b)) - t2_arr = np.asarray(cast(np.ndarray, t2)) - if t2_arr.ndim != 4: - raise ValueError( - "The t2 amplitudes must be 4-dimensional (nocc, nocc, nvrt, nvrt); got shape " - f"{t2_arr.shape}." - ) - nocc, _, nvrt, _ = t2_arr.shape - return np.zeros((nocc, nvrt)) - - @staticmethod - def _t1_from_parameters( - params: np.ndarray, index: int, nocc: int, nvrt: int - ) -> tuple[np.ndarray, int]: - """Reads an ``(nocc, nvrt)`` ``t1`` block off ``params``, returning it and the new index.""" - n = nocc * nvrt - t1 = np.asarray(params[index : index + n], dtype=float).reshape((nocc, nvrt)) - return t1, index + n - - @staticmethod - def _t1_to_parameters(params: np.ndarray, index: int, t1: np.ndarray) -> int: - """Writes ``t1`` into ``params`` at ``index``, returning the new index.""" - n = int(t1.size) - params[index : index + n] = t1.real.ravel() - return index + n - - @staticmethod - def _occ_vrt_pairs(nocc: int, nvrt: int) -> list[tuple[int, int]]: - """Returns the ``(i, a)`` occupied-virtual index pairs, in ffsim's parameter order.""" - return list(itertools.product(range(nocc), range(nvrt))) - - @classmethod - def _same_spin_t2_from_parameters( - cls, params: np.ndarray, index: int, nocc: int, nvrt: int, antisymmetric: bool - ) -> tuple[np.ndarray, int]: - """Reads one same-spin ``t2`` block off ``params``, in the basis ``antisymmetric`` selects.""" - if antisymmetric: - return cls._antisymmetric_t2_from_parameters(params, index, nocc, nvrt) - return cls._exchange_symmetric_t2_from_parameters(params, index, nocc, nvrt) - - @classmethod - def _same_spin_t2_to_parameters( - cls, params: np.ndarray, index: int, t2: np.ndarray, antisymmetric: bool - ) -> int: - """Writes one same-spin ``t2`` block into ``params``, inverting :meth:`._same_spin_t2_from_parameters`.""" - if antisymmetric: - return cls._antisymmetric_t2_to_parameters(params, index, t2) - return cls._exchange_symmetric_t2_to_parameters(params, index, t2) - - @classmethod - def _exchange_symmetric_t2_from_parameters( - cls, params: np.ndarray, index: int, nocc: int, nvrt: int - ) -> tuple[np.ndarray, int]: - """Reads a ``t2`` block off ``params``, filling in its ``(i, a) <-> (j, b)`` symmetry. - - Every ``t2`` block this ansatz parameterizes -- the restricted and spinless single block, and - each same-spin block of the unrestricted variant -- is symmetric under the *simultaneous* - exchange of the occupied and virtual indices, i.e. ``t2[i, j, a, b] == t2[j, i, b, a]``. It is - therefore determined by the upper triangle of the occupied-virtual pair index. This mirrors - ffsim's ``UCCSDOpRestrictedReal``/``UCCSDOpUnrestrictedReal`` ``from_parameters`` ordering - exactly: the parameters run over ``combinations_with_replacement`` of those pairs, and each - parameter sets both ``t2[i, j, a, b]`` and its ``(j, i, b, a)`` partner. - """ - t2 = np.zeros((nocc, nocc, nvrt, nvrt)) - for (i, a), (j, b) in itertools.combinations_with_replacement( - cls._occ_vrt_pairs(nocc, nvrt), 2 - ): - t2[i, j, a, b] = params[index] - t2[j, i, b, a] = params[index] - index += 1 - return t2, index - - @classmethod - def _exchange_symmetric_t2_to_parameters( - cls, params: np.ndarray, index: int, t2: np.ndarray - ) -> int: - """Inverts :meth:`._exchange_symmetric_t2_from_parameters`, reading the same entries back.""" - nocc, _, nvrt, _ = t2.shape - for (i, a), (j, b) in itertools.combinations_with_replacement( - cls._occ_vrt_pairs(nocc, nvrt), 2 - ): - params[index] = t2[i, j, a, b].real - index += 1 - return index - - @classmethod - def _antisymmetric_t2_from_parameters( - cls, params: np.ndarray, index: int, nocc: int, nvrt: int - ) -> tuple[np.ndarray, int]: - """Expands the independent ``(i < j, a < b)`` entries into a fully antisymmetric ``t2``. - - Each parameter fixes four entries of the block via - ``t2[i, j, a, b] == -t2[j, i, a, b] == -t2[i, j, b, a] == t2[j, i, b, a]``, so the result - satisfies the separate occupied and virtual antisymmetry by construction (see - :attr:`.antisymmetric`). Entries with ``i == j`` or ``a == b`` are left zero, as that - antisymmetry requires. - """ - t2 = np.zeros((nocc, nocc, nvrt, nvrt)) - for i, j, a, b in cls._antisymmetric_double_indices(nocc, nvrt): - value = params[index] - t2[i, j, a, b] = value - t2[j, i, b, a] = value - t2[j, i, a, b] = -value - t2[i, j, b, a] = -value - index += 1 - return t2, index - - @classmethod - def _antisymmetric_t2_to_parameters(cls, params: np.ndarray, index: int, t2: np.ndarray) -> int: - """Inverts :meth:`._antisymmetric_t2_from_parameters`, reading the independent entries back.""" - nocc, _, nvrt, _ = t2.shape - for i, j, a, b in cls._antisymmetric_double_indices(nocc, nvrt): - params[index] = t2[i, j, a, b].real - index += 1 - return index - - def _nocc(self) -> int | tuple[int, int]: - """Returns the number of occupied orbitals, per spin sector for the unrestricted variant.""" - if self._variant is UCC.Variant.UNRESTRICTED: - t1a, t1b = cast("tuple[np.ndarray, np.ndarray]", self.t1) - return int(t1a.shape[0]), int(t1b.shape[0]) - return int(cast(np.ndarray, self.t1).shape[0]) - - def _validate_shapes(self) -> int: - """Validates the amplitude shapes, returning the implied ``norb``. - - Raises: - ValueError: if the amplitude shapes are inconsistent with each other, or if the - unrestricted variant's two spin sectors imply different orbital counts. - """ - variant = self._variant - - if variant is not UCC.Variant.UNRESTRICTED: - t1_arr = cast(np.ndarray, self.t1) - t2_arr = cast(np.ndarray, self.t2) - if t1_arr.ndim != 2: - raise ValueError( - "The t1 amplitudes must be 2-dimensional (nocc, nvrt); got shape " - f"{t1_arr.shape}." - ) - nocc, nvrt = t1_arr.shape - if t2_arr.shape != (nocc, nocc, nvrt, nvrt): - raise ValueError( - f"Inconsistent {variant.value} UCC amplitude shapes: t2 should have shape " - f"{(nocc, nocc, nvrt, nvrt)} (implied by t1's shape {t1_arr.shape}) but got " - f"{t2_arr.shape}." - ) - return int(nocc + nvrt) - - t1_tuple = cast("tuple[np.ndarray, ...]", self.t1) - t2_tuple = cast("tuple[np.ndarray, ...]", self.t2) - if len(t1_tuple) != 2: - raise ValueError( - "The 'unrestricted' variant expects a (t1a, t1b) pair of t1 amplitudes; got " - f"{len(t1_tuple)} array(s)." - ) - if len(t2_tuple) != 3: - raise ValueError( - "The 'unrestricted' variant expects a (t2aa, t2ab, t2bb) triple of t2 amplitudes; " - f"got {len(t2_tuple)} array(s)." - ) - t1a, t1b = t1_tuple - t2aa, t2ab, t2bb = t2_tuple - if t1a.ndim != 2 or t1b.ndim != 2: - raise ValueError( - "The unrestricted t1 amplitudes must each be 2-dimensional (nocc, nvrt); got " - f"shapes {t1a.shape} and {t1b.shape}." - ) - nocc_a, nvrt_a = t1a.shape - nocc_b, nvrt_b = t1b.shape - # both spin sectors span the same set of spatial orbitals, so their occupied/virtual splits - # must add up to the same norb -- otherwise the block-spin register is ill-defined - if nocc_a + nvrt_a != nocc_b + nvrt_b: - raise ValueError( - "The unrestricted alpha and beta t1 amplitudes imply different numbers of spatial " - f"orbitals: {nocc_a + nvrt_a} (alpha) vs {nocc_b + nvrt_b} (beta)." - ) - expected = { - "t2aa": ((nocc_a, nocc_a, nvrt_a, nvrt_a), t2aa.shape), - "t2ab": ((nocc_a, nocc_b, nvrt_a, nvrt_b), t2ab.shape), - "t2bb": ((nocc_b, nocc_b, nvrt_b, nvrt_b), t2bb.shape), - } - for name, (want, got) in expected.items(): - if got != want: - raise ValueError( - f"Inconsistent unrestricted UCC amplitude shapes: {name} should have shape " - f"{want} (implied by the t1 amplitudes) but got {got}." - ) - return int(nocc_a + nvrt_a) - def _excitation_operator(self) -> FermionOperator: r"""Builds the (non-Hermitian) cluster operator :math:`T = T_1 + T_2`. @@ -834,29 +237,22 @@ def _excitation_operator(self) -> FermionOperator: restricted same-spin doubles, where ``(i, j, a, b)`` and ``(j, i, b, a)`` produce the same actions); their coefficients must add rather than overwrite each other. """ - variant = self._variant norb = self.norb - # spinless acts on a single register, so its (single) sector carries no spin offset - spinless = variant is UCC.Variant.SPINLESS - terms: dict[tuple, complex] = {} def mode(sigma: int, orb: int) -> int: """Maps a ``(spin sector, orbital)`` pair onto its block-spin mode index.""" - return orb if spinless else orb + sigma * norb + return orb + sigma * norb def add(actions: tuple, coeff: complex) -> None: if coeff == 0.0: return terms[actions] = terms.get(actions, 0.0) + coeff - match variant: - case UCC.Variant.UNRESTRICTED: - self._add_unrestricted_terms(add, mode) - case UCC.Variant.RESTRICTED: - self._add_restricted_terms(add, mode) - case _: # SPINLESS - self._add_spinless_terms(add, mode) + if self._unrestricted: + self._add_unrestricted_terms(add, mode) + else: + self._add_restricted_terms(add, mode) return FermionOperator.from_dict(terms) # type: ignore[arg-type] @@ -867,8 +263,8 @@ def _add_restricted_terms(self, add, mode) -> None: convention: the singles are applied identically in both spin sectors, and the doubles carry a factor ``1/2`` on the two same-spin blocks and ``1`` on the alpha-beta block. """ - t1 = cast(np.ndarray, self.t1) - t2 = cast(np.ndarray, self.t2) + t1 = cast(np.ndarray, self.uccsd_op.t1) + t2 = cast(np.ndarray, self.uccsd_op.t2) nocc, nvrt = t1.shape for i, a in itertools.product(range(nocc), range(nvrt)): @@ -906,8 +302,8 @@ def _add_unrestricted_terms(self, add, mode) -> None: block. Note that each spin sector uses its *own* occupied count as the virtual-orbital offset, so the two sectors may split the same ``norb`` orbitals differently. """ - t1a, t1b = cast("tuple[np.ndarray, np.ndarray]", self.t1) - t2aa, t2ab, t2bb = cast("tuple[np.ndarray, np.ndarray, np.ndarray]", self.t2) + t1a, t1b = cast("tuple[np.ndarray, np.ndarray]", self.uccsd_op.t1) + t2aa, t2ab, t2bb = cast("tuple[np.ndarray, np.ndarray, np.ndarray]", self.uccsd_op.t2) nocc = (t1a.shape[0], t1b.shape[0]) nvrt = (t1a.shape[1], t1b.shape[1]) @@ -941,30 +337,6 @@ def _add_unrestricted_terms(self, add, mode) -> None: complex(t2ab[i, j, a, b]), ) - def _add_spinless_terms(self, add, mode) -> None: - r"""Adds the spinless :math:`T_1` and :math:`T_2` terms on the single ``norb``-mode register. - - The doubles carry the same ``1/4`` factor as an unrestricted same-spin block, since a - spinless register is exactly one such sector. - """ - t1 = cast(np.ndarray, self.t1) - t2 = cast(np.ndarray, self.t2) - nocc, nvrt = t1.shape - - for i, a in itertools.product(range(nocc), range(nvrt)): - add((cre(mode(0, nocc + a)), ann(mode(0, i))), complex(t1[i, a])) - - for i, j, a, b in itertools.product(range(nocc), range(nocc), range(nvrt), range(nvrt)): - add( - ( - cre(mode(0, nocc + a)), - cre(mode(0, nocc + b)), - ann(mode(0, j)), - ann(mode(0, i)), - ), - 0.25 * complex(t2[i, j, a, b]), - ) - def _apply_unitary_placed_( self, vec: np.ndarray, @@ -1017,6 +389,15 @@ def _build_definition(self) -> FermionicCircuit: definition.modes, ) + # ffsim's UCCSD operators can carry a trailing orbital rotation; it acts per spin sector + if self.uccsd_op.final_orbital_rotation is not None: + norb = self.norb + rotation = OrbitalRotation(self.uccsd_op.final_orbital_rotation) + definition.append(rotation, definition.modes[:norb]) + definition.append( + OrbitalRotation(self.uccsd_op.final_orbital_rotation), definition.modes[norb:] + ) + return definition def _define(self) -> None: diff --git a/releasenotes/notes/ucc-from-ffsim-op-88e7c005f5390543.yaml b/releasenotes/notes/ucc-from-ffsim-op-88e7c005f5390543.yaml new file mode 100644 index 000000000..302b754dc --- /dev/null +++ b/releasenotes/notes/ucc-from-ffsim-op-88e7c005f5390543.yaml @@ -0,0 +1,35 @@ +--- +upgrade: + - | + The :class:`.UCC` gate is now constructed from one of ffsim's UCCSD operators, which is its only + constructor argument. The spin variant and the number of modes are read off that operator, so the + ``variant`` argument and the ``UCC.Variant`` enum are gone: + + .. code-block:: python + + # before + ansatz = UCC("restricted", t1, t2) + + # after + ansatz = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) + + Accordingly, ``UCC.from_t_amplitudes``, ``UCC.num_parameters``, ``UCC.from_parameters`` and + ``UCC.to_parameters`` have been removed, as have the ``spinless`` variant and the opt-in + ``antisymmetric`` parameterization, which have no ffsim equivalent. ffsim's operators provide + :external:meth:`~ffsim.UCCSDOpRestrictedReal.n_params`, + :external:meth:`~ffsim.UCCSDOpRestrictedReal.from_parameters` and + :external:meth:`~ffsim.UCCSDOpRestrictedReal.to_parameters` with identical conventions. To build + an ansatz outside that family (a spinless one, or an antisymmetrized :math:`t_2`), construct an + :class:`.Evolution` over your own cluster operator directly, which also gives you control over + the Trotter ordering of its terms. + + The wrapped operator is available as ``UCC.uccsd_op``, so its amplitudes are reachable as + ``gate.uccsd_op.t1`` and ``gate.uccsd_op.t2``; the gate no longer mirrors them as attributes of + its own. :meth:`.UCC.cluster_operator` is unchanged. A + ``final_orbital_rotation`` carried by the ffsim operator is now appended as a closing + :class:`.OrbitalRotation`. + + Since ffsim is now the gate's input type rather than an optional accelerator, constructing a + :class:`.UCC` requires the ``ffsim`` extra (``pip install "qiskit-fermions[ffsim]"``) and raises + :class:`~qiskit.exceptions.MissingOptionalLibraryError` without it. ffsim does not support + Windows, so the gate is unavailable there; use WSL. diff --git a/tests/python/circuit/library/test_ucc.py b/tests/python/circuit/library/test_ucc.py index 9a67822e0..0cbc5e52b 100644 --- a/tests/python/circuit/library/test_ucc.py +++ b/tests/python/circuit/library/test_ucc.py @@ -10,7 +10,13 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -"""Structural tests for the UCC ansatz gate.""" +"""Structural tests for the UCC ansatz gate. + +The amplitudes and their conventions belong to ffsim, so the tests here cover only what this package +adds: reading the spin variant off the ffsim operator type, the cluster generator built from the +amplitudes, and the conjugate-paired grouping that keeps every factor of the product formula +unitary. +""" from __future__ import annotations @@ -19,216 +25,72 @@ from qiskit_fermions.circuit import FermionicCircuit from qiskit_fermions.circuit.library import UCC - -def _restricted_amplitudes(nocc, nvrt, *, seed): - """Returns random ``(t1, t2)`` amplitudes with the restricted ``(i, a) <-> (j, b)`` symmetry.""" - rng = np.random.default_rng(seed) - t1 = rng.standard_normal((nocc, nvrt)) - t2 = rng.standard_normal((nocc, nocc, nvrt, nvrt)) - return t1, t2 + t2.transpose(1, 0, 3, 2) +ffsim = pytest.importorskip("ffsim") -def _unrestricted_amplitudes(nocc_a, nocc_b, nvrt_a, nvrt_b, *, seed): - """Returns random unrestricted ``((t1a, t1b), (t2aa, t2ab, t2bb))`` amplitudes.""" +def _restricted_op(nocc, nvrt, *, seed, with_final=False): + """Returns an ffsim restricted UCCSD operator with random real amplitudes.""" + rng = np.random.default_rng(seed) + t1 = rng.standard_normal((nocc, nvrt)) * 0.05 + t2 = rng.standard_normal((nocc, nocc, nvrt, nvrt)) * 0.05 + t2 = 0.5 * (t2 + t2.transpose(1, 0, 3, 2)) # the symmetry the cluster operator sees + final = None + if with_final: + norb = nocc + nvrt + final = np.eye(norb) + return ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2, final_orbital_rotation=final) + + +def _unrestricted_op(nocc_a, nocc_b, nvrt_a, nvrt_b, *, seed): + """Returns an ffsim unrestricted UCCSD operator with random real amplitudes.""" rng = np.random.default_rng(seed) - t1a = rng.standard_normal((nocc_a, nvrt_a)) - t1b = rng.standard_normal((nocc_b, nvrt_b)) - t2aa = rng.standard_normal((nocc_a, nocc_a, nvrt_a, nvrt_a)) - t2bb = rng.standard_normal((nocc_b, nocc_b, nvrt_b, nvrt_b)) - t2ab = rng.standard_normal((nocc_a, nocc_b, nvrt_a, nvrt_b)) - t2aa = t2aa + t2aa.transpose(1, 0, 3, 2) - t2bb = t2bb + t2bb.transpose(1, 0, 3, 2) - return (t1a, t1b), (t2aa, t2ab, t2bb) - - -def test_ucc_restricted_variant(): - """Restricted amplitudes act on ``2 * norb`` block-spin modes.""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=0) - gate = UCC("restricted", t1, t2) - assert gate.norb == nocc + nvrt - assert gate.num_modes == 2 * (nocc + nvrt) - assert gate._variant is UCC.Variant.RESTRICTED - - -def test_ucc_unrestricted_variant(): - """Unrestricted amplitudes act on ``2 * norb`` block-spin modes.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=1) - gate = UCC("unrestricted", t1, t2) + t1 = ( + rng.standard_normal((nocc_a, nvrt_a)) * 0.05, + rng.standard_normal((nocc_b, nvrt_b)) * 0.05, + ) + t2aa = rng.standard_normal((nocc_a, nocc_a, nvrt_a, nvrt_a)) * 0.05 + t2bb = rng.standard_normal((nocc_b, nocc_b, nvrt_b, nvrt_b)) * 0.05 + t2ab = rng.standard_normal((nocc_a, nocc_b, nvrt_a, nvrt_b)) * 0.05 + t2aa = 0.5 * (t2aa + t2aa.transpose(1, 0, 3, 2)) + t2bb = 0.5 * (t2bb + t2bb.transpose(1, 0, 3, 2)) + return ffsim.UCCSDOpUnrestrictedReal(t1=t1, t2=(t2aa, t2ab, t2bb)) + + +def test_ucc_restricted_reads_norb_off_the_operator(): + """A restricted operator gives a ``2 * norb``-mode gate.""" + op = _restricted_op(2, 2, seed=13) + gate = UCC(op) + assert gate.uccsd_op is op assert gate.norb == 4 assert gate.num_modes == 8 - assert gate._variant is UCC.Variant.UNRESTRICTED - - -def test_ucc_spinless_variant(): - """Spinless amplitudes act on a single ``norb``-mode register.""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=2) - gate = UCC("spinless", t1, t2) - assert gate.norb == nocc + nvrt - assert gate.num_modes == nocc + nvrt - assert gate._variant is UCC.Variant.SPINLESS - - -def test_ucc_variant_accepts_enum_and_string(): - """The variant may be passed either as the enum or as its string value.""" - t1, t2 = _restricted_amplitudes(1, 2, seed=3) - assert UCC(UCC.Variant.RESTRICTED, t1, t2)._variant is UCC.Variant.RESTRICTED - assert UCC("restricted", t1, t2)._variant is UCC.Variant.RESTRICTED - - -def test_ucc_unknown_variant_raises(): - """An unrecognized variant is rejected with a helpful message.""" - t1, t2 = _restricted_amplitudes(1, 1, seed=4) - with pytest.raises(ValueError, match="Unknown UCC variant"): - UCC("balanced", t1, t2) - - -def test_ucc_inconsistent_t2_shape_raises(): - """A ``t2`` whose shape is not implied by ``t1`` is rejected.""" - t1 = np.zeros((2, 2)) - with pytest.raises(ValueError, match="t2 should have shape"): - UCC("restricted", t1, np.zeros((2, 2, 3, 3))) - - -def test_ucc_non_2d_t1_raises(): - """A ``t1`` that is not a 2-dimensional ``(nocc, nvrt)`` matrix is rejected.""" - with pytest.raises(ValueError, match="must be 2-dimensional"): - UCC("restricted", np.zeros(4), np.zeros((2, 2, 2, 2))) - - -def test_ucc_unrestricted_mismatched_norb_raises(): - """Unrestricted alpha/beta amplitudes implying different orbital counts are rejected.""" - t1a = np.zeros((2, 2)) # norb = 4 - t1b = np.zeros((1, 3)) # norb = 4 - t1b_bad = np.zeros((1, 4)) # norb = 5 - t2 = (np.zeros((2, 2, 2, 2)), np.zeros((2, 1, 2, 4)), np.zeros((1, 1, 4, 4))) - with pytest.raises(ValueError, match="different numbers of spatial orbitals"): - UCC("unrestricted", (t1a, t1b_bad), t2) - # the consistent counterpart is accepted (guards against the test passing for the wrong reason) - good_t2 = (np.zeros((2, 2, 2, 2)), np.zeros((2, 1, 2, 3)), np.zeros((1, 1, 3, 3))) - assert UCC("unrestricted", (t1a, t1b), good_t2).norb == 4 - - -def test_ucc_unrestricted_wrong_arity_raises(): - """The unrestricted variant requires a ``t1`` pair and a ``t2`` triple.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=5) - with pytest.raises(ValueError, match=r"\(t2aa, t2ab, t2bb\) triple"): - UCC("unrestricted", t1, t2[:2]) - with pytest.raises(ValueError, match=r"\(t1a, t1b\) pair"): - UCC("unrestricted", (*t1, t1[0]), t2) - - -def test_ucc_unrestricted_wrong_block_shape_raises(): - """An unrestricted ``t2`` block inconsistent with the ``t1`` amplitudes is rejected.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=6) - bad = (t2[0], np.zeros((2, 1, 2, 2)), t2[2]) - with pytest.raises(ValueError, match="t2ab should have shape"): - UCC("unrestricted", t1, bad) - - -def test_ucc_from_t_amplitudes_defaults_t1_to_zero(): - """Omitting ``t1`` yields a doubles-only (UCCD) ansatz with zero singles amplitudes.""" - _, t2 = _restricted_amplitudes(2, 2, seed=7) - gate = UCC.from_t_amplitudes(t2) - assert gate._variant is UCC.Variant.RESTRICTED - np.testing.assert_allclose(gate.t1, 0.0) - np.testing.assert_allclose(gate.t2, t2) - - -def test_ucc_from_t_amplitudes_unrestricted_defaults_t1_to_zero(): - """Omitting ``t1`` for the unrestricted variant yields per-spin zero singles amplitudes.""" - _, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=8) - gate = UCC.from_t_amplitudes(t2, variant="unrestricted") - t1a, t1b = gate.t1 - assert t1a.shape == (2, 2) - assert t1b.shape == (1, 3) - np.testing.assert_allclose(t1a, 0.0) - np.testing.assert_allclose(t1b, 0.0) - - -@pytest.mark.parametrize( - ("variant", "nocc", "nvrt"), - [("restricted", 2, 2), ("spinless", 2, 2), ("restricted", 1, 3)], -) -def test_ucc_num_parameters_matches_to_parameters(variant, nocc, nvrt): - """``num_parameters`` agrees with the length of the vector ``to_parameters`` produces.""" - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=9) - gate = UCC(variant, t1, t2) - expected = UCC.num_parameters(nocc + nvrt, nocc, variant) - assert len(gate.to_parameters()) == expected - - -def test_ucc_num_parameters_unrestricted_matches_to_parameters(): - """``num_parameters`` agrees with ``to_parameters`` for the unrestricted variant.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=10) - gate = UCC("unrestricted", t1, t2) - assert len(gate.to_parameters()) == UCC.num_parameters(4, (2, 1), "unrestricted") - - -@pytest.mark.parametrize("variant", ["restricted", "spinless"]) -def test_ucc_parameters_round_trip(variant): - """``from_parameters`` inverts ``to_parameters`` for the single-tensor variants.""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=11) - gate = UCC(variant, t1, t2) - params = gate.to_parameters() - rebuilt = UCC.from_parameters(params, nocc + nvrt, nocc, variant) - np.testing.assert_allclose(rebuilt.to_parameters(), params, atol=1e-12) - np.testing.assert_allclose(rebuilt.t1, t1, atol=1e-12) - np.testing.assert_allclose(rebuilt.t2, t2, atol=1e-12) - - -def test_ucc_parameters_round_trip_unrestricted(): - """``from_parameters`` inverts ``to_parameters`` for the unrestricted variant.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=12) - gate = UCC("unrestricted", t1, t2) - params = gate.to_parameters() - rebuilt = UCC.from_parameters(params, 4, (2, 1), "unrestricted") - np.testing.assert_allclose(rebuilt.to_parameters(), params, atol=1e-12) - for got, want in zip(rebuilt.t2, t2, strict=True): - np.testing.assert_allclose(got, want, atol=1e-12) - - -@pytest.mark.parametrize( - ("variant", "nocc", "antisymmetric"), - [ - ("restricted", 2, False), - ("unrestricted", (3, 2), False), - ("unrestricted", (3, 2), True), - ("spinless", 3, False), - ("spinless", 3, True), - ], -) -def test_ucc_parameters_round_trip_is_two_sided_at_any_scale(variant, nocc, antisymmetric): - """``to_parameters(from_parameters(p)) == p`` exactly, however large ``p`` is. - - The amplitudes are this ansatz's parameters directly, so both directions are a plain re-indexing - and the round-trip is scale-free -- pinned here at a scale far beyond the other round-trip tests', - and to exact equality rather than a tolerance, since no arithmetic is performed on the values. - """ - norb = 5 - expected = UCC.num_parameters(norb, nocc, variant, antisymmetric=antisymmetric) - rng = np.random.default_rng(abs(hash((variant, antisymmetric))) % (2**32)) - params = rng.standard_normal(expected) * 20.0 - - gate = UCC.from_parameters(params, norb, nocc, variant, antisymmetric=antisymmetric) - - np.testing.assert_allclose(gate.to_parameters(), params, rtol=0, atol=0) + assert gate.uccsd_op.final_orbital_rotation is None -def test_ucc_from_parameters_wrong_length_raises(): - """A parameter vector of the wrong length is rejected.""" - with pytest.raises(ValueError, match="did not match the number expected"): - UCC.from_parameters(np.zeros(3), 4, 2, "restricted") +def test_ucc_unrestricted_reads_norb_off_the_operator(): + """An unrestricted operator gives a ``2 * norb``-mode gate with per-spin amplitudes.""" + gate = UCC(_unrestricted_op(2, 1, 2, 3, seed=19)) + assert gate.norb == 4 + assert len(gate.uccsd_op.t1) == 2 + assert len(gate.uccsd_op.t2) == 3 + + +def test_ucc_accepts_the_complex_operator_flavors(): + """The non-``Real`` operator flavors are accepted too.""" + t1 = np.zeros((1, 1)) + t2 = np.zeros((1, 1, 1, 1)) + assert UCC(ffsim.UCCSDOpRestricted(t1=t1, t2=t2)).num_modes == 4 + assert ( + UCC( + ffsim.UCCSDOpUnrestricted(t1=(t1, t1), t2=(t2, t2, t2)), + ).num_modes + == 4 + ) -def test_ucc_nocc_arity_is_validated(): - """``nocc`` must be a pair for the unrestricted variant and an integer otherwise.""" - with pytest.raises(ValueError, match="requires a \\(nocc_a, nocc_b\\) pair"): - UCC.num_parameters(4, 2, "unrestricted") - with pytest.raises(ValueError, match="only valid for the 'unrestricted' variant"): - UCC.num_parameters(4, (2, 2), "restricted") +def test_ucc_rejects_a_non_ffsim_operator(): + """Anything other than one of ffsim's four UCCSD operators is rejected up front.""" + with pytest.raises(TypeError, match="requires one of ffsim's UCCSD operators"): + UCC(np.zeros((2, 2))) def test_ucc_cluster_operator_generator_is_anti_hermitian(): @@ -237,8 +99,7 @@ def test_ucc_cluster_operator_generator_is_anti_hermitian(): This is exactly the property :meth:`.UCC._build_definition` relies on to express the ansatz as an :class:`.Evolution`, whose operator must be Hermitian for the evolution to be unitary. """ - t1, t2 = _restricted_amplitudes(2, 2, seed=13) - generator = UCC("restricted", t1, t2).cluster_operator() + generator = UCC(_restricted_op(2, 2, seed=13)).cluster_operator() assert (generator * 1j).is_hermitian() @@ -246,29 +107,26 @@ def test_ucc_cluster_operator_conserves_sector(): """The cluster generator conserves the particle number of each spin species. Every excitation replaces an occupied orbital with a virtual one *within* a spin sector, so the - generator must preserve both the total particle number and the z-component of spin -- the + generator must preserve both the total particle number and the z-component of spin, the condition :meth:`.Evolution._apply_unitary_placed_` enforces before simulating. """ - t1, t2 = _restricted_amplitudes(2, 2, seed=14) - gate = UCC("restricted", t1, t2) + gate = UCC(_restricted_op(2, 2, seed=14)) generator = gate.cluster_operator() assert generator.conserves_particle_number() assert generator.conserves_sector([gate.norb, gate.norb]) -@pytest.mark.parametrize("variant", ["restricted", "spinless"]) -def test_ucc_generator_groups_are_individually_hermitian(variant): - """Every group of the Hermitian generator is itself Hermitian -- a regression guard. +def test_ucc_generator_groups_are_individually_hermitian(): + """Every group of the Hermitian generator is itself Hermitian, a regression guard. :class:`.Evolution` decomposes group-by-group, so each group becomes one factor ``exp(-i H_k)`` of the product formula. A factor is unitary only if its ``H_k`` is Hermitian. Splitting ``i (T - T^dagger)`` *term*-by-term instead yields non-Hermitian factors (each - excitation is separated from its conjugate), which makes the synthesized circuit non-unitary -- + excitation is separated from its conjugate), which makes the synthesized circuit non-unitary: it does not even preserve the norm. The generator therefore pairs every excitation with its conjugate in a shared group, which this test locks in. """ - t1, t2 = _restricted_amplitudes(2, 2, seed=18) - generator = UCC(variant, t1, t2).cluster_operator() * 1j + generator = UCC(_restricted_op(2, 2, seed=18)).cluster_operator() * 1j assert generator.is_hermitian() assert generator.has_groups() @@ -279,16 +137,14 @@ def test_ucc_generator_groups_are_individually_hermitian(variant): def test_ucc_generator_groups_are_individually_hermitian_unrestricted(): """The conjugate-pairing group invariant also holds for the unrestricted variant.""" - t1, t2 = _unrestricted_amplitudes(2, 1, 2, 3, seed=19) - generator = UCC("unrestricted", t1, t2).cluster_operator() * 1j + generator = UCC(_unrestricted_op(2, 1, 2, 3, seed=19)).cluster_operator() * 1j assert generator.is_hermitian() for group in generator.split_out_groups(): assert group.is_hermitian() -@pytest.mark.parametrize("variant", ["restricted", "spinless"]) -def test_ucc_generator_group_layout_is_canonical(variant): +def test_ucc_generator_group_layout_is_canonical(): """The group layout is sorted by mode support, so it cannot depend on term-iteration order. The group index decides where in the product formula :class:`.Evolution` places that group's @@ -299,8 +155,7 @@ def test_ucc_generator_group_layout_is_canonical(variant): failure of ``test_ucc_trotterized_circuit_converges_to_the_exact_gate``. Asserting the layout is sorted pins that down without needing to vary the hash seed, which a single process cannot do. """ - t1, t2 = _restricted_amplitudes(2, 2, seed=23) - generator = UCC(variant, t1, t2).cluster_operator() + generator = UCC(_restricted_op(2, 2, seed=23)).cluster_operator() support_by_group: dict[int, tuple[int, ...]] = {} for (actions, _), group in zip( @@ -314,14 +169,6 @@ def test_ucc_generator_group_layout_is_canonical(variant): assert layout == sorted(layout), layout -def test_ucc_spinless_cluster_operator_acts_only_on_norb_modes(): - """The spinless generator stays within the single ``norb``-mode register (no spin offset).""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=15) - gate = UCC("spinless", t1, t2) - assert max(gate.cluster_operator().get_support()) < nocc + nvrt - - def test_ucc_definition_is_a_single_evolution(): """The gate's definition is one :class:`.Evolution` carrying the whole cluster generator. @@ -329,8 +176,7 @@ def test_ucc_definition_is_a_single_evolution(): simulation path exponentiate it exactly while leaving the Trotter decomposition to the transpiler. """ - t1, t2 = _restricted_amplitudes(1, 2, seed=16) - gate = UCC("restricted", t1, t2) + gate = UCC(_restricted_op(1, 2, seed=16)) circuit = FermionicCircuit(gate.num_modes) circuit.append(gate, circuit.modes) assert dict(circuit.decompose().count_ops()) == {"Evolution": 1} @@ -338,8 +184,7 @@ def test_ucc_definition_is_a_single_evolution(): def test_ucc_definition_decomposes_into_group_evolutions(): """Decomposing the definition's ``Evolution`` splits it into one evolution per group.""" - t1, t2 = _restricted_amplitudes(1, 2, seed=17) - gate = UCC("restricted", t1, t2) + gate = UCC(_restricted_op(1, 2, seed=17)) circuit = FermionicCircuit(gate.num_modes) circuit.append(gate, circuit.modes) counts = dict(circuit.decompose().decompose().count_ops()) @@ -351,11 +196,10 @@ def test_ucc_stays_hermitian_at_any_decomposition_depth(): The cluster generator carries conjugate-paired groups precisely so that each factor is Hermitian. Decomposing past those groups used to split them term-by-term, and an individual excitation is - *not* Hermitian on its own -- which produced complex Pauli coefficients that the transpiler + *not* Hermitian on its own, which produced complex Pauli coefficients that the transpiler rejected, and a non-normalized state vector in simulation. """ - t1, t2 = _restricted_amplitudes(1, 2, seed=18) - gate = UCC("restricted", t1, t2) + gate = UCC(_restricted_op(1, 2, seed=18)) circuit = FermionicCircuit(gate.num_modes) circuit.append(gate, circuit.modes) @@ -368,182 +212,11 @@ def test_ucc_stays_hermitian_at_any_decomposition_depth(): assert operator.is_hermitian(), f"non-Hermitian factor at reps={reps}" -def _antisymmetrize(t2): - """Projects a same-spin ``t2`` block onto the standard coupled-cluster antisymmetric subspace.""" - t2 = t2 - t2.transpose(1, 0, 2, 3) - return t2 - t2.transpose(0, 1, 3, 2) - - -def test_ucc_antisymmetric_num_parameters_is_smaller_unrestricted(): - """``antisymmetric=True`` shrinks the unrestricted parameter count to the smaller subspace. - - Both same-spin blocks drop from the full exchange-symmetric basis - (``n_pairs * (n_pairs + 1) / 2`` entries) to the strictly-upper-triangular - ``(i < j, a < b)`` basis, while the singles and the cross-spin ``t2ab`` block are untouched. - """ - norb, nocc_a, nocc_b = 4, 2, 1 - nvrt_a, nvrt_b = norb - nocc_a, norb - nocc_b - - full = UCC.num_parameters(norb, (nocc_a, nocc_b), "unrestricted") - reduced = UCC.num_parameters(norb, (nocc_a, nocc_b), "unrestricted", antisymmetric=True) - - unconstrained = nocc_a * nvrt_a + nocc_b * nvrt_b + nocc_a * nocc_b * nvrt_a * nvrt_b - n_pairs_a, n_pairs_b = nocc_a * nvrt_a, nocc_b * nvrt_b - assert full == unconstrained + sum(n * (n + 1) // 2 for n in (n_pairs_a, n_pairs_b)) - # only i < j, a < b survives: 1 occupied pair x 1 virtual pair for alpha, none at all for beta - assert reduced == unconstrained + 1 - assert reduced < full - - -def test_ucc_antisymmetric_num_parameters_is_smaller_spinless(): - """``antisymmetric=True`` shrinks the spinless parameter count to the smaller subspace.""" - norb, nocc = 4, 2 - nvrt = norb - nocc - - full = UCC.num_parameters(norb, nocc, "spinless") - reduced = UCC.num_parameters(norb, nocc, "spinless", antisymmetric=True) - - n_pairs = nocc * nvrt - assert full == nocc * nvrt + n_pairs * (n_pairs + 1) // 2 - assert reduced == nocc * nvrt + 1 # one (i < j) x (a < b) entry - assert reduced < full - - -@pytest.mark.parametrize( - ("variant", "norb", "nocc"), - [("unrestricted", 5, (3, 2)), ("spinless", 5, 3)], -) -def test_ucc_antisymmetric_parameters_round_trip(variant, norb, nocc): - """``to_parameters`` inverts ``from_parameters`` in the antisymmetric basis. - - The flag is stored on the instance, so ``to_parameters`` writes out the *same* (smaller) basis - ``from_parameters`` read -- a round-trip through the full exchange-symmetric basis would return a - longer vector. - """ - expected = UCC.num_parameters(norb, nocc, variant, antisymmetric=True) - rng = np.random.default_rng(20) - params = rng.standard_normal(expected) - - gate = UCC.from_parameters(params, norb, nocc, variant, antisymmetric=True) - assert gate.antisymmetric - assert len(gate.to_parameters()) == expected - np.testing.assert_allclose(gate.to_parameters(), params, atol=1e-12) - - -@pytest.mark.parametrize( - ("variant", "norb", "nocc"), - [("unrestricted", 5, (3, 2)), ("spinless", 5, 3)], -) -def test_ucc_antisymmetric_from_parameters_builds_antisymmetric_amplitudes(variant, norb, nocc): - """The amplitudes ``from_parameters(antisymmetric=True)`` builds really are antisymmetric. - - Construction and validation are deliberately independent: ``from_parameters`` expands each - independent entry into its four sign copies, and the constructor then re-checks the result. This - test pins down the property itself rather than trusting that composition. - """ - expected = UCC.num_parameters(norb, nocc, variant, antisymmetric=True) - rng = np.random.default_rng(21) - gate = UCC.from_parameters( - rng.standard_normal(expected), norb, nocc, variant, antisymmetric=True - ) - - blocks = [gate.t2[0], gate.t2[2]] if variant == "unrestricted" else [gate.t2] - for t2 in blocks: - np.testing.assert_allclose(t2, -t2.transpose(1, 0, 2, 3), atol=1e-12) - np.testing.assert_allclose(t2, -t2.transpose(0, 1, 3, 2), atol=1e-12) - - -def test_ucc_antisymmetric_rejects_non_antisymmetric_amplitudes(): - """A ``t2`` outside the antisymmetric subspace is rejected when the flag is set -- but not else. - - The exchange-symmetric amplitudes ``t2[i, j, a, b] == t2[j, i, b, a]`` used everywhere else are a - strictly larger family, so the *same* tensor must be accepted with the flag off. Asserting both - halves keeps the test from passing for the wrong reason (e.g. an unrelated shape error). - """ - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=22) - - with pytest.raises(ValueError, match="are not antisymmetric under occupied"): - UCC("spinless", t1, t2, antisymmetric=True) - - assert UCC("spinless", t1, t2).norb == nocc + nvrt - - -def test_ucc_antisymmetric_accepts_antisymmetric_amplitudes(): - """Genuinely antisymmetric amplitudes pass the opt-in validation.""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=23) - gate = UCC("spinless", t1, _antisymmetrize(t2), antisymmetric=True) - assert gate.antisymmetric - - -def test_ucc_antisymmetric_validates_both_unrestricted_same_spin_blocks(): - """Both unrestricted same-spin blocks are validated, and the cross-spin block is exempt. - - ``t2ab`` carries no antisymmetry (exchanging an alpha with a beta excitation is not an exchange of - identical operators), so leaving it unsymmetrized must not trip the check. - """ - (t1a, t1b), (t2aa, t2ab, t2bb) = _unrestricted_amplitudes(2, 2, 2, 2, seed=24) - - # a violation in either same-spin block is caught, named by block - with pytest.raises(ValueError, match="the t2aa amplitudes are not antisymmetric"): - UCC("unrestricted", (t1a, t1b), (t2aa, t2ab, _antisymmetrize(t2bb)), antisymmetric=True) - with pytest.raises(ValueError, match="the t2bb amplitudes are not antisymmetric"): - UCC("unrestricted", (t1a, t1b), (_antisymmetrize(t2aa), t2ab, t2bb), antisymmetric=True) - - # ... while the raw (unsymmetrized) cross-spin block is accepted - gate = UCC( - "unrestricted", - (t1a, t1b), - (_antisymmetrize(t2aa), t2ab, _antisymmetrize(t2bb)), - antisymmetric=True, - ) - assert gate.antisymmetric - - -def test_ucc_antisymmetric_atol_is_honored(): - """A violation below ``atol`` is tolerated, and the default tolerance is tight enough to catch it.""" - nocc, nvrt = 2, 2 - t1, t2 = _restricted_amplitudes(nocc, nvrt, seed=25) - t2 = _antisymmetrize(t2) - t2[0, 0, 0, 0] = 1e-6 # a small, deliberate antisymmetry violation - - with pytest.raises(ValueError, match="are not antisymmetric"): - UCC("spinless", t1, t2, antisymmetric=True) - assert UCC("spinless", t1, t2, antisymmetric=True, atol=1e-4).antisymmetric - - -@pytest.mark.parametrize("variant", ["restricted", UCC.Variant.RESTRICTED]) -def test_ucc_antisymmetric_is_refused_for_restricted(variant): - """``antisymmetric=True`` is refused for the restricted variant at every entry point. - - The restricted ``t2`` supplies both the same-spin *and* the alpha-beta amplitudes, and the - antisymmetry only applies to the former; honoring the flag would silently over-constrain the - cross-spin channel. Refusing explicitly is preferred over a flag that means something else here. - """ - t1, t2 = _restricted_amplitudes(2, 2, seed=26) - match = "not supported for the 'restricted' variant" - - with pytest.raises(ValueError, match=match): - UCC(variant, t1, t2, antisymmetric=True) - with pytest.raises(ValueError, match=match): - UCC.from_t_amplitudes(t2, t1=t1, variant=variant, antisymmetric=True) - with pytest.raises(ValueError, match=match): - UCC.num_parameters(4, 2, variant, antisymmetric=True) - with pytest.raises(ValueError, match=match): - UCC.from_parameters(np.zeros(1), 4, 2, variant, antisymmetric=True) - - -def test_ucc_antisymmetric_flag_does_not_change_the_operator(): - """The flag only restricts *which* amplitudes are allowed, never what they mean. - - Given amplitudes that satisfy the antisymmetry, the generator must be identical whether or not the - flag was set -- it selects a parameterization and switches on a validation, and must not sneak an - extra symmetrization or prefactor into the operator itself. - """ - t1, t2 = _restricted_amplitudes(2, 2, seed=27) - t2 = _antisymmetrize(t2) - - checked = UCC("spinless", t1, t2, antisymmetric=True).cluster_operator() - unchecked = UCC("spinless", t1, t2).cluster_operator() - assert checked.equiv(unchecked) +def test_ucc_final_orbital_rotation_is_appended_per_spin_sector(): + """A wrapped operator's final orbital rotation becomes a closing per-spin OrbitalRotation.""" + gate = UCC(_restricted_op(1, 1, seed=21, with_final=True)) + circuit = FermionicCircuit(gate.num_modes) + circuit.append(gate, circuit.modes) + counts = dict(circuit.decompose().count_ops()) + assert counts["Evolution"] == 1 + assert counts["OrbitalRotation"] == 2 diff --git a/tests/python/circuit/library/test_ucc_apply_unitary.py b/tests/python/circuit/library/test_ucc_apply_unitary.py index 882128940..f7e4375aa 100644 --- a/tests/python/circuit/library/test_ucc_apply_unitary.py +++ b/tests/python/circuit/library/test_ucc_apply_unitary.py @@ -16,8 +16,8 @@ ``UCCSDOpUnrestrictedReal``: we take the amplitudes from a random ffsim UCCSD operator, rebuild the ansatz with our :class:`.UCC` gate, and require the resulting state vectors to agree. This pins down both the mode convention (ffsim's interleaved ``(orb, spin)`` versus our block-spin register) and the -per-block prefactors of the cluster operator. The spinless variant has no ffsim counterpart, so it is -validated against a directly exponentiated cluster generator instead. +per-block prefactors of the cluster operator. A separate group of tests uses a directly exponentiated +cluster generator as the oracle, which is independent of the gate's ``Evolution``-based definition. """ from __future__ import annotations @@ -67,7 +67,7 @@ def test_ucc_restricted_matches_ffsim(): reference = ffsim.hartree_fock_state(norb, nelec) expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - gate = UCC("restricted", t1, t2) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) np.testing.assert_allclose(result, expected, atol=1e-10) @@ -88,7 +88,7 @@ def test_ucc_unrestricted_matches_ffsim(): reference = ffsim.hartree_fock_state(norb, nelec) expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - gate = UCC("unrestricted", t1, t2) + gate = UCC(ffsim.UCCSDOpUnrestrictedReal(t1=t1, t2=t2)) result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) np.testing.assert_allclose(result, expected, atol=1e-10) @@ -103,7 +103,7 @@ def test_ucc_restricted_doubles_only_matches_ffsim(): reference = ffsim.hartree_fock_state(norb, nelec) expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - gate = UCC.from_t_amplitudes(t2) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=np.zeros((t2.shape[0], t2.shape[2])), t2=t2)) result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) np.testing.assert_allclose(result, expected, atol=1e-10) @@ -122,7 +122,7 @@ def test_ucc_gate_apply_unitary_matches_ffsim(): reference = ffsim.hartree_fock_state(norb, nelec) expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - gate = UCC("restricted", t1, t2) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) result = gate._apply_unitary_(reference.copy(), norb, nelec, copy=True) np.testing.assert_allclose(result, expected, atol=1e-10) @@ -137,73 +137,26 @@ def test_ucc_gate_through_circuit_matches_ffsim(): reference = ffsim.hartree_fock_state(norb, nelec) expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - gate = UCC("restricted", t1, t2) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) circ = FermionicCircuit(2 * norb) circ.append(gate, circ.modes) result = ffsim.apply_unitary(reference, circ, norb=norb, nelec=nelec) np.testing.assert_allclose(result, expected, atol=1e-10) -def test_ucc_from_parameters_matches_ffsim_parameter_ordering(): - """UCC.from_parameters, fed ffsim's own to_parameters() vector, reproduces ffsim's state. +def test_ucc_matches_directly_exponentiated_generator(): + """The gate matches ``exp(T - T^dagger)`` applied via scipy directly. - This validates that the parameter *ordering* :meth:`.UCC.from_parameters` / - :meth:`.UCC.to_parameters` use matches ffsim's ``UCCSDOpRestrictedReal`` convention exactly, not - just internal round-trip self-consistency: a vector produced by ffsim is handed unmodified to our - gate. + The oracle here is the gate's own cluster generator exponentiated independently of the gate's + ``Evolution``-based definition, which confirms that the ``e^{T - T^dagger} == e^{-i H}`` rewrite + with ``H = i (T - T^dagger)`` carries the right sign. """ norb, nocc = 4, 2 - nelec = (nocc, nocc) - t1, t2 = _restricted_amplitudes(nocc, norb - nocc, seed=4200) - - ucc_op = ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2) - params = ucc_op.to_parameters() - assert len(params) == UCC.num_parameters(norb, nocc, "restricted") - - gate = UCC.from_parameters(params, norb, nocc, "restricted") - reference = ffsim.hartree_fock_state(norb, nelec) - expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) - np.testing.assert_allclose(result, expected, atol=1e-10) - - -def test_ucc_from_parameters_unrestricted_matches_ffsim_parameter_ordering(): - """The unrestricted parameter ordering matches ffsim's, including the t2ab block's position. - - ffsim orders the unrestricted doubles as ``t2aa, t2ab, t2bb`` -- with the cross-spin block - *between* the two same-spin blocks. Feeding ffsim's own vector in unmodified would silently swap - the ``t2ab``/``t2bb`` slices if that order were misread. - """ - norb = 4 - nocc_a, nocc_b = 2, 1 - nelec = (nocc_a, nocc_b) - t1, t2 = _unrestricted_amplitudes(nocc_a, nocc_b, norb - nocc_a, norb - nocc_b, seed=99) - - ucc_op = ffsim.UCCSDOpUnrestrictedReal(t1=t1, t2=t2) - params = ucc_op.to_parameters() - assert len(params) == UCC.num_parameters(norb, (nocc_a, nocc_b), "unrestricted") - - gate = UCC.from_parameters(params, norb, (nocc_a, nocc_b), "unrestricted") - reference = ffsim.hartree_fock_state(norb, nelec) - expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) - np.testing.assert_allclose(result, expected, atol=1e-10) - - -def test_ucc_spinless_matches_directly_exponentiated_generator(): - """The spinless gate matches ``exp(T - T^dagger)`` applied via scipy directly. - - ffsim has no spinless UCCSD operator, so the oracle here is the gate's own cluster generator - exponentiated independently of the gate's ``Evolution``-based definition. This confirms the - spinless variant stays on its single ``norb``-mode register (no spin offset) and that the - ``e^{T - T^dagger} == e^{-i H}`` rewrite with ``H = i (T - T^dagger)`` carries the right sign. - """ - norb, nocc = 4, 2 - nelec = 3 + nelec = (2, 2) t1, t2 = _restricted_amplitudes(nocc, norb - nocc, seed=11) - gate = UCC("spinless", t1, t2) - reference = ffsim.slater_determinant(norb, list(range(nelec))) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) + reference = ffsim.hartree_fock_state(norb, nelec) linop = ffsim.linear_operator(gate.cluster_operator(), norb, nelec) expected = scipy.sparse.linalg.expm_multiply(linop, reference, traceA=0.0) @@ -224,7 +177,7 @@ def test_ucc_restricted_matches_directly_exponentiated_generator(): nelec = (nocc, nocc) t1, t2 = _restricted_amplitudes(nocc, norb - nocc, seed=12) - gate = UCC("restricted", t1, t2) + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) reference = ffsim.hartree_fock_state(norb, nelec) linop = ffsim.linear_operator(gate.cluster_operator(), norb, nelec) @@ -237,61 +190,34 @@ def test_ucc_restricted_matches_directly_exponentiated_generator(): def test_ucc_gate_subset_placement_matches_global_embedding(): """A UCC placed on a subset of a larger register acts on its absolute (global) modes. - A spinless UCC on ``norb_local`` modes is appended to a ``norb_global``-mode register on the - non-identity subset ``placement``. The oracle is the *same* cluster generator relabeled onto those - global modes and exponentiated directly. This guards the ``_apply_unitary_placed_`` routing: with - the placement ignored, the gate would act on modes ``0..norb_local`` instead of ``placement``. + A 4-mode UCC is appended to a larger register on the non-identity subset ``placement``. The + oracle is the *same* cluster generator relabeled onto those global modes and exponentiated + directly. This guards the ``_apply_unitary_placed_`` routing: with the placement ignored, the + gate would act on the leading modes instead of ``placement``. """ - norb_local, norb_global = 2, 4 - nelec = 2 # spinless integer nelec on the global register - placement = [1, 3] # global modes the local UCC modes map onto (non-identity) + norb_global = 6 + nelec = (2, 1) + # the gate's 4 block-spin modes land on global alpha modes {1, 2} and beta modes {norb+0, norb+2} + placement = [1, 2, norb_global + 0, norb_global + 2] t1, t2 = _restricted_amplitudes(1, 1, seed=909) - gate = UCC("spinless", t1, t2) - assert gate.num_modes == norb_local + gate = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)) + assert gate.num_modes == 4 - reference = ffsim.slater_determinant(norb_global, placement) + reference = ffsim.hartree_fock_state(norb_global, nelec) # oracle: the same generator, relabeled from local modes onto the placed global modes relabeled = gate.cluster_operator().relabel_modes(placement) linop = ffsim.linear_operator(relabeled, norb_global, nelec) expected = scipy.sparse.linalg.expm_multiply(linop, reference, traceA=0.0) - circ = FermionicCircuit(norb_global) + circ = FermionicCircuit(2 * norb_global) circ.append(gate, [circ.modes[p] for p in placement]) result = ffsim.apply_unitary(reference, circ, norb=norb_global, nelec=nelec) np.testing.assert_allclose(result, expected, atol=1e-10) -def test_ucc_from_t_amplitudes_restricted_matches_ffsim_with_ccsd_amplitudes(): - """UCC built from real CCSD amplitudes matches ffsim's operator on the same amplitudes. - - Exercises the ansatz on physically meaningful amplitudes (a genuine CCSD calculation) rather than - random tensors, and confirms ``from_t_amplitudes`` threads ``t1``/``t2`` through unchanged. - """ - pyscf = pytest.importorskip("pyscf") - import pyscf.cc as _pyscf_cc - - mol = pyscf.gto.Mole() - mol.build( - atom=[["H", (0, 0, 0)], ["H", (0, 0, 0.74)]], basis="6-31g", symmetry="Dooh", verbose=0 - ) - scf = pyscf.scf.RHF(mol).run() - mol_data = ffsim.MolecularData.from_scf(scf) - norb, nelec = mol_data.norb, mol_data.nelec - ccsd = _pyscf_cc.CCSD(scf).run() - t1, t2 = ccsd.t1, ccsd.t2 - - reference = ffsim.hartree_fock_state(norb, nelec) - ucc_op = ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2) - expected = ffsim.apply_unitary(reference, ucc_op, norb=norb, nelec=nelec) - - gate = UCC.from_t_amplitudes(t2, t1=t1, variant="restricted") - result = ffsim.apply_unitary(reference, gate, norb=norb, nelec=nelec) - np.testing.assert_allclose(result, expected, atol=1e-10) - - def test_ucc_trotterized_circuit_converges_to_the_exact_gate(): """The synthesized (Trotterized) circuit approaches the exactly-applied gate as reps increase. @@ -305,12 +231,14 @@ def test_ucc_trotterized_circuit_converges_to_the_exact_gate(): t1, t2 = _restricted_amplitudes(nocc, norb - nocc, seed=17) reference = ffsim.hartree_fock_state(norb, nelec) - exact = ffsim.apply_unitary(reference, UCC("restricted", t1, t2), norb=norb, nelec=nelec) + exact = ffsim.apply_unitary( + reference, UCC(ffsim.UCCSDOpRestrictedReal(t1=t1, t2=t2)), norb=norb, nelec=nelec + ) errors = [] for n_reps in (1, 4, 16, 64): # n_reps repetitions of a 1/n_reps-strength ansatz, each Trotterized group-by-group - step = UCC("restricted", t1 / n_reps, t2 / n_reps) + step = UCC(ffsim.UCCSDOpRestrictedReal(t1=t1 / n_reps, t2=t2 / n_reps)) circ = FermionicCircuit(2 * norb) for _ in range(n_reps): circ.append(step, circ.modes)