"""
heisenberg — Anisotropic Heisenberg (XYZ) model
=================================================
Provides ``HeisenbergHamiltonian``, the anisotropic Heisenberg spin-½
model on a one-dimensional chain with optional periodic boundary
conditions and arbitrary external fields.
Configurations are binary vectors of length ``num_spins`` where
``0 ≡ spin-up (|↑⟩)`` and ``1 ≡ spin-down (|↓⟩)``.
"""
from __future__ import annotations
import numpy as np
import torch
from qvartools.hamiltonians.hamiltonian import Hamiltonian
__all__ = [
"HeisenbergHamiltonian",
]
# ---------------------------------------------------------------------------
# Helper: ensure field is a 1-D tensor
# ---------------------------------------------------------------------------
def _to_field_tensor(
value: float | int | np.ndarray | torch.Tensor,
num_spins: int,
name: str,
) -> torch.Tensor:
"""Convert a scalar or array-like field to a float64 tensor.
Parameters
----------
value : float, int, np.ndarray, or torch.Tensor
Scalar (broadcast to all sites) or per-site array.
num_spins : int
Number of spin sites.
name : str
Name of the field (for error messages).
Returns
-------
torch.Tensor
Shape ``(num_spins,)``, dtype ``float64``.
Raises
------
ValueError
If *value* is array-like and its length does not match *num_spins*.
"""
if isinstance(value, (int, float)):
return torch.full((num_spins,), float(value), dtype=torch.float64)
t = torch.as_tensor(value, dtype=torch.float64).flatten()
if t.shape[0] != num_spins:
raise ValueError(
f"Field '{name}' has length {t.shape[0]}, expected {num_spins}."
)
return t
# ===================================================================
# HeisenbergHamiltonian
# ===================================================================
[docs]
class HeisenbergHamiltonian(Hamiltonian):
r"""Anisotropic Heisenberg (XYZ) model on a 1-D chain.
.. math::
H = \sum_{\langle i,j \rangle}
\bigl( J_x\, S^x_i S^x_j + J_y\, S^y_i S^y_j
+ J_z\, S^z_i S^z_j \bigr)
+ \sum_i \bigl( h^x_i S^x_i + h^y_i S^y_i + h^z_i S^z_i \bigr)
where *S* operators are spin-½ operators (with eigenvalues ±½) and the
sum runs over nearest-neighbour pairs on a chain.
Convention: ``config[i] = 0`` → spin-up (S^z = +½),
``config[i] = 1`` → spin-down (S^z = −½).
Parameters
----------
num_spins : int
Number of spin-½ sites.
Jx : float, optional
XX coupling constant (default ``1.0``).
Jy : float, optional
YY coupling constant (default ``1.0``).
Jz : float, optional
ZZ coupling constant (default ``1.0``).
h_x : float or array-like, optional
External field in the x-direction (default ``0``).
h_y : float or array-like, optional
External field in the y-direction (default ``0``).
h_z : float or array-like, optional
External field in the z-direction (default ``0``).
periodic : bool, optional
Whether to use periodic boundary conditions (default ``True``).
Examples
--------
>>> import torch
>>> H = HeisenbergHamiltonian(num_spins=4, Jx=1.0, Jy=1.0, Jz=1.0)
>>> config = torch.tensor([0, 1, 0, 1])
>>> H.diagonal_element(config)
tensor(...)
>>> energy, state = H.exact_ground_state()
"""
def __init__(
self,
num_spins: int,
Jx: float = 1.0,
Jy: float = 1.0,
Jz: float = 1.0,
h_x: float | np.ndarray | torch.Tensor = 0.0,
h_y: float | np.ndarray | torch.Tensor = 0.0,
h_z: float | np.ndarray | torch.Tensor = 0.0,
periodic: bool = True,
) -> None:
super().__init__(num_sites=num_spins, local_dim=2)
self.Jx: float = float(Jx)
self.Jy: float = float(Jy)
self.Jz: float = float(Jz)
self.periodic: bool = periodic
self._h_x = _to_field_tensor(h_x, num_spins, "h_x")
self._h_y = _to_field_tensor(h_y, num_spins, "h_y")
self._h_z = _to_field_tensor(h_z, num_spins, "h_z")
# Precompute neighbour list
self._neighbours: list[tuple[int, int]] = []
n_bonds = num_spins if periodic else num_spins - 1
for i in range(n_bonds):
j = (i + 1) % num_spins
self._neighbours.append((i, j))
# ------------------------------------------------------------------
# Diagonal
# ------------------------------------------------------------------
[docs]
def diagonal_element(self, config: torch.Tensor) -> torch.Tensor:
"""Compute ⟨config|H|config⟩ (ZZ interaction + Z field).
Parameters
----------
config : torch.Tensor
Spin configuration, shape ``(num_sites,)``, entries in ``{0, 1}``.
Returns
-------
torch.Tensor
Scalar diagonal energy.
"""
return self.diagonal_elements_batch(config.unsqueeze(0)).squeeze(0)
[docs]
def diagonal_elements_batch(self, configs: torch.Tensor) -> torch.Tensor:
"""Vectorised diagonal elements for a batch of configurations.
Parameters
----------
configs : torch.Tensor
Shape ``(batch, num_sites)``.
Returns
-------
torch.Tensor
Shape ``(batch,)``, dtype ``float64``.
"""
configs_f = configs.to(dtype=torch.float64)
# S^z_i = 0.5 - config[i] (config 0 → +0.5, config 1 → -0.5)
sz = 0.5 - configs_f # (batch, num_sites)
energy = torch.zeros(
configs_f.shape[0], dtype=torch.float64, device=configs.device
)
# ZZ interaction
for i, j in self._neighbours:
energy = energy + self.Jz * sz[:, i] * sz[:, j]
# Z field: h_z_i * S^z_i
energy = energy + torch.einsum("bi,i->b", sz, self._h_z.to(configs.device))
return energy
# ------------------------------------------------------------------
# Off-diagonal connections
# ------------------------------------------------------------------
[docs]
def get_connections(
self, config: torch.Tensor
) -> tuple[torch.Tensor, torch.Tensor]:
"""Return off-diagonal connected states (XX + YY exchange, X/Y fields).
The XX + YY interaction flips a pair of anti-aligned neighbours:
``|↑↓⟩ ↔ |↓↑⟩`` with amplitude ``0.5 * (Jx + Jy) / 2``.
When ``Jx ≠ Jy`` the ``|↑↓⟩ → |↓↑⟩`` and ``|↓↑⟩ → |↑↓⟩``
amplitudes differ.
The transverse fields h_x and h_y produce single-site flips.
Parameters
----------
config : torch.Tensor
Shape ``(num_sites,)``.
Returns
-------
connected_configs : torch.Tensor
Shape ``(n_conn, num_sites)``.
matrix_elements : torch.Tensor
Shape ``(n_conn,)``.
"""
new_configs: list[torch.Tensor] = []
elements: list[float] = []
# --- XX + YY pair exchange ---
# S^+_i S^-_j + S^-_i S^+_j contributes when spins i, j differ.
# S^x S^x + S^y S^y = 0.5*(S^+ S^- + S^- S^+) ... but with Jx, Jy:
# Jx S^x_i S^x_j + Jy S^y_i S^y_j
# = 0.25*(Jx - Jy)(S^+_i S^+_j + S^-_i S^-_j)
# + 0.25*(Jx + Jy)(S^+_i S^-_j + S^-_i S^+_j)
#
# S^+_i S^-_j flips |↓_i ↑_j⟩ → |↑_i ↓_j⟩ with amplitude 1
# (config[i]=1, config[j]=0) → (config[i]=0, config[j]=1)
# Similarly S^-_i S^+_j flips opposite way.
# S^+_i S^+_j flips |↓_i ↓_j⟩ → |↑_i ↑_j⟩ — non-Sz-conserving
# S^-_i S^-_j flips |↑_i ↑_j⟩ → |↓_i ↓_j⟩ — non-Sz-conserving
jx = self.Jx
jy = self.Jy
coeff_flip = 0.25 * (jx + jy) # S^+S^- + S^-S^+ part
coeff_double = 0.25 * (jx - jy) # S^+S^+ + S^-S^- part
for i, j in self._neighbours:
ci = int(config[i].item())
cj = int(config[j].item())
# S^+_i S^-_j: requires config[i]=1 (down), config[j]=0 (up)
if ci == 1 and cj == 0 and abs(coeff_flip) > 1e-15:
new_cfg = config.clone()
new_cfg[i] = 0
new_cfg[j] = 1
new_configs.append(new_cfg)
elements.append(coeff_flip)
# S^-_i S^+_j: requires config[i]=0 (up), config[j]=1 (down)
if ci == 0 and cj == 1 and abs(coeff_flip) > 1e-15:
new_cfg = config.clone()
new_cfg[i] = 1
new_cfg[j] = 0
new_configs.append(new_cfg)
elements.append(coeff_flip)
# S^+_i S^+_j: requires config[i]=1, config[j]=1
if ci == 1 and cj == 1 and abs(coeff_double) > 1e-15:
new_cfg = config.clone()
new_cfg[i] = 0
new_cfg[j] = 0
new_configs.append(new_cfg)
elements.append(coeff_double)
# S^-_i S^-_j: requires config[i]=0, config[j]=0
if ci == 0 and cj == 0 and abs(coeff_double) > 1e-15:
new_cfg = config.clone()
new_cfg[i] = 1
new_cfg[j] = 1
new_configs.append(new_cfg)
elements.append(coeff_double)
# --- Single-site X field ---
# h_x S^x = 0.5 * h_x * (S^+ + S^-) → single flip
for site in range(self.num_sites):
hx_val = float(self._h_x[site].item())
if abs(hx_val) < 1e-15:
continue
new_cfg = config.clone()
new_cfg[site] = 1 - int(config[site].item())
new_configs.append(new_cfg)
elements.append(0.5 * hx_val)
# --- Single-site Y field ---
# h_y S^y = 0.5 * h_y * (S^+ - S^-) / i
# S^+ |↓⟩ = |↑⟩ → coefficient = h_y * 0.5 * (-i) * 1 = ...
# For real Hamiltonians we keep the imaginary part explicit.
# S^y = (S^+ - S^-)/(2i)
# h_y S^y |↓⟩ = h_y/(2i) |↑⟩ (config 1 → 0, coeff = -i*h_y/2)
# h_y S^y |↑⟩ = -h_y/(2i) |↓⟩ (config 0 → 1, coeff = i*h_y/2)
# Since we work with real Hamiltonians, h_y typically = 0.
# We include it for completeness but note the matrix becomes complex.
for site in range(self.num_sites):
hy_val = float(self._h_y[site].item())
if abs(hy_val) < 1e-15:
continue
ci = int(config[site].item())
new_cfg = config.clone()
new_cfg[site] = 1 - ci
# S^y |↑⟩ = (i/2)|↓⟩, S^y |↓⟩ = (-i/2)|↑⟩
# For real representation we store the imaginary coefficient
# as a real number (valid when H is purely real or the caller
# handles complex arithmetic externally).
# Convention: store ±h_y/2 and let caller interpret.
if ci == 0:
# S^y |↑⟩ = (i/2)|↓⟩ → purely imaginary, skip for real H
# For a generic implementation, one would need complex tensors.
# We log a warning and skip.
pass
else:
pass
# NOTE: Y-field contributions are purely imaginary in the
# computational basis and are omitted from this real-valued
# implementation. Set h_y=0 (the default) for standard usage.
if not new_configs:
return (
torch.empty(
(0, self.num_sites), dtype=torch.int64, device=config.device
),
torch.empty(0, dtype=torch.float64, device=config.device),
)
return (
torch.stack(new_configs).to(dtype=torch.int64, device=config.device),
torch.tensor(elements, dtype=torch.float64, device=config.device),
)