"""Plaquette interaction terms and plaquette observables on lattice models.
This module provides :class:`PlaquetteTerm`, a ``QMBTerm``
subclass for constructing four-body plaquette Hamiltonian contributions and for
measuring plaquette expectation values on quantum many-body states.
"""
import logging
from itertools import product
from math import prod
import numpy as np
from scipy.sparse import csr_matrix, isspmatrix
from edlgt.symmetries import nbody_term
from edlgt.tools import validate_parameters
from .lattice_geometry import get_plaquette_neighbors
from .lattice_mappings import zig_zag
from .qmb_operations import four_body_op
from .qmb_state import QMB_state
from .qmb_term import QMBTerm
logger = logging.getLogger(__name__)
__all__ = ["PlaquetteTerm"]
[docs]
class PlaquetteTerm(QMBTerm):
"""Four-body plaquette term on a lattice.
The class supports both:
- direct sparse-matrix construction (no symmetry-sector reduction), and
- symmetry-sector workflows where Hamiltonian contributions are returned as
``(rows, cols, vals)`` triplets.
"""
def __init__(self, axes, op_list, op_names_list, print_plaq=True, **kwargs):
"""Initialize a plaquette term definition.
Parameters
----------
axes : list
Two lattice axes defining the plaquette plane (for example
``["x", "y"]``).
op_list : list
Operators used to build the plaquette term.
op_names_list : list
Human-readable names corresponding to ``op_list``.
print_plaq : bool, optional
If ``True``, print plaquette values when calling
:meth:`get_expval`.
**kwargs
Additional keyword arguments forwarded to :class:`QMBTerm`.
Returns
-------
None
"""
validate_parameters(
axes=axes,
op_list=op_list,
op_names_list=op_names_list,
print_plaq=print_plaq,
)
# Preprocess arguments
super().__init__(op_list=op_list, op_names_list=op_names_list, **kwargs)
self.axes = axes
self.print_plaq = print_plaq
# Number of lattice sites
self.n_sites = prod(self.lvals)
self.obs = None
self.var = None
self.avg = None
self.std = None
axes_name = "".join(self.axes)
logger.info("PlaqTerm %s: %s", axes_name, op_names_list)
[docs]
def get_hamiltonian(self, strength, add_dagger=False, mask=None):
"""Build the plaquette Hamiltonian contribution.
The plaquette term is summed over all lattice sites where a plaquette is
defined (and where ``mask`` allows it), then multiplied by ``strength``.
Parameters
----------
strength : scalar
Coupling constant multiplying the plaquette term.
add_dagger : bool, optional
If ``True``, add the Hermitian conjugate of the constructed term.
mask : numpy.ndarray, optional
Boolean mask controlling where the local term is applied.
Returns
-------
scipy.sparse.csr_matrix or tuple
Return type depends on the current workflow:
- if ``self.sector_configs is None``: sparse matrix Hamiltonian term;
- otherwise: ``(row_list, col_list, val_list)`` as three NumPy arrays in the
symmetry-reduced basis.
Raises
------
TypeError
If ``strength`` is not scalar or ``add_dagger`` is invalid.
"""
if not np.isscalar(strength):
raise TypeError(f"strength must be scalar, not {type(strength)}")
validate_parameters(add_dagger=add_dagger)
# --------------------------------------------------------------------
# 1) No sector_configs ⇒ build one big sparse matrix
if self.sector_configs is None:
h_plaquette = 0
for ii in range(self.n_sites):
coords = zig_zag(self.lvals, ii)
_, sites = get_plaquette_neighbors(
coords, self.lvals, self.axes, self.has_obc
)
if sites is None or not self.get_mask_conditions(coords, mask):
continue
h_plaquette += strength * four_body_op(
op_list=self.op_list, op_sites_list=sites, **self.def_params
)
if not isspmatrix(h_plaquette):
h_plaquette = csr_matrix(h_plaquette)
if add_dagger:
h_plaquette = h_plaquette + csr_matrix(h_plaquette.conj().T)
return h_plaquette
# --------------------------------------------------------------------
# 2) With sector_configs ⇒ collect (r,c,v) arrays
all_r = []
all_c = []
all_v = []
for ii in range(self.n_sites):
coords = zig_zag(self.lvals, ii)
_, sites = get_plaquette_neighbors(
coords, self.lvals, self.axes, self.has_obc
)
# logger.info(f"coords: {ii} {coords}, sites: {sites}")
if sites is None or not self.get_mask_conditions(coords, mask):
continue
if len(self.sector_configs) > 2**24.5:
logger.info("Sites %s", sites)
# this gives three 1D arrays for this plaquette
sites_array = np.array(sites, dtype=np.int32)
row_triplet, col_triplet, value_triplet = nbody_term(
op_list=self.sym_ops,
op_sites_list=sites_array,
sector_configs=self.sector_configs,
momentum_basis=self.momentum_basis,
)
all_r.append(row_triplet)
all_c.append(col_triplet)
all_v.append(value_triplet)
# merge them
row_list = np.concatenate(all_r)
col_list = np.concatenate(all_c)
val_list = np.concatenate(all_v) * strength
if add_dagger:
# Add Hermitian conjugate after the full term construction
dagger_row_list = col_list
dagger_col_list = row_list
dagger_val_list = np.conjugate(val_list)
# Concatenate the dagger part to the original term
row_list = np.concatenate([row_list, dagger_row_list])
col_list = np.concatenate([col_list, dagger_col_list])
val_list = np.concatenate([val_list, dagger_val_list])
return row_list, col_list, val_list
get_Hamiltonian = get_hamiltonian # pylint: disable=invalid-name
[docs]
def get_expval(self, psi, component: str = "real", stag_label: str | None = None):
"""Compute plaquette expectation values site by site and aggregate statistics.
Parameters
----------
psi : QMB_state
Quantum many-body state on which the expectation values are
evaluated.
component : str, optional
Component of the plaquette operator to measure. Allowed values are
``"real"`` (Hermitian part) and ``"imag"`` (anti-Hermitian part).
stag_label : str, optional
Optional staggered-site selector passed through the common validation
logic. Allowed values are ``"even"`` and ``"odd"``. If ``None``,
all plaquettes are considered.
Returns
-------
None
Results are stored on the instance attributes ``obs``, ``var``,
``avg``, and ``std``.
Raises
------
TypeError
If ``psi`` is not a ``QMB_state`` instance.
ValueError
If ``component`` is not ``"real"`` or ``"imag"``.
Notes
-----
In symmetry-sector mode, variances are currently not computed.
"""
# Check on parameters
if not isinstance(psi, QMB_state):
raise TypeError(f"psi must be instance of class:QMB_state not {type(psi)}")
validate_parameters(stag_label=stag_label)
if component not in ["real", "imag"]:
raise ValueError(f"component must be 'real' or 'imag': got {component}")
# ADVERTISE OF THE CHOSEN PART OF THE PLAQUETTE YOU WANT TO COMPUTE
if self.print_plaq:
logger.info("----------------------------------------------------")
logger.info("PLAQUETTE: %s", " ".join(self.op_names_list))
logger.info("----------------------------------------------------")
plaquettes = list(self.iter_plaquettes())
self.obs = self.get_empty_obs_array()
self.var = self.get_empty_obs_array()
# IN CASE OF NO SYMMETRY SECTOR
if self.sector_configs is None:
for origin_coords, coords_list, sites_list in plaquettes:
plaq_op = four_body_op(
op_list=self.op_list,
op_sites_list=sites_list,
**self.def_params,
)
if component == "real":
obs_op = 0.5 * (plaq_op + plaq_op.getH())
elif component == "imag":
obs_op = (-0.5j) * (plaq_op - plaq_op.getH())
else:
raise ValueError(
f"component must be 'real' or 'imag', not {component}"
)
self.obs[origin_coords] = np.real(np.vdot(psi.psi, obs_op.dot(psi.psi)))
tmp = obs_op.dot(psi.psi)
self.var[origin_coords] = np.real(np.vdot(tmp, tmp))
self.var[origin_coords] -= self.obs[origin_coords] ** 2
if self.print_plaq:
plaq_strings = [f"{coord_label}" for coord_label in coords_list]
self.print_plaquette(plaq_strings, self.obs[origin_coords])
else:
# GET THE EXPVAL ON THE SYMMETRY SECTOR
for origin_coords, coords_list, sites_list in plaquettes:
rows, cols, vals = nbody_term(
op_list=self.sym_ops,
op_sites_list=np.array(sites_list),
sector_configs=self.sector_configs,
momentum_basis=self.momentum_basis,
)
self.obs[origin_coords] = psi.expectation_value(
(rows, cols, vals), component
)
# for the moment, variance is not computed in symmetry sectors
if self.print_plaq:
plaq_strings = [f"{coord_label}" for coord_label in coords_list]
self.print_plaquette(plaq_strings, self.obs[origin_coords])
# GET STATISTICS
self.avg = np.mean(self.obs)
self.std = np.sqrt(np.mean(np.abs(self.var)))
if self.print_plaq:
logger.info("%.10f +/- %.10f", self.avg, self.std)
[docs]
def print_plaquette(self, sites_list, value):
"""Log a formatted ASCII representation of a plaquette value.
Parameters
----------
sites_list : list
List of four coordinate labels describing the plaquette corners.
value : float
Plaquette value to print.
Returns
-------
None
Raises
------
TypeError
If ``sites_list`` is not a list or ``value`` is not a float.
ValueError
If ``sites_list`` does not contain exactly four entries.
"""
if not isinstance(sites_list, list):
raise TypeError(f"sites_list should be a LIST, not a {type(sites_list)}")
if len(sites_list) != 4:
raise ValueError(f"sites_list has 4 elements, not {str(len(sites_list))}")
if not isinstance(value, float):
raise TypeError(f"sites_list should be FLOAT, not a {type(value)}")
if value > 0:
value = format(value, ".10f")
else:
if np.abs(value) < 10 ** (-10):
value = format(np.abs(value), ".10f")
else:
value = format(value, ".9f")
logger.info("%s------------%s", sites_list[2], sites_list[3])
logger.info(" | |")
logger.info(" | %s |", value)
logger.info(" | |")
logger.info("%s------------%s", sites_list[0], sites_list[1])
logger.info("")
print_Plaquette = print_plaquette # pylint: disable=invalid-name
[docs]
def get_plaquette_shape(self):
"""Return the resolved plaquette-array shape for this term.
Plane axes have one fewer plaquette origin under open boundary
conditions and one origin per lattice site under periodic boundary
conditions. Spectator axes keep the full lattice extent.
"""
dims = "xyz"[: len(self.lvals)]
axes_idx = {dims.index(axis) for axis in self.axes}
shape = []
for idx, length in enumerate(self.lvals):
if idx in axes_idx and self.has_obc[idx]:
shape.append(max(length - 1, 0))
else:
shape.append(length)
return tuple(shape)
[docs]
def iter_plaquettes(self):
"""Enumerate valid plaquette origins, corner coordinates, and sites."""
shape_ranges = (range(length) for length in self.get_plaquette_shape())
for origin_coords in product(*shape_ranges):
coords_list, sites_list = get_plaquette_neighbors(
origin_coords, self.lvals, self.axes, self.has_obc
)
if sites_list is not None:
yield origin_coords, coords_list, sites_list
[docs]
def get_empty_obs_array(self, fill_value=0.0):
"""Return a plaquette-resolved observable array with the proper shape."""
return np.full(self.get_plaquette_shape(), fill_value, dtype=float)