diff --git a/cpp/src/qdk/chemistry/data/pauli_operator.cpp b/cpp/src/qdk/chemistry/data/pauli_operator.cpp index 72f39d394..a874de3b0 100644 --- a/cpp/src/qdk/chemistry/data/pauli_operator.cpp +++ b/cpp/src/qdk/chemistry/data/pauli_operator.cpp @@ -687,20 +687,26 @@ std::unique_ptr SumPauliOperatorExpression::clone() std::unique_ptr SumPauliOperatorExpression::distribute() const { - // Short-circuit: if already distributed, just clone - if (this->is_distributed()) { - return std::make_unique(*this); - } - auto result = std::make_unique(); result->reserve_capacity(terms_.size() * 2); - for (const auto& term : terms_) { - auto distributed_term = term->distribute(); - for (const auto& dist_term : distributed_term->get_terms()) { - result->add_term(dist_term->clone()); + // Recursively flatten nested sums and distribute product terms. + std::function flatten; + flatten = [&result, &flatten](const SumPauliOperatorExpression* sum) { + for (const auto& term : sum->get_terms()) { + if (const auto* nested_sum = term->as_sum_expression()) { + // Flatten nested sums recursively + flatten(nested_sum); + } else { + auto distributed_term = term->distribute(); + for (const auto& dist_term : distributed_term->get_terms()) { + result->add_term(dist_term->clone()); + } + } } - } + }; + + flatten(this); return result; } diff --git a/docs/source/user/comprehensive/model_hamiltonians.rst b/docs/source/user/comprehensive/model_hamiltonians.rst index 0c1015d75..025995dfa 100644 --- a/docs/source/user/comprehensive/model_hamiltonians.rst +++ b/docs/source/user/comprehensive/model_hamiltonians.rst @@ -20,10 +20,11 @@ Fermionic models * **Pariser-Parr-Pople (PPP)** — extends Hubbard with long-range intersite Coulomb interactions Spin models - Operate on spin-½ degrees of freedom and produce :class:`~qdk_chemistry.data.QubitHamiltonian` objects expressed as sums of Pauli operators. + Operate on spin-S degrees of freedom and produce :class:`~qdk_chemistry.data.QubitHamiltonian` objects expressed as sums of Pauli operators. + Support arbitrary spin quantum numbers and mixed-spin lattices. - * **Heisenberg** — anisotropic spin-spin coupling with external magnetic fields - * **Ising** — special case of Heisenberg with ZZ coupling and transverse X field + * **Heisenberg** — anisotropic spin-spin coupling with biquadratic terms and external magnetic fields + * **Ising** — special case of Heisenberg with :math:`S^z S^z` coupling and transverse :math:`S^x` field All model Hamiltonian builders take a :doc:`LatticeGraph ` as their first argument, which defines the site connectivity and hopping structure. For a brief description of the available model Hamiltonian builders, see the table below. @@ -52,11 +53,11 @@ For a more detailed description of each model Hamiltonian and their parameters, * - ``create_heisenberg_hamiltonian`` - Spin - QubitHamiltonian - - Anisotropic spin coupling + external fields + - Anisotropic spin-S coupling + biquadratic + fields * - ``create_ising_hamiltonian`` - Spin - QubitHamiltonian - - ZZ coupling + transverse X field + - :math:`S^z S^z` coupling + transverse :math:`S^x` field Fermionic models ---------------- @@ -216,60 +217,202 @@ The ``pairwise_potential`` function accepts a user-defined callable ``func(i, j, Spin models ----------- +The spin model builders construct Hamiltonians for quantum spin systems on a lattice. +They use physical spin operators :math:`\hat{S}^\alpha` with eigenvalues +:math:`-S, -S+1, \ldots, S`, and produce :class:`~qdk_chemistry.data.QubitHamiltonian` +objects expressed as sums of Pauli operators. + +Spin models support **arbitrary spin quantum numbers** (S = 1/2, 1, 3/2, 2, ...), +including **mixed-spin lattices** where different sites carry different spin values. +Each site is encoded into :math:`\lceil \log_2(2S+1) \rceil` qubits. + +.. note:: + + The spin model builders use physical spin operators :math:`\hat{S}^\alpha` + (not Pauli matrices :math:`\hat{\sigma}^\alpha`). + For spin-1/2, :math:`\hat{S}^\alpha = \hat{\sigma}^\alpha / 2`, so coupling + constants differ by a factor of 4 (two-body) or 2 (one-body) compared to the + Pauli-matrix convention. + .. _model-heisenberg: Heisenberg model ~~~~~~~~~~~~~~~~ -The anisotropic Heisenberg model describes spin-½ particles interacting on a lattice with external magnetic fields: +The anisotropic Heisenberg model describes spin-S particles interacting on a lattice +with optional biquadratic coupling and external magnetic fields: .. math:: \hat{H}_\text{Heisenberg} = \sum_{\langle i,j \rangle} w_{ij}\,\bigl( - J_x^{ij}\,\hat{\sigma}_i^x \hat{\sigma}_j^x - + J_y^{ij}\,\hat{\sigma}_i^y \hat{\sigma}_j^y - + J_z^{ij}\,\hat{\sigma}_i^z \hat{\sigma}_j^z + J_x^{ij}\,\hat{S}_i^x \hat{S}_j^x + + J_y^{ij}\,\hat{S}_i^y \hat{S}_j^y + + J_z^{ij}\,\hat{S}_i^z \hat{S}_j^z \bigr) + + \sum_{\langle i,j \rangle} w_{ij}\,J_\text{bq}^{ij}\, + (\hat{\mathbf{S}}_i \cdot \hat{\mathbf{S}}_j)^2 + \sum_i \bigl( - h_x^{i}\,\hat{\sigma}_i^x - + h_y^{i}\,\hat{\sigma}_i^y - + h_z^{i}\,\hat{\sigma}_i^z + h_x^{i}\,\hat{S}_i^x + + h_y^{i}\,\hat{S}_i^y + + h_z^{i}\,\hat{S}_i^z \bigr) -where :math:`J_x, J_y, J_z` are the spin-spin coupling constants, :math:`h_x, h_y, h_z` are external magnetic field components, and :math:`w_{ij}` is the edge weight from the lattice adjacency matrix. -Each qubit corresponds to a lattice site. +where :math:`J_x, J_y, J_z` are the bilinear spin-spin coupling constants, +:math:`J_\text{bq}` is the biquadratic coupling constant, +:math:`h_x, h_y, h_z` are external magnetic field components, +and :math:`w_{ij}` is the edge weight from the lattice adjacency matrix. -.. tab:: Python API - - .. literalinclude:: ../../_static/examples/python/model_hamiltonians.py - :language: python - :start-after: # start-cell-create-heisenberg - :end-before: # end-cell-create-heisenberg +The ``spins`` parameter controls the spin quantum number per site (default 0.5). +Pass a scalar for uniform spin or a list for mixed-spin lattices. Special cases of the Heisenberg model include: - **Isotropic (XXX)**: :math:`J_x = J_y = J_z` - **XXZ**: :math:`J_x = J_y \ne J_z` - **XY**: :math:`J_z = 0` +- **Bilinear-biquadratic (BLBQ)**: :math:`J_x = J_y = J_z,\; J_\text{bq} \ne 0` (requires S > 1/2 for non-trivial effects) + +.. tab:: Python API + + .. code-block:: python + + from qdk_chemistry.data import LatticeGraph + from qdk_chemistry.utils.model_hamiltonians import create_heisenberg_hamiltonian + + lattice = LatticeGraph.chain(4) + + # Spin-1/2 (default) — isotropic Heisenberg + qh_half = create_heisenberg_hamiltonian(lattice, jx=1.0, jy=1.0, jz=1.0) + + # Spin-3/2 — no penalty needed (dim 4 = 2^2) + qh_3half = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.5 + ) + + # Spin-1 — requires penalty (dim 3 ≠ 2^k) + qh_one = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.0, penalty_strength=10.0 + ) + + # Bilinear-biquadratic spin-1 model + qh_blbq = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, + j_biquadratic=0.5, spins=1.0, penalty_strength=10.0 + ) + + # Mixed-spin lattice + qh_mixed = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, + spins=[0.5, 1.5, 0.5, 1.5] + ) .. _model-ising: Ising model ~~~~~~~~~~~ -The transverse-field Ising model is a special case of the Heisenberg model with ZZ coupling and a transverse X field: +The transverse-field Ising model is a special case of the Heisenberg model with +:math:`\hat{S}^z \hat{S}^z` coupling and a transverse :math:`\hat{S}^x` field: .. math:: - \hat{H}_\text{Ising} = \sum_{\langle i,j \rangle} w_{ij}\,J^{ij}\,\hat{\sigma}_i^z \hat{\sigma}_j^z - + \sum_i h^{i}\,\hat{\sigma}_i^x + \hat{H}_\text{Ising} = \sum_{\langle i,j \rangle} w_{ij}\,J^{ij}\,\hat{S}_i^z \hat{S}_j^z + + \sum_i h^{i}\,\hat{S}_i^x + +The ``spins`` and ``penalty_strength`` parameters work identically to the Heisenberg builder. .. tab:: Python API - .. literalinclude:: ../../_static/examples/python/model_hamiltonians.py - :language: python - :start-after: # start-cell-create-ising - :end-before: # end-cell-create-ising + .. code-block:: python + + from qdk_chemistry.data import LatticeGraph + from qdk_chemistry.utils.model_hamiltonians import create_ising_hamiltonian + + lattice = LatticeGraph.chain(4) + + # Spin-1/2 Ising (default) + qh = create_ising_hamiltonian(lattice, j=1.0, h=0.5) + + # Spin-3/2 Ising + qh_3half = create_ising_hamiltonian(lattice, j=1.0, h=0.5, spins=1.5) + +.. _qubit-encoding: + +Qubit encoding for higher spins +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Each lattice site with spin quantum number :math:`S` is encoded into +:math:`n = \lceil \log_2(2S+1) \rceil` qubits, giving :math:`2^n` computational +basis states per site. When :math:`2S+1` is a power of 2 (e.g. S = 1/2, 3/2, 7/2), +all computational basis states are physical and no special treatment is needed. + +When :math:`2S+1` is **not** a power of 2 (e.g. S = 1, 2, 5/2), the encoding +introduces :math:`2^n - (2S+1)` unphysical states per site. The standard approach +is to add a **penalty Hamiltonian** that raises the energy of these states: + +.. math:: + + \hat{H}_\text{total} = \hat{H}_\text{model} + + \lambda \sum_i \hat{P}_{\text{unphys},i} + +where :math:`\hat{P}_{\text{unphys},i}` projects onto the unphysical subspace of +site :math:`i` and :math:`\lambda` is a large positive constant (the +``penalty_strength`` parameter). A sufficiently large :math:`\lambda` ensures that +unphysical states lie well above the physical spectrum. + +.. list-table:: Qubit encoding summary + :header-rows: 1 + :widths: 10 15 15 15 20 + + * - Spin S + - States (2S+1) + - Qubits + - 2\ :sup:`n` + - Penalty needed? + * - 1/2 + - 2 + - 1 + - 2 + - No + * - 1 + - 3 + - 2 + - 4 + - **Yes** (1 unphysical) + * - 3/2 + - 4 + - 2 + - 4 + - No + * - 2 + - 5 + - 3 + - 8 + - **Yes** (3 unphysical) + * - 5/2 + - 6 + - 3 + - 8 + - **Yes** (2 unphysical) + * - 7/2 + - 8 + - 3 + - 8 + - No + +The builders raise a ``ValueError`` if any site requires a penalty and +``penalty_strength`` is not set, preventing accidental use of unphysical states. + +The :class:`~qdk_chemistry.utils.spin_operators.SpinEncoding` helper tracks the +qubit-to-site mapping for mixed-spin lattices: + +.. code-block:: python + + from qdk_chemistry.utils.spin_operators import SpinEncoding + + enc = SpinEncoding([0.5, 1.5, 1.0, 0.5]) + print(enc.total_qubits) # 6 (1 + 2 + 2 + 1) + print(enc.site_qubits(1)) # range(1, 3) — site 1 uses qubits 1 and 2 Parameter flexibility --------------------- diff --git a/python/src/qdk_chemistry/utils/__init__.py b/python/src/qdk_chemistry/utils/__init__.py index 29eb6f9c8..f89cb2f48 100644 --- a/python/src/qdk_chemistry/utils/__init__.py +++ b/python/src/qdk_chemistry/utils/__init__.py @@ -8,7 +8,7 @@ from qdk_chemistry._core.utils import Logger, compute_valence_space_parameters, rotate_orbitals from qdk_chemistry.utils.enum import CaseInsensitiveStrEnum -from . import model_hamiltonians +from . import model_hamiltonians, spin_operators __all__ = [ "CaseInsensitiveStrEnum", @@ -16,4 +16,5 @@ "compute_valence_space_parameters", "model_hamiltonians", "rotate_orbitals", + "spin_operators", ] diff --git a/python/src/qdk_chemistry/utils/model_hamiltonians.py b/python/src/qdk_chemistry/utils/model_hamiltonians.py index fac2f63eb..505f1323d 100644 --- a/python/src/qdk_chemistry/utils/model_hamiltonians.py +++ b/python/src/qdk_chemistry/utils/model_hamiltonians.py @@ -5,6 +5,10 @@ # Licensed under the MIT License. See LICENSE.txt in the project root for license information. # -------------------------------------------------------------------------------------------- +from __future__ import annotations + +from typing import Sequence + import numpy as np from qdk_chemistry._core.utils.model_hamiltonians import ( @@ -18,6 +22,7 @@ to_site_param, ) from qdk_chemistry.data import LatticeGraph, PauliOperator, QubitHamiltonian +from qdk_chemistry.utils.spin_operators import SpinEncoding __all__ = [ "create_heisenberg_hamiltonian", @@ -48,38 +53,58 @@ def create_heisenberg_hamiltonian( hx: np.ndarray | float = 0.0, hy: np.ndarray | float = 0.0, hz: np.ndarray | float = 0.0, + spins: float | Sequence[float] = 0.5, + j_biquadratic: np.ndarray | float = 0.0, + penalty_strength: float | None = None, ) -> QubitHamiltonian: r"""Create the anisotropic Heisenberg model Hamiltonian on a lattice. + For spin-1/2 (default), the Hamiltonian is: + .. math:: H = \sum_{\langle i,j \rangle} w_{ij}\,\bigl[ - J_x^{ij}\,\sigma_i^x \sigma_j^x - + J_y^{ij}\,\sigma_i^y \sigma_j^y - + J_z^{ij}\,\sigma_i^z \sigma_j^z + J_x^{ij}\,S_i^x S_j^x + + J_y^{ij}\,S_i^y S_j^y + + J_z^{ij}\,S_i^z S_j^z \bigr] + + \sum_{\langle i,j \rangle} w_{ij}\,J_\mathrm{bq}^{ij}\, + (\mathbf{S}_i \cdot \mathbf{S}_j)^2 + \sum_i \bigl[ - h_x^{i}\,\sigma_i^x - + h_y^{i}\,\sigma_i^y - + h_z^{i}\,\sigma_i^z + h_x^{i}\,S_i^x + h_y^{i}\,S_i^y + h_z^{i}\,S_i^z \bigr] - where :math:`w_{ij}` is the edge weight from the lattice adjacency matrix. + where :math:`w_{ij}` is the edge weight from the lattice adjacency matrix and + :math:`S_i^\alpha` are spin-S operators with eigenvalues :math:`-S, -S+1, \ldots, S`. + + Each lattice site is encoded into :math:`\lceil \log_2(2S+1) \rceil` qubits. + When :math:`2S+1` is not a power of 2, the encoding introduces unphysical + computational basis states. Use ``penalty_strength`` to add an energy penalty + that pushes these states out of the low-energy spectrum. - Each qubit corresponds to a lattice site. + .. note:: + + This function uses physical spin operators :math:`S^\alpha` (not Pauli matrices + :math:`\sigma^\alpha`). For spin-1/2, :math:`S^\alpha = \sigma^\alpha / 2`. Args: graph: Lattice graph defining the connectivity. - jx: Coupling constant for XX interactions. Scalar (uniform) or ``(n, n)`` array for per-pair values. - jy: Coupling constant for YY interactions (same format as *jx*). - jz: Coupling constant for ZZ interactions (same format as *jx*). + jx: Coupling constant for :math:`S^x S^x` interactions. Scalar (uniform) or ``(n, n)`` array for per-pair values. + jy: Coupling constant for :math:`S^y S^y` interactions (same format as *jx*). + jz: Coupling constant for :math:`S^z S^z` interactions (same format as *jx*). hx: External magnetic field in the x direction. Scalar or length-n array. Defaults to 0. hy: External magnetic field in the y direction. Defaults to 0. hz: External magnetic field in the z direction. Defaults to 0. + spins: Spin quantum number per site. Scalar for uniform spin, or a sequence of length ``n`` for mixed-spin lattices. Defaults to 0.5. + j_biquadratic: Biquadratic coupling constant for :math:`(\mathbf{S}_i \cdot \mathbf{S}_j)^2`. Scalar or ``(n, n)`` array. Defaults to 0. + penalty_strength: Energy penalty applied to unphysical qubit states. Required when any site has :math:`2S+1 \neq 2^k`. Set to ``None`` (default) for power-of-2 encodings. Returns: QubitHamiltonian: The Heisenberg model as a qubit Hamiltonian. + Raises: + ValueError: If any site has unphysical states and ``penalty_strength`` is not set. + """ if not graph.is_symmetric: raise ValueError("Lattice graph must be symmetric for a valid Hamiltonian.") @@ -94,54 +119,101 @@ def create_heisenberg_hamiltonian( hy_vec = to_site_param(hy, graph, "hy") hz_vec = to_site_param(hz, graph, "hz") + # Build qubit layout from spin assignments + encoding = SpinEncoding(spins, num_sites=n) + + # Check for unphysical states + if encoding.has_unphysical_states and penalty_strength is None: + raise ValueError( + "Some sites have 2S+1 ≠ 2^k, introducing unphysical qubit states. " + "Set penalty_strength to a positive value to penalize these states, " + "e.g. penalty_strength=10.0." + ) + + # Build site operators + site_ops = [encoding.site_operators(i) for i in range(n)] + h = 0.0 * PauliOperator.I(0) # seed with zero expression + # Parse biquadratic coupling + has_biquadratic = not (isinstance(j_biquadratic, (int, float)) and j_biquadratic == 0.0) + if has_biquadratic: + jbq_mat = to_pair_param(j_biquadratic, graph, "j_biquadratic") + # Two-body interaction terms (each edge counted once: i < j), scaled by the lattice edge weight. for i in range(n): for j in range(i + 1, n): edge_weight = adj[i, j] if edge_weight == 0.0: continue + + # Bilinear terms: Jα * Sᵢα * Sⱼα if jx_mat[i, j] != 0.0: - h = h + jx_mat[i, j] * edge_weight * PauliOperator.X(i) * PauliOperator.X(j) + h = h + jx_mat[i, j] * edge_weight * site_ops[i].sx * site_ops[j].sx if jy_mat[i, j] != 0.0: - h = h + jy_mat[i, j] * edge_weight * PauliOperator.Y(i) * PauliOperator.Y(j) + h = h + jy_mat[i, j] * edge_weight * site_ops[i].sy * site_ops[j].sy if jz_mat[i, j] != 0.0: - h = h + jz_mat[i, j] * edge_weight * PauliOperator.Z(i) * PauliOperator.Z(j) + h = h + jz_mat[i, j] * edge_weight * site_ops[i].sz * site_ops[j].sz + + # Biquadratic term: Jbq * (Sᵢ·Sⱼ)² + if has_biquadratic and jbq_mat[i, j] != 0.0: + dot_ij = ( + site_ops[i].sx * site_ops[j].sx + + site_ops[i].sy * site_ops[j].sy + + site_ops[i].sz * site_ops[j].sz + ) + h = h + jbq_mat[i, j] * edge_weight * dot_ij * dot_ij # Single-body field terms for i in range(n): if hx_vec[i] != 0.0: - h = h + hx_vec[i] * PauliOperator.X(i) + h = h + hx_vec[i] * site_ops[i].sx if hy_vec[i] != 0.0: - h = h + hy_vec[i] * PauliOperator.Y(i) + h = h + hy_vec[i] * site_ops[i].sy if hz_vec[i] != 0.0: - h = h + hz_vec[i] * PauliOperator.Z(i) + h = h + hz_vec[i] * site_ops[i].sz + + # Penalty terms for unphysical states + if penalty_strength is not None and penalty_strength != 0.0: + for i in range(n): + pen = site_ops[i].penalty_projector() + if pen is not None: + h = h + penalty_strength * pen - return _pauli_expr_to_qubit_hamiltonian(h, n) + return _pauli_expr_to_qubit_hamiltonian(h, encoding.total_qubits) def create_ising_hamiltonian( graph: LatticeGraph, j: np.ndarray | float, h: np.ndarray | float = 0.0, + spins: float | Sequence[float] = 0.5, + penalty_strength: float | None = None, ) -> QubitHamiltonian: r"""Create the Ising model Hamiltonian on a lattice. .. math:: - H = \sum_{\langle i,j \rangle} w_{ij}\,J^{ij}\,\sigma_i^z \sigma_j^z - + \sum_i h^{i}\,\sigma_i^x + H = \sum_{\langle i,j \rangle} w_{ij}\,J^{ij}\,S_i^z S_j^z + + \sum_i h^{i}\,S_i^x + + where :math:`w_{ij}` is the edge weight from the lattice adjacency matrix and + :math:`S_i^\alpha` are spin-S operators. - where :math:`w_{ij}` is the edge weight from the lattice adjacency matrix. + This is a special case of the Heisenberg model with only :math:`S^z S^z` + coupling and a transverse :math:`S^x` field. Args: graph: Lattice graph defining the connectivity. - j: Coupling constant for ZZ interactions. Scalar or ``(n, n)`` array. + j: Coupling constant for :math:`S^z S^z` interactions. Scalar or ``(n, n)`` array. h: Transverse field strength (x direction). Scalar or length-n array. Defaults to 0. + spins: Spin quantum number per site. Scalar or sequence. Defaults to 0.5. + penalty_strength: Energy penalty for unphysical qubit states. See :func:`create_heisenberg_hamiltonian`. Returns: QubitHamiltonian: The Ising model as a qubit Hamiltonian. """ - return create_heisenberg_hamiltonian(graph, jx=0.0, jy=0.0, jz=j, hx=h) + return create_heisenberg_hamiltonian( + graph, jx=0.0, jy=0.0, jz=j, hx=h, spins=spins, penalty_strength=penalty_strength + ) diff --git a/python/src/qdk_chemistry/utils/spin_operators.py b/python/src/qdk_chemistry/utils/spin_operators.py new file mode 100644 index 000000000..a9313c822 --- /dev/null +++ b/python/src/qdk_chemistry/utils/spin_operators.py @@ -0,0 +1,417 @@ +"""Spin-S operator construction and qubit encoding for higher-spin lattice models. + +This module provides utilities for constructing spin-S operators (Sˣ, Sʸ, Sᶻ) +decomposed into Pauli operator expressions, enabling the construction of +arbitrary-spin Heisenberg and Ising Hamiltonians on qubit hardware. + +Each lattice site with spin quantum number S is encoded into ⌈log₂(2S+1)⌉ qubits. +When 2S+1 is not a power of 2, the embedding introduces unphysical computational +basis states that must be handled via energy penalties. +""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +from __future__ import annotations + +import itertools +import math +from functools import lru_cache +from typing import Sequence + +import numpy as np + +from qdk_chemistry.data import PauliOperator + +__all__ = [ + "SpinEncoding", + "SpinSOperators", +] + +# Pauli matrices used in decomposition (column-major-friendly complex128). +_PAULI = { + "I": np.eye(2, dtype=np.complex128), + "X": np.array([[0, 1], [1, 0]], dtype=np.complex128), + "Y": np.array([[0, -1j], [1j, 0]], dtype=np.complex128), + "Z": np.array([[1, 0], [0, -1]], dtype=np.complex128), +} + +_PAULI_LABELS = ("I", "X", "Y", "Z") + + +def _validate_spin(spin: float) -> int: + """Validate spin and return 2*spin as an integer. + + Args: + spin: Spin quantum number (must be a positive integer or half-integer). + + Returns: + Integer value of 2*spin. + + Raises: + ValueError: If spin is not a positive integer or half-integer. + + """ + two_s = round(2 * spin) + if two_s < 1 or not math.isclose(2 * spin, two_s): + raise ValueError(f"spin must be a positive integer or half-integer, got {spin}") + return two_s + + +def _spin_matrices(spin: float) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Build the (2S+1)-dimensional spin matrices Sˣ, Sʸ, Sᶻ. + + Uses standard angular momentum algebra: + - Sᶻ is diagonal with eigenvalues m = S, S-1, ..., -S (descending) + - S⁺|m⟩ = √(S(S+1) - m(m+1)) |m+1⟩ + - Sˣ = (S⁺ + S⁻)/2, Sʸ = (S⁺ - S⁻)/(2i) + + Args: + spin: Spin quantum number S. + + Returns: + Tuple (Sˣ, Sʸ, Sᶻ) as complex128 matrices of shape (2S+1, 2S+1). + + """ + two_s = _validate_spin(spin) + dim = two_s + 1 + s = spin + + # Sᶻ: diagonal, eigenvalues descending from S to -S + m_values = np.arange(s, -s - 0.5, -1.0) + sz = np.diag(m_values).astype(np.complex128) + + # S⁺ (raising operator): superdiagonal + s_plus = np.zeros((dim, dim), dtype=np.complex128) + for i in range(dim - 1): + m = m_values[i + 1] # m of the ket (lower state) + s_plus[i, i + 1] = np.sqrt(s * (s + 1) - m * (m + 1)) + + s_minus = s_plus.T.copy() + + sx = (s_plus + s_minus) / 2.0 + sy = (s_plus - s_minus) / 2.0j + + return sx, sy, sz + + +def _embed_matrix(mat: np.ndarray, num_qubits: int) -> np.ndarray: + """Embed a (d x d) matrix into a (2ⁿ x 2ⁿ) matrix by zero-padding. + + The physical states occupy the first d rows/columns; unphysical states are zero. + + Args: + mat: Matrix of shape (d, d) where d <= 2^num_qubits. + num_qubits: Number of qubits for the embedding. + + Returns: + Matrix of shape (2^num_qubits, 2^num_qubits). + + """ + dim_full = 1 << num_qubits + d = mat.shape[0] + if d == dim_full: + return mat.copy() + result = np.zeros((dim_full, dim_full), dtype=np.complex128) + result[:d, :d] = mat + return result + + +@lru_cache(maxsize=64) +def _pauli_decomposition_cached( + spin: float, +) -> tuple[dict[str, complex], dict[str, complex], dict[str, complex]]: + """Compute the Pauli decomposition of spin-S operators (cached by spin value). + + Decomposes Sˣ, Sʸ, Sᶻ into sums of Pauli strings: + M = Σ_P (Tr(P · M) / 2ⁿ) · P + + Args: + spin: Spin quantum number. + + Returns: + Tuple of dicts (sx_terms, sy_terms, sz_terms) mapping Pauli string labels to complex coefficients. + + """ + two_s = _validate_spin(spin) + dim = two_s + 1 + num_qubits = math.ceil(math.log2(dim)) if dim > 1 else 1 + dim_full = 1 << num_qubits + + sx, sy, sz = _spin_matrices(spin) + sx_emb = _embed_matrix(sx, num_qubits) + sy_emb = _embed_matrix(sy, num_qubits) + sz_emb = _embed_matrix(sz, num_qubits) + + atol = 1e-14 + + results = [] + for mat in (sx_emb, sy_emb, sz_emb): + terms: dict[str, complex] = {} + for pauli_labels in itertools.product(_PAULI_LABELS, repeat=num_qubits): + # Build the n-qubit Pauli matrix via Kronecker product + p_mat = _PAULI[pauli_labels[0]] + for lbl in pauli_labels[1:]: + p_mat = np.kron(p_mat, _PAULI[lbl]) + + coeff = np.trace(p_mat @ mat) / dim_full + if abs(coeff) > atol: + # Clean up near-real and near-imaginary coefficients + if abs(coeff.imag) < atol: + coeff = complex(coeff.real, 0.0) + elif abs(coeff.real) < atol: + coeff = complex(0.0, coeff.imag) + label = "".join(pauli_labels) + terms[label] = coeff + results.append(terms) + + return tuple(results) # type: ignore[return-value] + + +@lru_cache(maxsize=64) +def _penalty_decomposition_cached(spin: float) -> dict[str, complex] | None: + """Compute the Pauli decomposition of the penalty projector for unphysical states. + + The projector is P_unphys = Σ_{k=dim}^{2ⁿ-1} |k⟩⟨k| = I - Σ_{k=0}^{dim-1} |k⟩⟨k|. + + Args: + spin: Spin quantum number. + + Returns: + Dict mapping Pauli string labels to coefficients, or None if dim is a power of 2. + + """ + two_s = _validate_spin(spin) + dim = two_s + 1 + num_qubits = math.ceil(math.log2(dim)) if dim > 1 else 1 + dim_full = 1 << num_qubits + + if dim == dim_full: + return None + + # Build projector onto unphysical subspace + projector = np.zeros((dim_full, dim_full), dtype=np.complex128) + for k in range(dim, dim_full): + projector[k, k] = 1.0 + + atol = 1e-14 + terms: dict[str, complex] = {} + for pauli_labels in itertools.product(_PAULI_LABELS, repeat=num_qubits): + p_mat = _PAULI[pauli_labels[0]] + for lbl in pauli_labels[1:]: + p_mat = np.kron(p_mat, _PAULI[lbl]) + + coeff = np.trace(p_mat @ projector) / dim_full + if abs(coeff) > atol: + if abs(coeff.imag) < atol: + coeff = complex(coeff.real, 0.0) + elif abs(coeff.real) < atol: + coeff = complex(0.0, coeff.imag) + terms["".join(pauli_labels)] = coeff + + return terms + + +def _build_pauli_expr(terms: dict[str, complex], qubit_offset: int) -> PauliOperator: + """Convert a dict of {pauli_label: coeff} into a PauliOperator expression. + + The label convention from the Pauli decomposition uses Kronecker-product order + (label[0] = MSB qubit). PauliOperator uses qubit indices where lower index = LSB, + so we reverse the mapping: label[k] → qubit (offset + num_qubits - 1 - k). + + Args: + terms: Mapping from Pauli string labels (e.g. "XZ") to complex coefficients. + qubit_offset: Global qubit index offset for this site. + + Returns: + PauliOperator expression representing the sum. + + """ + _pauli_factory = {"I": PauliOperator.I, "X": PauliOperator.X, "Y": PauliOperator.Y, "Z": PauliOperator.Z} + + num_qubits = len(next(iter(terms))) # all labels have the same length + expr = None + for label, coeff in terms.items(): + # Build the product of single-qubit Paulis for this term. + # label[0] acts on the MSB → highest qubit index. + term = coeff + for k, ch in enumerate(label): + term = term * _pauli_factory[ch](qubit_offset + num_qubits - 1 - k) + expr = term if expr is None else expr + term + return expr + + +class SpinSOperators: + """Spin-S operators for a single lattice site, decomposed into Pauli operators. + + Provides Sˣ, Sʸ, Sᶻ as :class:`~qdk_chemistry.data.PauliOperator` expressions + acting on the qubits assigned to this site. Operator matrices follow the standard + angular momentum algebra with eigenvalues of Sᶻ being -S, -S+1, ..., S. + + Args: + spin: Spin quantum number S (positive integer or half-integer). + qubit_offset: Index of the first qubit for this site in the global qubit register. + + Raises: + ValueError: If spin is not a positive integer or half-integer. + + """ + + def __init__(self, spin: float, qubit_offset: int) -> None: + two_s = _validate_spin(spin) + self._spin = spin + self._two_s = two_s + self._dim = two_s + 1 + self._num_qubits = math.ceil(math.log2(self._dim)) if self._dim > 1 else 1 + self._qubit_offset = qubit_offset + + sx_terms, sy_terms, sz_terms = _pauli_decomposition_cached(spin) + self._sx = _build_pauli_expr(sx_terms, qubit_offset) + self._sy = _build_pauli_expr(sy_terms, qubit_offset) + self._sz = _build_pauli_expr(sz_terms, qubit_offset) + + @property + def spin(self) -> float: + """Spin quantum number S.""" + return self._spin + + @property + def dim(self) -> int: + """Number of physical spin states (2S+1).""" + return self._dim + + @property + def num_qubits(self) -> int: + """Number of qubits used to encode this site.""" + return self._num_qubits + + @property + def qubit_offset(self) -> int: + """Index of the first qubit for this site.""" + return self._qubit_offset + + @property + def has_unphysical_states(self) -> bool: + """Whether the qubit encoding has unphysical (unused) computational basis states.""" + return self._dim < (1 << self._num_qubits) + + @property + def sx(self) -> PauliOperator: + """Sˣ operator as a PauliOperator expression.""" + return self._sx + + @property + def sy(self) -> PauliOperator: + """Sʸ operator as a PauliOperator expression.""" + return self._sy + + @property + def sz(self) -> PauliOperator: + """Sᶻ operator as a PauliOperator expression.""" + return self._sz + + def penalty_projector(self) -> PauliOperator | None: + """Projector onto the unphysical subspace as a PauliOperator expression. + + Returns: + PauliOperator expression, or None if the encoding has no unphysical states. + + """ + terms = _penalty_decomposition_cached(self._spin) + if terms is None: + return None + return _build_pauli_expr(terms, self._qubit_offset) + + +class SpinEncoding: + """Qubit layout for a lattice with per-site spin assignments. + + Tracks which qubits encode which lattice site, supporting mixed-spin lattices + where different sites may have different spin quantum numbers. + + Args: + spins: Per-site spin quantum numbers. A scalar assigns the same spin to all sites (use with ``num_sites``). A sequence assigns spins individually. + num_sites: Number of lattice sites (required when ``spins`` is a scalar). + + Raises: + ValueError: If spins is a scalar and num_sites is not provided, or if any spin value is invalid. + + Examples: + >>> enc = SpinEncoding(0.5, num_sites=4) # uniform spin-1/2 + >>> enc = SpinEncoding([0.5, 1.5, 0.5, 1.5]) # mixed spins + + """ + + def __init__(self, spins: float | Sequence[float], num_sites: int | None = None) -> None: + if isinstance(spins, (int, float)): + if num_sites is None: + raise ValueError("num_sites is required when spins is a scalar") + self._spins = [float(spins)] * num_sites + else: + self._spins = [float(s) for s in spins] + + # Validate all spins and compute qubit layout + self._offsets: list[int] = [] + self._num_qubits_per_site: list[int] = [] + offset = 0 + for s in self._spins: + two_s = _validate_spin(s) + dim = two_s + 1 + nq = math.ceil(math.log2(dim)) if dim > 1 else 1 + self._offsets.append(offset) + self._num_qubits_per_site.append(nq) + offset += nq + self._total_qubits = offset + + @property + def num_sites(self) -> int: + """Number of lattice sites.""" + return len(self._spins) + + @property + def spins(self) -> list[float]: + """Per-site spin quantum numbers.""" + return list(self._spins) + + @property + def total_qubits(self) -> int: + """Total number of qubits across all sites.""" + return self._total_qubits + + @property + def has_unphysical_states(self) -> bool: + """Whether any site has unphysical states in its qubit encoding.""" + return any( + int(2 * s + 1) < (1 << nq) for s, nq in zip(self._spins, self._num_qubits_per_site) + ) + + def site_qubits(self, site: int) -> range: + """Return the qubit index range for a given lattice site. + + Args: + site: Lattice site index (0-based). + + Returns: + Range of qubit indices assigned to this site. + + """ + start = self._offsets[site] + return range(start, start + self._num_qubits_per_site[site]) + + def site_operators(self, site: int) -> SpinSOperators: + """Create a SpinSOperators instance for a given lattice site. + + Args: + site: Lattice site index (0-based). + + Returns: + SpinSOperators for the specified site. + + """ + return SpinSOperators(self._spins[site], self._offsets[site]) + + def __repr__(self) -> str: + spin_str = ", ".join(f"{s}" for s in self._spins) + return f"SpinEncoding(spins=[{spin_str}], total_qubits={self._total_qubits})" diff --git a/python/tests/test_model_hamiltonians.py b/python/tests/test_model_hamiltonians.py index e7ceaf502..5086a16ba 100644 --- a/python/tests/test_model_hamiltonians.py +++ b/python/tests/test_model_hamiltonians.py @@ -123,18 +123,20 @@ def test_ising_chain(self): j = 1.0 h = 0.5 + # Convention: S^α = σ^α/2 for spin-1/2, so two-body coefficients are J/4 + # and one-body coefficients are h/2. expected = {} - # edges + # edges: J * S^z_i S^z_j → J/4 * Z_i Z_j for edge in [(0, 1), (1, 2), (2, 3)]: pauli = ["I"] * n pauli[n - 1 - edge[0]] = "Z" pauli[n - 1 - edge[1]] = "Z" - expected["".join(pauli)] = j - # sites + expected["".join(pauli)] = j / 4 + # sites: h * S^x_i → h/2 * X_i for i in range(n): pauli = ["I"] * n pauli[n - 1 - i] = "X" - expected["".join(pauli)] = h + expected["".join(pauli)] = h / 2 # scalar qh = create_ising_hamiltonian(lattice, j=j, h=h) @@ -160,21 +162,23 @@ def test_ising_chain(self): h_vec[2] = 0.9 qh_modified = create_ising_hamiltonian(lattice, j=j_mat, h=h_vec) terms_mod = _get_terms_dict(qh_modified) - assert terms_mod["IIZZ"] == pytest.approx(2.5, abs=float_comparison_absolute_tolerance) - assert terms_mod["IZZI"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_mod["ZZII"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_mod["IIIX"] == pytest.approx(0.5, abs=float_comparison_absolute_tolerance) - assert terms_mod["IXII"] == pytest.approx(0.9, abs=float_comparison_absolute_tolerance) + assert terms_mod["IIZZ"] == pytest.approx(2.5 / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IZZI"] == pytest.approx(1.0 / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["ZZII"] == pytest.approx(1.0 / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IIIX"] == pytest.approx(0.5 / 2, abs=float_comparison_absolute_tolerance) + assert terms_mod["IXII"] == pytest.approx(0.9 / 2, abs=float_comparison_absolute_tolerance) - # weighted edges + # weighted edges: edge weight 0.5 scales two-body terms lattice_w = LatticeGraph.chain(3, t=0.5) qh_w = create_ising_hamiltonian(lattice_w, j=2.0, h=1.0) terms_w = _get_terms_dict(qh_w) - assert terms_w["IZZ"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_w["ZZI"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_w["IIX"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_w["IXI"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) - assert terms_w["XII"] == pytest.approx(1.0, abs=float_comparison_absolute_tolerance) + # 2.0 * 0.5 / 4 = 0.25 + assert terms_w["IZZ"] == pytest.approx(0.25, abs=float_comparison_absolute_tolerance) + assert terms_w["ZZI"] == pytest.approx(0.25, abs=float_comparison_absolute_tolerance) + # fields are not edge-weighted: 1.0 / 2 = 0.5 + assert terms_w["IIX"] == pytest.approx(0.5, abs=float_comparison_absolute_tolerance) + assert terms_w["IXI"] == pytest.approx(0.5, abs=float_comparison_absolute_tolerance) + assert terms_w["XII"] == pytest.approx(0.5, abs=float_comparison_absolute_tolerance) def test_heisenberg_chain(self): n = 4 @@ -187,6 +191,8 @@ def test_heisenberg_chain(self): hy = -2.0 hz = -3.0 + # Convention: S^α = σ^α/2 for spin-1/2, so two-body coefficients are J/4 + # and one-body coefficients are h/2. expected = {} # edges for edge in edges: @@ -194,13 +200,13 @@ def test_heisenberg_chain(self): ps = ["I"] * n ps[n - 1 - edge[0]] = pauli_char ps[n - 1 - edge[1]] = pauli_char - expected["".join(ps)] = j_val + expected["".join(ps)] = j_val / 4 # sites for i in range(n): for pauli_char, h_val in [("X", hx), ("Y", hy), ("Z", hz)]: ps = ["I"] * n ps[n - 1 - i] = pauli_char - expected["".join(ps)] = h_val + expected["".join(ps)] = h_val / 2 # scalar qh = create_heisenberg_hamiltonian(lattice, jx=jx, jy=jy, jz=jz, hx=hx, hy=hy, hz=hz) @@ -235,16 +241,16 @@ def test_heisenberg_chain(self): lattice, jx=jx_mat, jy=jy_mat, jz=jz_mat, hx=hx_vec, hy=hy_vec, hz=hz_vec ) terms_mod = _get_terms_dict(qh_modified) - assert terms_mod["IIXX"] == pytest.approx(2.5, abs=float_comparison_absolute_tolerance) - assert terms_mod["IIYY"] == pytest.approx(jy, abs=float_comparison_absolute_tolerance) - assert terms_mod["IIZZ"] == pytest.approx(jz, abs=float_comparison_absolute_tolerance) - assert terms_mod["IXXI"] == pytest.approx(jx, abs=float_comparison_absolute_tolerance) - assert terms_mod["IYYI"] == pytest.approx(jy, abs=float_comparison_absolute_tolerance) - assert terms_mod["IZZI"] == pytest.approx(0.3, abs=float_comparison_absolute_tolerance) - assert terms_mod["IIIX"] == pytest.approx(hx, abs=float_comparison_absolute_tolerance) - assert terms_mod["IXII"] == pytest.approx(0.9, abs=float_comparison_absolute_tolerance) - - # weighted edges + assert terms_mod["IIXX"] == pytest.approx(2.5 / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IIYY"] == pytest.approx(jy / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IIZZ"] == pytest.approx(jz / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IXXI"] == pytest.approx(jx / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IYYI"] == pytest.approx(jy / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IZZI"] == pytest.approx(0.3 / 4, abs=float_comparison_absolute_tolerance) + assert terms_mod["IIIX"] == pytest.approx(hx / 2, abs=float_comparison_absolute_tolerance) + assert terms_mod["IXII"] == pytest.approx(0.9 / 2, abs=float_comparison_absolute_tolerance) + + # weighted edges: edge weight 2.0 scales two-body terms lattice_w = LatticeGraph.chain(3, t=2.0) qh_w = create_heisenberg_hamiltonian(lattice_w, jx=1.0, jy=1.0, jz=1.0) terms_w = _get_terms_dict(qh_w) @@ -253,4 +259,5 @@ def test_heisenberg_chain(self): ps = ["I"] * 3 ps[3 - 1 - edge[0]] = pauli_char ps[3 - 1 - edge[1]] = pauli_char - assert terms_w["".join(ps)] == pytest.approx(2.0, abs=float_comparison_absolute_tolerance) + # 1.0 * 2.0 / 4 = 0.5 + assert terms_w["".join(ps)] == pytest.approx(0.5, abs=float_comparison_absolute_tolerance) diff --git a/python/tests/test_spin_operators.py b/python/tests/test_spin_operators.py new file mode 100644 index 000000000..c6811b44f --- /dev/null +++ b/python/tests/test_spin_operators.py @@ -0,0 +1,358 @@ +"""Tests for spin-S operator construction and higher-spin model Hamiltonians.""" + +# -------------------------------------------------------------------------------------------- +# Copyright (c) Microsoft Corporation. All rights reserved. +# Licensed under the MIT License. See LICENSE.txt in the project root for license information. +# -------------------------------------------------------------------------------------------- + +import math + +import numpy as np +import pytest + +from qdk_chemistry.data import LatticeGraph, QubitHamiltonian +from qdk_chemistry.utils.model_hamiltonians import create_heisenberg_hamiltonian, create_ising_hamiltonian +from qdk_chemistry.utils.pauli_matrix import pauli_to_dense_matrix +from qdk_chemistry.utils.spin_operators import SpinEncoding, SpinSOperators, _spin_matrices + + +class TestSpinMatrices: + """Test angular momentum matrix construction.""" + + @pytest.mark.parametrize("spin", [0.5, 1.0, 1.5, 2.0, 2.5, 3.0]) + def test_commutation_relations(self, spin): + """[Sˣ, Sʸ] = iSᶻ (and cyclic permutations).""" + sx, sy, sz = _spin_matrices(spin) + atol = 1e-12 + + assert np.allclose(sx @ sy - sy @ sx, 1j * sz, atol=atol) + assert np.allclose(sy @ sz - sz @ sy, 1j * sx, atol=atol) + assert np.allclose(sz @ sx - sx @ sz, 1j * sy, atol=atol) + + @pytest.mark.parametrize("spin", [0.5, 1.0, 1.5, 2.0]) + def test_hermiticity(self, spin): + """Spin matrices must be Hermitian.""" + sx, sy, sz = _spin_matrices(spin) + atol = 1e-14 + assert np.allclose(sx, sx.conj().T, atol=atol) + assert np.allclose(sy, sy.conj().T, atol=atol) + assert np.allclose(sz, sz.conj().T, atol=atol) + + @pytest.mark.parametrize("spin", [0.5, 1.0, 1.5, 2.0]) + def test_casimir(self, spin): + """S² = Sˣ² + Sʸ² + Sᶻ² = S(S+1)·I.""" + sx, sy, sz = _spin_matrices(spin) + s_squared = sx @ sx + sy @ sy + sz @ sz + expected = spin * (spin + 1) * np.eye(int(2 * spin + 1)) + assert np.allclose(s_squared, expected, atol=1e-12) + + def test_spin_half_matches_pauli(self): + """For S=1/2, spin operators are Pauli/2.""" + sx, sy, sz = _spin_matrices(0.5) + assert np.allclose(sx, np.array([[0, 0.5], [0.5, 0]])) + assert np.allclose(sy, np.array([[0, -0.5j], [0.5j, 0]])) + assert np.allclose(sz, np.array([[0.5, 0], [0, -0.5]])) + + def test_spin_one_eigenvalues(self): + """Sᶻ for S=1 has eigenvalues {-1, 0, 1}.""" + _, _, sz = _spin_matrices(1.0) + eigenvalues = np.sort(np.linalg.eigvalsh(sz)) + assert np.allclose(eigenvalues, [-1.0, 0.0, 1.0], atol=1e-14) + + +class TestSpinSOperators: + """Test Pauli decomposition of spin-S operators.""" + + @pytest.mark.parametrize("spin", [0.5, 1.0, 1.5, 2.0]) + def test_pauli_decomposition_roundtrip(self, spin): + """Reconstruct matrix from Pauli decomposition and compare to original.""" + dim = int(2 * spin + 1) + num_qubits = math.ceil(math.log2(dim)) if dim > 1 else 1 + total_qubits = num_qubits + + ops = SpinSOperators(spin, qubit_offset=0) + + for op_name, op_expr in [("sx", ops.sx), ("sy", ops.sy), ("sz", ops.sz)]: + simplified = op_expr.simplify() + terms = simplified.to_canonical_terms(total_qubits) + pauli_strings = [t[1][::-1] for t in terms] + coefficients = np.array([complex(t[0]) for t in terms]) + + mat_reconstructed = pauli_to_dense_matrix(pauli_strings, coefficients) + + # Compare against the original spin matrix (embedded) + sx_orig, sy_orig, sz_orig = _spin_matrices(spin) + orig = {"sx": sx_orig, "sy": sy_orig, "sz": sz_orig}[op_name] + dim_full = 1 << num_qubits + embedded = np.zeros((dim_full, dim_full), dtype=np.complex128) + embedded[:dim, :dim] = orig + + assert np.allclose(mat_reconstructed, embedded, atol=1e-12), ( + f"Pauli decomposition roundtrip failed for {op_name} with spin={spin}" + ) + + def test_spin_half_gives_pauli(self): + """S=1/2 decomposition gives σ/2 coefficients.""" + ops = SpinSOperators(0.5, qubit_offset=0) + sx_terms = ops.sx.simplify().to_canonical_terms(1) + sy_terms = ops.sy.simplify().to_canonical_terms(1) + sz_terms = ops.sz.simplify().to_canonical_terms(1) + + # Sˣ = 0.5 * X + assert len(sx_terms) == 1 + assert sx_terms[0][1] == "X" + assert abs(complex(sx_terms[0][0]) - 0.5) < 1e-14 + + # Sʸ = 0.5 * Y + assert len(sy_terms) == 1 + assert sy_terms[0][1] == "Y" + assert abs(complex(sy_terms[0][0]) - 0.5) < 1e-14 + + # Sᶻ = 0.5 * Z + assert len(sz_terms) == 1 + assert sz_terms[0][1] == "Z" + assert abs(complex(sz_terms[0][0]) - 0.5) < 1e-14 + + def test_qubit_offset(self): + """Operators with offset act on the correct qubits.""" + ops = SpinSOperators(0.5, qubit_offset=3) + sz_terms = ops.sz.simplify().to_canonical_terms(5) + # Should be Z on qubit 3, I everywhere else + assert len(sz_terms) == 1 + # to_canonical_terms returns little-endian: index k = qubit k + label_le = sz_terms[0][1] + assert label_le[3] == "Z" + for k in range(5): + if k != 3: + assert label_le[k] == "I" + + @pytest.mark.parametrize("spin", [1.0, 2.0, 2.5]) + def test_has_unphysical_states(self, spin): + """Non-power-of-2 dims flag unphysical states.""" + ops = SpinSOperators(spin, qubit_offset=0) + dim = int(2 * spin + 1) + num_qubits = math.ceil(math.log2(dim)) + expected = dim < (1 << num_qubits) + assert ops.has_unphysical_states == expected + + def test_no_unphysical_states_power_of_two(self): + """S=1/2 (dim=2), S=3/2 (dim=4), S=7/2 (dim=8) have no unphysical states.""" + for spin in [0.5, 1.5, 3.5]: + ops = SpinSOperators(spin, qubit_offset=0) + assert not ops.has_unphysical_states + assert ops.penalty_projector() is None + + @pytest.mark.parametrize("spin", [1.0, 2.0, 2.5]) + def test_penalty_projector(self, spin): + """Penalty projector projects onto unphysical subspace.""" + dim = int(2 * spin + 1) + num_qubits = math.ceil(math.log2(dim)) + dim_full = 1 << num_qubits + + ops = SpinSOperators(spin, qubit_offset=0) + pen = ops.penalty_projector() + assert pen is not None + + # Reconstruct penalty matrix + pen_terms = pen.simplify().to_canonical_terms(num_qubits) + pauli_strings = [t[1][::-1] for t in pen_terms] + coefficients = np.array([complex(t[0]) for t in pen_terms]) + pen_mat = pauli_to_dense_matrix(pauli_strings, coefficients) + + # Should be projector onto states dim..dim_full-1 + expected = np.zeros((dim_full, dim_full), dtype=np.complex128) + for k in range(dim, dim_full): + expected[k, k] = 1.0 + + assert np.allclose(pen_mat, expected, atol=1e-12) + + def test_invalid_spin_raises(self): + """Invalid spin values raise ValueError.""" + with pytest.raises(ValueError): + SpinSOperators(0.3, qubit_offset=0) + with pytest.raises(ValueError): + SpinSOperators(-0.5, qubit_offset=0) + with pytest.raises(ValueError): + SpinSOperators(0.0, qubit_offset=0) + + +class TestSpinEncoding: + """Test qubit layout for mixed-spin lattices.""" + + def test_uniform_spin_half(self): + enc = SpinEncoding(0.5, num_sites=4) + assert enc.num_sites == 4 + assert enc.total_qubits == 4 + assert enc.spins == [0.5, 0.5, 0.5, 0.5] + assert enc.site_qubits(0) == range(0, 1) + assert enc.site_qubits(3) == range(3, 4) + assert not enc.has_unphysical_states + + def test_uniform_spin_three_half(self): + enc = SpinEncoding(1.5, num_sites=3) + assert enc.total_qubits == 6 # 2 qubits per site + assert enc.site_qubits(0) == range(0, 2) + assert enc.site_qubits(1) == range(2, 4) + assert enc.site_qubits(2) == range(4, 6) + assert not enc.has_unphysical_states # dim=4 = 2^2 + + def test_mixed_spins(self): + enc = SpinEncoding([0.5, 1.5, 1.0, 0.5]) + assert enc.num_sites == 4 + # site 0: S=1/2 → 1 qubit, site 1: S=3/2 → 2 qubits, + # site 2: S=1 → 2 qubits, site 3: S=1/2 → 1 qubit + assert enc.total_qubits == 6 + assert enc.site_qubits(0) == range(0, 1) + assert enc.site_qubits(1) == range(1, 3) + assert enc.site_qubits(2) == range(3, 5) + assert enc.site_qubits(3) == range(5, 6) + assert enc.has_unphysical_states # site 2 has S=1 (dim=3, 2 qubits=4 states) + + def test_scalar_requires_num_sites(self): + with pytest.raises(ValueError): + SpinEncoding(1.0) + + def test_site_operators_returns_correct_type(self): + enc = SpinEncoding(1.5, num_sites=2) + ops = enc.site_operators(0) + assert isinstance(ops, SpinSOperators) + assert ops.spin == 1.5 + assert ops.qubit_offset == 0 + + +class TestHigherSpinHamiltonians: + """Test Heisenberg/Ising builders with higher spin values.""" + + def test_spin_three_half_heisenberg_hermitian(self): + """Spin-3/2 Heisenberg on a 2-site chain produces a Hermitian QubitHamiltonian.""" + lattice = LatticeGraph.chain(2) + qh = create_heisenberg_hamiltonian(lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.5) + assert isinstance(qh, QubitHamiltonian) + assert qh.num_qubits == 4 # 2 sites × 2 qubits/site + assert qh.is_hermitian() + + def test_spin_one_heisenberg_hermitian(self): + """Spin-1 Heisenberg on a 3-site chain is Hermitian and requires penalty.""" + lattice = LatticeGraph.chain(3) + qh = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.0, penalty_strength=10.0 + ) + assert isinstance(qh, QubitHamiltonian) + assert qh.num_qubits == 6 # 3 sites × 2 qubits/site + assert qh.is_hermitian() + + def test_spin_one_requires_penalty(self): + """Spin-1 (dim=3, not power of 2) raises ValueError without penalty_strength.""" + lattice = LatticeGraph.chain(2) + with pytest.raises(ValueError, match="penalty_strength"): + create_heisenberg_hamiltonian(lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.0) + + def test_spin_three_half_no_penalty_needed(self): + """Spin-3/2 (dim=4 = 2^2) does not require penalty.""" + lattice = LatticeGraph.chain(2) + # Should not raise + qh = create_heisenberg_hamiltonian(lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.5) + assert isinstance(qh, QubitHamiltonian) + + def test_spin_one_heisenberg_matrix_correctness(self): + """Verify spin-1 Heisenberg 2-site matrix against direct construction.""" + lattice = LatticeGraph.chain(2) + jx, jy, jz = 1.0, 1.0, 1.0 + + qh = create_heisenberg_hamiltonian( + lattice, jx=jx, jy=jy, jz=jz, spins=1.0, penalty_strength=100.0 + ) + + # Build the matrix from the QubitHamiltonian + pauli_strs, coeffs = qh.pauli_strings, qh.coefficients + pauli_labels = ["".join(ps) for ps in pauli_strs] + h_mat = pauli_to_dense_matrix(pauli_labels, coeffs) + + # Build expected matrix from spin-1 operators directly + sx, sy, sz = _spin_matrices(1.0) + I3 = np.eye(3, dtype=np.complex128) + + # Physical subspace is 3×3 per site = 9-dim, embedded in 4×4 per site = 16-dim + # S_1^α ⊗ I + I ⊗ S_2^α for the full 4×4 ⊗ 4×4 space + def embed(m): + result = np.zeros((4, 4), dtype=np.complex128) + result[:3, :3] = m + return result + + sx_emb = embed(sx) + sy_emb = embed(sy) + sz_emb = embed(sz) + I4 = np.eye(4, dtype=np.complex128) + + h_expected = ( + jx * np.kron(sx_emb, sx_emb) + + jy * np.kron(sy_emb, sy_emb) + + jz * np.kron(sz_emb, sz_emb) + ) + + # Add penalty for unphysical states + pen1 = np.zeros((4, 4), dtype=np.complex128) + pen1[3, 3] = 1.0 + h_expected += 100.0 * np.kron(pen1, I4) + h_expected += 100.0 * np.kron(I4, pen1) + + # Compare only the physical subspace (first 3×3 block per site = rows/cols 0-8) + assert np.allclose(h_mat, h_expected, atol=1e-10) + + def test_mixed_spin_lattice(self): + """Mixed spin lattice (S=1/2 and S=3/2) produces valid Hamiltonian.""" + lattice = LatticeGraph.chain(3) + spins = [0.5, 1.5, 0.5] # site 0: 1 qubit, site 1: 2 qubits, site 2: 1 qubit + + qh = create_heisenberg_hamiltonian(lattice, jx=1.0, jy=1.0, jz=1.0, spins=spins) + assert isinstance(qh, QubitHamiltonian) + assert qh.num_qubits == 4 # 1 + 2 + 1 + assert qh.is_hermitian() + + def test_ising_with_higher_spin(self): + """Ising model with spin-3/2.""" + lattice = LatticeGraph.chain(2) + qh = create_ising_hamiltonian(lattice, j=1.0, h=0.5, spins=1.5) + assert isinstance(qh, QubitHamiltonian) + assert qh.num_qubits == 4 + assert qh.is_hermitian() + + def test_biquadratic_heisenberg(self): + """Bilinear-biquadratic model for spin-1.""" + lattice = LatticeGraph.chain(2) + qh = create_heisenberg_hamiltonian( + lattice, + jx=1.0, + jy=1.0, + jz=1.0, + j_biquadratic=0.5, + spins=1.0, + penalty_strength=10.0, + ) + assert isinstance(qh, QubitHamiltonian) + assert qh.is_hermitian() + + # The biquadratic term should add additional Pauli strings beyond the bilinear ones + qh_no_bq = create_heisenberg_hamiltonian( + lattice, jx=1.0, jy=1.0, jz=1.0, spins=1.0, penalty_strength=10.0 + ) + # Biquadratic model should have at least as many terms (typically more) + assert len(qh.pauli_strings) >= len(qh_no_bq.pauli_strings) + + def test_biquadratic_spin_half_eigenvalues(self): + """For spin-1/2, (S⃗ᵢ·S⃗ⱼ)² has eigenvalues 9/16 (singlet) and 1/16 (triplet).""" + lattice = LatticeGraph.chain(2) + + # Pure biquadratic model: H = (S⃗₁·S⃗₂)² + qh = create_heisenberg_hamiltonian( + lattice, jx=0.0, jy=0.0, jz=0.0, j_biquadratic=1.0, spins=0.5 + ) + + labels = ["".join(ps) for ps in qh.pauli_strings] + mat = pauli_to_dense_matrix(labels, qh.coefficients) + eigenvalues = np.sort(np.linalg.eigvalsh(mat.real)) + + # Singlet: (S⃗·S⃗)² = (-3/4)² = 9/16 + # Triplet (3-fold): (S⃗·S⃗)² = (1/4)² = 1/16 + expected = np.sort([9 / 16, 1 / 16, 1 / 16, 1 / 16]) + assert np.allclose(eigenvalues, expected, atol=1e-12)