"""SU(2) lattice gauge-theory model helpers."""
import logging
import numpy as np
from numba import typed
from scipy.sparse import csr_matrix, identity
from scipy.sparse.linalg import norm as spnorm
from sympy import Rational, S
from edlgt.modeling import (
LocalTerm,
NBodyTerm,
PlaquetteTerm,
QMB_hamiltonian,
QMB_state,
TwoBodyTerm,
check_link_symmetry,
get_origin_surfaces,
staggered_mask,
)
from edlgt.operators import (
SU2_dressed_site_operators,
SU2_gauge_integrated_operators,
SU2_gauge_invariant_states,
SU2_gen_dressed_site_operators,
)
from edlgt.symmetries import build_parity_operator
from .quantum_model import (
QuantumModel,
format_loc_dims,
validate_staggered_matter_lattice,
)
logger = logging.getLogger(__name__)
__all__ = ["SU2_Model"]
[docs]
class SU2_Model(QuantumModel):
"""SU(2) lattice gauge model with hardcoded and generalized operator sets."""
def __init__(
self,
spin,
pure_theory,
bg_list=None,
sectors=None,
use_generic_model=False,
n_flavors=None,
**kwargs,
):
"""Initialize the SU(2) model and construct its symmetry sector.
Parameters
----------
spin : float or str = integrated
Gauge-link spin representation.
pure_theory : bool
If ``True``, exclude matter fields.
bg_list : list, optional
Optional background-charge specification.
sectors : list, optional
Global matter-sector labels. In the current SU(2) matter models
this selects the total particle-number sector. For the integrated
1D model, the exact global ``Delta_Z_tot = 0`` constraint is added
automatically on top of that particle sector.
use_generic_model : bool, optional
If ``True``, force the generalized operator construction.
n_flavors : int, optional
Number of matter flavors. Defaults to ``1`` for matter theories.
Multi-flavor Hamiltonians are not yet implemented; passing
``n_flavors > 1`` raises :class:`NotImplementedError`. The
flavored gauge-invariant local basis is still available via
:func:`~edlgt.operators.SU2_gauge_invariant_states`.
**kwargs
Arguments forwarded to :class:`~edlgt.models.quantum_model.QuantumModel`.
"""
# Initialize base class with the common parameters
super().__init__(**kwargs)
n_flavors_eff = 0 if pure_theory else (1 if n_flavors is None else int(n_flavors))
if n_flavors_eff > 1:
raise NotImplementedError(
"SU2_Model does not yet support n_flavors > 1. The flavored "
"gauge-invariant local basis is available via "
"SU2_gauge_invariant_states(..., n_flavors=N); the matching "
"fermionic operators and Hamiltonian terms are tracked in "
"docs/su2_singlets_optimization/follow_ups.md (F2)."
)
if not pure_theory:
validate_staggered_matter_lattice(self.lvals, "SU2", self.has_obc)
# SU(2) dressed-site operators are intrinsically complex.
self.configure_dtype_mode(dtype_mode="complex", auto_mode="complex")
self.spin = spin
self.pure_theory = pure_theory
self.use_generic_model = use_generic_model
if isinstance(spin, str):
if spin != "integrated":
raise ValueError(f"spin must be numeric or 'integrated': got {spin}")
if pure_theory or self.dim != 1:
msg = "Integrated SU2 is defined only in 1D with dynamical matter."
raise ValueError(msg)
if bg_list is not None:
raise ValueError("Integrated SU2 currently does not support bg_list")
self.background = 0
self.bg_list = None
self.gauge_basis = None
self.singlet_states = None
self.gauge_states = None
logger.info("----------------------------------------------------")
logger.info("(%s+1)D SU(2) LGT with gauge fields integrated", self.dim)
logger.info("----------------------------------------------------")
ops = SU2_gauge_integrated_operators()
self.project_operators(ops)
elif np.isscalar(spin):
if bg_list is None:
self.background = S(0)
self.bg_list = None
else:
S_bg_list = [Rational(str(bg)) for bg in bg_list]
self.background = max((bg for bg in S_bg_list), default=S(0))
self.bg_list = S_bg_list if self.background != 0 else None
pure_label = "pure" if self.pure_theory else "with matter"
logger.info("----------------------------------------------------")
msg = f"({self.dim}+1)D SU(2) LGT {pure_label} j={spin}"
logger.info(msg)
logger.info("----------------------------------------------------")
# ---------------------------------------------------------------------------
# Acquire gauge invariant basis and states
self.gauge_basis, self.singlet_states = SU2_gauge_invariant_states(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
)
self.gauge_states = {}
for key in self.singlet_states:
self.gauge_states[key] = np.array(
[singlet.J_config for singlet in self.singlet_states[key]],
dtype=object,
)
# ---------------------------------------------------------------------------
# Acquire operators
if self.spin < 1 and not self.use_generic_model:
ops = SU2_dressed_site_operators(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
)
else:
ops = SU2_gen_dressed_site_operators(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
)
# Initialize the operators, local dimension and lattice labels
self.project_operators(ops, bg_sector_list=self.bg_list)
if self.bg_list is not None:
logger.info("----------------------------------------------------")
logger.info("BACKGROUND CHARGES per SITE")
format_loc_dims(self.bg_list, self.lvals)
else:
raise ValueError(f"spin must be numeric or 'integrated': got {spin}")
# -------------------------------------------------------------------------------
# GLOBAL SYMMETRIES
if self.pure_theory:
global_ops = None
global_sectors = None
self.zero_density = None
else:
sector_values = np.atleast_1d(
[self.n_sites] if sectors is None else sectors
).astype(float)
particle_sector = float(sector_values[0])
global_ops = [self.ops["N_tot"]]
global_sectors = [particle_sector]
# Add an extra global constraint in the integrated model
if self.spin == "integrated":
global_ops.append(self.ops["Delta_Z"])
global_sectors.append(0.0)
elif sector_values.size != 1:
raise ValueError("SU2 matter sectors expect a single particle sector.")
self.zero_density = bool(np.isclose(particle_sector, self.n_sites))
# -------------------------------------------------------------------------------
# LINK SYMMETRIES
if self.spin == "integrated":
link_ops = None
link_sectors = None
else:
dirs = self.directions
link_ops = [[self.ops[f"T2_p{d}"], -self.ops[f"T2_m{d}"]] for d in dirs]
link_sectors = [0 for _ in dirs]
# -------------------------------------------------------------------------------
# SU2 ELECTRIC-FLUX “NBODY” SYMMETRIES
# only in the pure (no-matter) theory, more than 1D, *and* PBC
# Constrain, for each direction, the parity of the face/line through the origin:
# 3D: (Ex → yz-face at x=0) (Ey → xz-face at y=0) (Ez → xy-face at z=0)
# 2D: (Ex → y-axis at x=0) (Ey → x-axis at y=0)
if self.pure_theory and not any(self.has_obc):
logger.info("fixing surface parity fluxes")
# one flux‐constraint per cartesian direction
nbody_sectors = list(np.ones(self.dim, dtype=float))
nbody_ops = []
nbody_sites_list = typed.List()
surfaces = get_origin_surfaces(self.lvals)
nbody_sym_type = "Z"
if self.dim == 2:
# in 2D we have two lines through (0,0):
line_of = {"x": "y", "y": "x"}
for direction in self.directions:
sites = np.array(surfaces[line_of[direction]][1], dtype=np.uint8)
nbody_sites_list.append(sites)
nbody_ops.append(self.ops[f"P_p{direction}"])
elif self.dim == 3:
# in 3D we have three faces through (0,0,0):
face_of = {"x": "yz", "y": "xz", "z": "xy"}
for direction in self.directions:
sites = np.array(surfaces[face_of[direction]][1], dtype=np.uint8)
nbody_sites_list.append(sites)
logger.debug(
"%s sites: %s %s",
direction,
sites,
surfaces[face_of[direction]][0],
)
nbody_ops.append(self.ops[f"P_p{direction}"])
else:
# no electric‐flux constraint in 1D, or in OBC or with matter
nbody_sectors = None
nbody_ops = None
nbody_sites_list = None
nbody_sym_type = None
# -------------------------------------------------------------------------------
# GET SYMMETRY SECTOR
self.get_abelian_symmetry_sector(
global_ops=global_ops,
global_sectors=global_sectors,
link_ops=link_ops,
link_sectors=link_sectors,
nbody_ops=nbody_ops,
nbody_sectors=nbody_sectors,
nbody_sites_list=nbody_sites_list,
nbody_sym_type=nbody_sym_type,
)
if self.sector_configs is None:
raise ValueError("No configurations found for the given symmetry sectors")
[docs]
def build_Hamiltonian(
self, g, m=None, theta=0.0, lambda_noise=0.0, dtype_mode="auto"
):
"""Dispatch to the hardcoded or generalized SU(2) Hamiltonian builder.
Parameters
----------
dtype_mode : str or bool, optional
``"auto"``, ``"real"``, ``"complex"``, or legacy bool flag.
"""
resolved_mode = self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=m, theta=theta),
)
if self.spin == "integrated":
return self.build_integrated_Hamiltonian(g, m, dtype_mode=resolved_mode)
else:
if self.spin < 1 and not self.use_generic_model:
self.build_base_Hamiltonian(
g, m, theta, lambda_noise, dtype_mode=resolved_mode
)
else:
self.build_gen_Hamiltonian(g, m, dtype_mode=resolved_mode)
def _get_auto_dtype_mode(self, m=None, theta=0.0):
"""Heuristic dtype mode for SU(2) Hamiltonian assembly."""
if self.pure_theory and np.abs(theta) <= 1e-12:
return "real"
return "complex"
[docs]
def build_base_Hamiltonian(
self, g, m=None, theta=0.0, lambda_noise=0.0, dtype_mode="auto"
):
"""Assemble the hardcoded low-spin SU(2) Hamiltonian.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare mass parameter.
theta : float, optional
Topological-angle parameter.
lambda_noise : float, optional
Optional Gauss-law-violating noise strength.
dtype_mode : str or bool, optional
``"auto"``, ``"real"``, ``"complex"``, or legacy bool flag.
"""
logger.info("----------------------------------------------------")
logger.info("BUILDING (H)ardocore (G)luon J=1/2 HAMILTONIAN")
logger.info("g=%s, m=%s, theta=%s, lambda_noise=%s", g, m, theta, lambda_noise)
logger.info("----------------------------------------------------")
self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=m, theta=theta),
)
# Hamiltonian Coefficients
self.SU2_Hamiltonian_couplings(g, m, theta)
h_terms = {}
# ---------------------------------------------------------------------------
# ELECTRIC ENERGY
op_name = "E_square"
h_terms[op_name] = LocalTerm(self.ops[op_name], op_name, **self.def_params)
self.H.add_term(h_terms[op_name].get_Hamiltonian(strength=self.coeffs["E"]))
# ---------------------------------------------------------------------------
# PLAQUETTE TERM: MAGNETIC INTERACTION
if self.dim > 1:
op_names_list = ["C_px,py", "C_py,mx", "C_my,px", "C_mx,my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_xy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_xy"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
if self.dim == 3:
# XZ Plane
op_names_list = ["C_px,pz", "C_pz,mx", "C_mz,px", "C_mx,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_xz"] = PlaquetteTerm(
["x", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_xz"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# YZ Plane
op_names_list = ["C_py,pz", "C_pz,my", "C_mz,py", "C_my,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["plaq_yz"] = PlaquetteTerm(
["y", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["plaq_yz"].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# ---------------------------------------------------------------------------
if not self.pure_theory:
# -----------------------------------------------------------------------
# STAGGERED MASS TERM
for site in ["even", "odd"]:
h_terms[f"N_{site}"] = LocalTerm(
self.ops["N-1"], "N-1", **self.def_params
)
self.H.add_term(
h_terms[f"N_{site}"].get_Hamiltonian(
self.coeffs[f"m_{site}"], staggered_mask(self.lvals, site)
)
)
# -----------------------------------------------------------------------
# HOPPING
for d in self.directions:
for site in ["even", "odd"]:
op_names_list = [f"Qp{d}_dag", f"Qm{d}"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms[f"{d}_hop_{site}"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
mask = staggered_mask(self.lvals, site)
self.H.add_term(
h_terms[f"{d}_hop_{site}"].get_Hamiltonian(
strength=self.coeffs[f"t{d}_{site}"],
add_dagger=True,
mask=mask,
)
)
# -------------------------------------------------------------------------------
# TOPOLOGICAL TERM
if self.dim == 3 and np.abs(self.coeffs["theta"]) > 10e-10:
logger.info("Adding topological term")
# XY Plane
op_names_list = ["EzC_px,py", "C_py,mx", "C_my,px", "C_mx,my"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ez_Bxy"] = PlaquetteTerm(
["x", "y"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ez_Bxy"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# XZ Plane
op_names_list = ["EyC_px,pz", "C_pz,mx", "C_mz,px", "C_mx,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ey_Bxz"] = PlaquetteTerm(
["x", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ey_Bxz"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# YZ Plane
op_names_list = ["ExC_py,pz", "C_pz,my", "C_mz,py", "C_my,mz"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["Ex_Byz"] = PlaquetteTerm(
["y", "z"], op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["Ex_Byz"].get_Hamiltonian(
strength=self.coeffs["theta"], add_dagger=True
)
)
# RANDOM NOISE VIOLATING GAUSS LAW
noise_requirements = [self.background > 0, lambda_noise > 0]
noise_requirements += [not obc for obc in self.has_obc]
if np.all(noise_requirements):
logger.info("SU2 GAUSS LAW VIOLATING TERM")
# HOPPING
op_names_list = ["Vpx_dag", "Vmx"]
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms["V_hop"] = TwoBodyTerm(
"x", op_list, op_names_list, **self.def_params
)
self.H.add_term(
h_terms["V_hop"].get_Hamiltonian(
strength=-lambda_noise, add_dagger=True
)
)
self.H.build(self.ham_format)
[docs]
def build_gen_Hamiltonian(self, g, m=None, dtype_mode="auto"):
"""Assemble the generalized SU(2) Hamiltonian.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare mass parameter.
dtype_mode : str or bool, optional
``"auto"``, ``"real"``, ``"complex"``, or legacy bool flag.
"""
logger.info("----------------------------------------------------")
logger.info("BUILDING generalized HAMILTONIAN with g=%s, m=%s", g, m)
self.configure_dtype_mode(
dtype_mode=dtype_mode,
auto_mode=self._get_auto_dtype_mode(m=m, theta=0.0),
)
# Hamiltonian Coefficients
self.SU2_Hamiltonian_couplings(g, m)
h_terms = {}
# ---------------------------------------------------------------------------
# ELECTRIC ENERGY
op_name = "E_square"
h_terms[op_name] = LocalTerm(self.ops[op_name], op_name, **self.def_params)
self.H.add_term(h_terms[op_name].get_Hamiltonian(strength=self.coeffs["E"]))
# ---------------------------------------------------------------------------
if not self.pure_theory:
# -----------------------------------------------------------------------
# STAGGERED MASS TERM
for site in ["even", "odd"]:
h_terms[f"N_{site}"] = LocalTerm(
self.ops["N_tot"], "N_tot", **self.def_params
)
self.H.add_term(
h_terms[f"N_{site}"].get_Hamiltonian(
self.coeffs[f"m_{site}"], staggered_mask(self.lvals, site)
)
)
# --------------------------------------------------------------------
# Generalized HOPPING
for d in self.directions:
for site in ["even", "odd"]:
hopping_terms = [
[f"Q1_p{d}_dag", f"Q2_m{d}"],
[f"Q2_p{d}_dag", f"Q1_m{d}"],
]
for ii, op_names_list in enumerate(hopping_terms):
op_list = [self.ops[op] for op in op_names_list]
# Define the Hamiltonian term
h_terms[f"{d}{ii}_hop_{site}"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
mask = staggered_mask(self.lvals, site)
self.H.add_term(
h_terms[f"{d}{ii}_hop_{site}"].get_Hamiltonian(
strength=self.coeffs[f"t{d}_{site}"],
add_dagger=True,
mask=mask,
)
)
# -------------------------------------------------------------------------------
# PLAQUETTE TERM: MAGNETIC INTERACTION
plaq_list = []
plaquette_directions = []
if self.dim == 3:
plaquette_directions = ["xy", "xz", "yz"]
elif self.dim == 2:
plaquette_directions = ["xy"]
plaquette_set = [
["AB", "AB", "AB", "AB"],
["AA", "AB", "BB", "AB"],
["AB", "AB", "AA", "BB"],
["AA", "AB", "BA", "BB"],
["AB", "BB", "AB", "AA"],
["AA", "BB", "BB", "AA"],
["AB", "BB", "AA", "BA"],
["AA", "BB", "BA", "BA"],
["BB", "AA", "AB", "AB"],
["BA", "AA", "BB", "AB"],
["BB", "AA", "AA", "BB"],
["BA", "AA", "BA", "BB"],
["BB", "BA", "AB", "AA"],
["BA", "BA", "BB", "AA"],
["BB", "BA", "AA", "BA"],
["BA", "BA", "BA", "BA"],
]
if self.dim > 1:
for pdir in plaquette_directions:
for p_set in plaquette_set:
# DEFINE THE LIST OF CORNER OPERATORS
op_names_list = [
f"C{p_set[0]}_p{pdir[0]},p{pdir[1]}",
f"C{p_set[1]}_p{pdir[1]},m{pdir[0]}",
f"C{p_set[2]}_m{pdir[1]},p{pdir[0]}",
f"C{p_set[3]}_m{pdir[0]},m{pdir[1]}",
]
# CORRESPONDING LIST OF OPERATORS
op_list = [self.ops[op] for op in op_names_list]
# DEFINE THE PLAQUETTE CLASS
plaq_name = f"P{pdir}_" + "".join(p_set)
h_terms[plaq_name] = PlaquetteTerm(
[pdir[0], pdir[1]],
op_list,
op_names_list,
print_plaq=False,
**self.def_params,
)
# ADD THE HAMILTONIAN TERM
self.H.add_term(
h_terms[plaq_name].get_Hamiltonian(
strength=self.coeffs["B"], add_dagger=True
)
)
# ADD THE PLAQUETTE TO THE LIST OF OBSERVABLES
plaq_list.append(plaq_name)
self.H.build(self.ham_format)
[docs]
def build_integrated_Hamiltonian(self, g, m, dtype_mode="auto"):
"""Assemble the integrated-gauge 1D SU2 Hamiltonian.
Parameters
----------
dtype_mode : str or bool, optional
``"auto"``, ``"real"``, ``"complex"``, or legacy bool flag.
"""
self.configure_dtype_mode(dtype_mode=dtype_mode, auto_mode="complex")
# Hamiltonian Coefficients
self.SU2_Hamiltonian_couplings(g, m)
logger.info("BUILDING integrated HAMILTONIAN")
if self.dim != 1 or self.pure_theory:
msg = "Integrated SU2 Hamiltonian is defined only in 1D with matter."
raise ValueError(msg)
h_terms = {}
n_links = self.n_sites - 1
# ---------------------------------------------------------------------------
# STAGGERED MASS TERM
h_terms["N"] = LocalTerm(self.ops["N_tot"], op_name="N_tot", **self.def_params)
for _, stag_label in enumerate(["even", "odd"]):
mask = staggered_mask(self.lvals, stag_label)
mass_coeff = self.coeffs[f"m_{stag_label}"]
self.H.add_term(h_terms["N"].get_Hamiltonian(mass_coeff, mask))
# ---------------------------------------------------------------------------
# HOPPING TERM
op_names_list = ["Q_r_dag", "Q_r"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["hop"] = TwoBodyTerm("x", op_list, op_names_list, **self.def_params)
self.H.add_term(
h_terms["hop"].get_Hamiltonian(strength=self.coeffs["t"], add_dagger=True)
)
op_names_list = ["Q_g_dag", "Q_g"]
op_list = [self.ops[op] for op in op_names_list]
h_terms["hop"] = TwoBodyTerm("x", op_list, op_names_list, **self.def_params)
self.H.add_term(
h_terms["hop"].get_Hamiltonian(strength=self.coeffs["t"], add_dagger=True)
)
# ---------------------------------------------------------------------------
# ELECTRIC ENERGY: long-range Hamiltonian
# Build and cache single-site masks once
site_masks = [self._single_site_mask(ii) for ii in range(self.n_sites)]
# TWOBODY TERM
twobody_terms = {
# Delta_Z_i Delta_Z_j part with J_ij = g^2/2 * (L-1-j)/8, i<j
"DeltaDelta": {
"op_list": [self.ops["Delta_Z"], self.ops["Delta_Z"]],
"op_names_list": ["Delta_Z", "Delta_Z"],
"prefactor": self.coeffs["E"] / 8,
},
# [SpSm, SmSp] part with J_ij = g^2/2 * (L-1-j), i<j
"SpSm_SmSp": {
"op_list": [self.ops["Sp_r_Sm_g"], self.ops["Sm_r_Sp_g"]],
"op_names_list": ["Sp_r_Sm_g", "Sm_r_Sp_g"],
"prefactor": self.coeffs["E"],
},
# [SmSp, SpSm] part with J_ij = g^2/2 * (L-1-j), i<j
"SmSp_SpSm": {
"op_list": [self.ops["Sm_r_Sp_g"], self.ops["Sp_r_Sm_g"]],
"op_names_list": ["Sm_r_Sp_g", "Sp_r_Sm_g"],
"prefactor": self.coeffs["E"],
},
}
for keyterm in twobody_terms.keys():
for dist in range(1, self.n_sites):
twobody_terms[keyterm][dist] = NBodyTerm(
twobody_terms[keyterm]["op_list"],
twobody_terms[keyterm]["op_names_list"],
distances=[(dist,)],
**self.def_params,
)
for kk in range(1, n_links):
# compute the coefficient (N-1-kk)
factor = self.n_sites - 1 - kk
weight = twobody_terms[keyterm]["prefactor"] * factor
if np.abs(weight) > 1e-14:
for ii in range(kk):
dist = kk - ii
self.H.add_term(
twobody_terms[keyterm][dist].get_Hamiltonian(
strength=weight, mask=site_masks[ii]
)
)
# ONEBODY TERM
h_terms["inv_Sz"] = LocalTerm(
self.ops["inv_Sz"], op_name="inv_Sz", **self.def_params
)
for kk in range(0, n_links):
mask = site_masks[kk]
factor = self.coeffs["E"] * (3 / 8) * (self.n_sites - 1 - kk)
self.H.add_term(
h_terms["inv_Sz"].get_Hamiltonian(strength=factor, mask=mask)
)
# ---------------------------------------------------------------------------
# TODO: implement static charges
self.H.build(format=self.ham_format)
[docs]
def reconstruct_Casimir_from_matter(
self,
single_obs_name: str = "N_single",
delta_corr_obs_name: str = "Delta_Z_Delta_Z",
flip_corr_obs_name: str = "Sp_r_Sm_g_Sm_r_Sp_g",
flop_corr_obs_name: str = "Sm_r_Sp_g_Sp_r_Sm_g",
state=None,
state_index: int | None = None,
dynamics: bool = False,
compute_single_obs: bool = True,
compute_pair_corr: bool = True,
print_values: bool = True,
) -> np.ndarray:
"""Reconstruct integrated 1D SU(2) link Casimirs from matter observables.
Notes
-----
For the current integrated 1D SU(2) model with OBC and no static
charges, the electric energy can be written in terms of the link
Casimirs
``T^2_n = (sum_{k=0}^n vec(T)_k)^2``.
In the local integrated basis, the single-site and pairwise pieces are
``T_k^2 = 3/4 * N_single(k)``
and
``2 vec(T)_i . vec(T)_j = Delta_Z(i) Delta_Z(j) / 8
+ Sp_r_Sm_g(i) Sm_r_Sp_g(j)
+ Sm_r_Sp_g(i) Sp_r_Sm_g(j)``.
"""
if self.spin != "integrated" or self.dim != 1 or self.pure_theory:
msg = "Integrated T2 reconstruction is defined only for 1D SU2 with matter."
raise ValueError(msg)
if not all(self.has_obc):
raise NotImplementedError(
"Integrated T2 reconstruction is currently implemented only with OBC."
)
if not hasattr(self, "def_params"):
self.default_params()
target_state = None
n_links = self.n_sites - 1
if n_links <= 0:
link_casimir = np.zeros(0, dtype=float)
self.res["T2"] = link_casimir
self.res["T2_avg"] = 0.0
self.res["E2"] = link_casimir
self.res["E2_avg"] = 0.0
self.res["T2_px"] = np.zeros(self.n_sites, dtype=float)
self.res["T2_mx"] = np.zeros(self.n_sites, dtype=float)
return link_casimir
reuse_cached_single = (
single_obs_name in self.res
and state is None
and state_index is None
and not dynamics
)
if reuse_cached_single:
single_values = np.asarray(self.res[single_obs_name], dtype=float)
elif compute_single_obs:
target_state = self._get_measurement_state(
state=state, state_index=state_index, dynamics=dynamics
)
single_term = LocalTerm(
self.ops["N_single"], op_name="N_single", **self.def_params
)
single_term.get_expval(target_state, print_values=False)
single_values = np.asarray(single_term.obs, dtype=float)
self.res[single_obs_name] = single_values
else:
raise KeyError(
f"Missing '{single_obs_name}' in self.res. Measure it first or "
"set compute_single_obs=True."
)
if single_values.shape != (self.n_sites,):
raise ValueError(
f"Expected '{single_obs_name}' shape ({self.n_sites},), got {single_values.shape}."
)
def get_pair_corr(obs_name: str, op_names_list: list[str]) -> np.ndarray:
reuse_cached_corr = (
obs_name in self.res
and not compute_pair_corr
and state is None
and state_index is None
and not dynamics
)
if reuse_cached_corr:
corr = np.asarray(self.res[obs_name], dtype=float)
else:
local_state = target_state
if local_state is None:
local_state = self._get_measurement_state(
state=state, state_index=state_index, dynamics=dynamics
)
corr_term = TwoBodyTerm(
axis="x",
op_list=[self.ops[name] for name in op_names_list],
op_names_list=op_names_list,
**self.def_params,
)
corr_term.get_expval(local_state)
corr = np.asarray(corr_term.corr, dtype=float)
self.res[obs_name] = corr
if corr.shape != (self.n_sites, self.n_sites):
raise ValueError(
f"Expected '{obs_name}' shape ({self.n_sites}, {self.n_sites}), got {corr.shape}."
)
return corr
delta_corr = get_pair_corr(delta_corr_obs_name, ["Delta_Z", "Delta_Z"])
flip_corr = get_pair_corr(flip_corr_obs_name, ["Sp_r_Sm_g", "Sm_r_Sp_g"])
flop_corr = get_pair_corr(flop_corr_obs_name, ["Sm_r_Sp_g", "Sp_r_Sm_g"])
prefix_single = np.cumsum(single_values, dtype=float)
pair_matrix = 0.125 * delta_corr + flip_corr + flop_corr
upper_pair = np.triu(pair_matrix, k=1)
pair_prefix = upper_pair.cumsum(axis=0).cumsum(axis=1)
link_casimir = 0.75 * prefix_single[:n_links] + np.diag(pair_prefix)[:n_links]
link_casimir = np.asarray(link_casimir, dtype=float)
t2_px = np.zeros(self.n_sites, dtype=float)
t2_mx = np.zeros(self.n_sites, dtype=float)
t2_px[:n_links] = link_casimir
t2_mx[1:] = link_casimir
self.res["T2"] = link_casimir
self.res["T2_avg"] = float(np.mean(link_casimir))
self.res["T2_px"] = t2_px
self.res["T2_mx"] = t2_mx
# Keep E2 aliases for workflow/result compatibility.
self.res["E2"] = link_casimir
self.res["E2_avg"] = self.res["T2_avg"]
if print_values:
logger.info("----------------------------------------------------")
logger.info("T2 ")
for ii in range(n_links):
logger.info("%s %.16f", (ii,), link_casimir[ii])
logger.info("%.16f", self.res["T2_avg"])
return link_casimir
[docs]
def get_fermionic_string_correlator(
self,
state=None,
state_index=None,
dynamics: bool = False,
print_values: bool = False,
):
"""Measure the gauge-invariant fermionic string correlator matrix.
Notes
-----
This implementation is currently defined for the 1D truncated SU(2)
dressed-site model, both in the hardcoded ``j=1/2`` construction and
in the generic truncated construction whenever the projected operator
set exposes the effective singlet operators ``Qpx_dag``, ``Qmx``, and
the scalar intermediate transporter ``W``.
A first attempt kept two explicit color labels at the endpoints and
produced a ``(2L, 2L)`` matrix. After projection to the gauge-invariant
dressed-site basis, however, the red and green endpoint channels become
identical. This is a direct consequence of working in the local
gauge-invariant dressed-site basis, where Gauss' law has already been
solved site by site. The physically meaningful covariance matrix
therefore has one effective fermionic mode per dressed site, i.e. size
``(L, L)``.
Pedagogically, the point is that the local dressed-site basis does not
keep bare color states such as ``|r>`` and ``|g>`` as independent
physical states. Instead, matter and rishons are first combined into
local SU(2)-invariant singlets. For a matter doublet and a rishon
doublet, the local color space decomposes as ``2 x 2 = 1 + 3``:
the projection keeps the singlet channel and discards the triplet. The
two endpoint color components are therefore just two representatives of
the same surviving singlet channel, and their difference is projected
out.
For ``i < j`` we measure
``C_ij = 1/4 <Qpx_dag(i) W(i+1) ... W(j-1) Qmx(j)>``
The factor ``1/4`` appears because, inside the projected basis, the
singlet endpoint operators satisfy ``Qpx_dag = 2 Fpx_eff_dag`` and
``Qmx = 2 Fmx_eff``. The intermediate color contraction is already
absorbed locally into the scalar dressed-site operator ``W``. On the
diagonal we store
``C_ii = <N_tot(i)> / 2``
because each effective site mode carries half of the total on-site
fermion number.
"""
if self.dim != 1 or self.pure_theory:
msg = "Fermionic string correlator is implemented only for 1D SU2 with matter."
raise NotImplementedError(msg)
if self.spin == "integrated":
msg = "Fermionic string correlator is implemented only in the truncated SU2"
raise NotImplementedError(msg)
required_ops = {"N_tot", "Qpx_dag", "Qmx", "W"}
if not required_ops.issubset(self.ops):
raise NotImplementedError(
"Fermionic string correlator currently requires 1D truncated "
"SU2 dressed-site operators exposing N_tot, Qpx_dag, Qmx, and W."
)
if not all(self.has_obc):
msg = "Fermionic string correlator is currently implemented only with OBC."
raise NotImplementedError(msg)
if not hasattr(self, "def_params"):
self.default_params()
target_state = self._get_measurement_state(
state=state, state_index=state_index, dynamics=dynamics
)
# The physical projected covariance matrix has one effective fermionic
# mode per dressed site, so it is LxL rather than 2Lx2L.
corr = np.zeros((self.n_sites, self.n_sites), dtype=np.complex128)
density_term = LocalTerm(self.ops["N_tot"], "N_tot", **self.def_params)
density_term.get_expval(target_state, print_values=False)
density_values = 0.5 * np.array(density_term.obs, dtype=float, copy=True)
for site_index in range(self.n_sites):
# The diagonal stores the effective site occupation C_ii=<N_tot>/2.
corr[site_index, site_index] = density_values[site_index]
for dist in range(1, self.n_sites):
distances = [(step,) for step in range(1, dist + 1)]
# In the projected dressed-site basis the physical endpoint channel
# is the surviving local singlet. In the hardcoded model these are
# the original Q operators, while in the generic truncated model
# they are reconstructed from the native Q1/Q2 half-link channels.
op_names_list = ["Qpx_dag"]
op_names_list += ["W" for _ in range(dist - 1)]
op_names_list += ["Qmx"]
op_list = [self.ops[name] for name in op_names_list]
# NBodyTerm evaluates the same ordered string on every interval
# (i, i+dist) in one pass.
corr_term = NBodyTerm(
op_list=op_list,
op_names_list=op_names_list,
distances=distances,
**self.def_params,
)
corr_term.get_expval(target_state, component="real", print_values=False)
real_values = np.array(corr_term.obs, dtype=float, copy=True)
corr_term.get_expval(target_state, component="imag", print_values=False)
imag_values = np.array(corr_term.obs, dtype=float, copy=True)
string_values = 0.25 * (real_values + 1j * imag_values)
for ii, value in enumerate(string_values):
jj = ii + dist
corr[ii, jj] = value
corr[jj, ii] = np.conjugate(value)
self.res["fermionic_string_correlator"] = corr
if print_values:
logger.info("----------------------------------------------------")
logger.info("FERMIONIC STRING CORRELATOR")
for ii in range(self.n_sites):
logger.info(
np.array2string(corr[ii], precision=3, suppress_small=False)
)
return corr
[docs]
def measure_fermionic_nongaussianity(
self,
state=None,
state_index=None,
dynamics: bool = False,
print_value: bool = True,
eig_tol: float = 1e-10,
):
"""Measure the fermionic non-Gaussianity from the string correlator."""
corr = self.get_fermionic_string_correlator(
state=state,
state_index=state_index,
dynamics=dynamics,
print_values=False,
)
herm_corr = 0.5 * (corr + corr.conj().T)
eigvals = np.linalg.eigvalsh(herm_corr)
eigvals = np.real_if_close(eigvals, tol=1000).astype(float)
min_eval = float(np.min(eigvals))
max_eval = float(np.max(eigvals))
if min_eval < -eig_tol or max_eval > 1 + eig_tol:
logger.warning(
"Fermionic string correlator has eigenvalues outside [0, 1]: "
"min=%.6e, max=%.6e. Clipping for NG.",
min_eval,
max_eval,
)
clipped = np.clip(eigvals, 0.0, 1.0)
entropy_vals = np.zeros_like(clipped)
valid = np.logical_and(clipped > eig_tol, clipped < 1.0 - eig_tol)
vv = clipped[valid]
entropy_vals[valid] = -(vv * np.log(vv) + (1.0 - vv) * np.log(1.0 - vv))
ng_value = float(np.sum(entropy_vals))
self.res["fermionic_string_correlator_eigvals"] = eigvals
self.res["fermionic_nongaussianity"] = ng_value
if print_value:
logger.info("----------------------------------------------------")
logger.info("FERMIONIC NON-GAUSSIANITY")
logger.info("eigvals")
logger.info(np.array2string(eigvals, precision=8, suppress_small=False))
logger.info("NG = %.16f", ng_value)
return ng_value, eigvals, corr
[docs]
def check_symmetries(self):
"""Check link-symmetry constraints on measured observables."""
# CHECK LINK SYMMETRIES
for ax in self.directions:
check_link_symmetry(
ax,
self.obs_list[f"T2_p{ax}"],
self.obs_list[f"T2_m{ax}"],
value=0,
sign=-1,
)
[docs]
def overlap_QMB_state(self, name):
"""Return predefined benchmark SU(2) basis configurations.
Parameters
----------
name : str
Label of a reference configuration.
Returns
-------
numpy.ndarray
Configuration in the model basis.
"""
# POLARIZED AND BARE VACUUM in 1D
if len(self.lvals) == 1:
L, R = 0, 0
if self.spin == "integrated":
if name == "V":
s1, s2, L, R = 0, 3, 0, 3
elif name == "PV":
s1, s2, L, R = 2, 1, 2, 1
else:
s1, s2, L, R = 0, 0, 0, 0
else:
if name == "V":
if self.spin < 1:
s1, s2, L, R = 0, 4, 0, 2
elif self.spin == 1:
s1, s2, L, R = 0, 7, 0, 2
elif self.spin > 1:
s1, s2, L, R = 0, 10, 0, 2
elif name == "PV":
if self.spin < 1:
s1, s2, L, R = 1, 5, 1, 1
elif self.spin == 1:
s1, s2, L, R = 1, 8, 1, 1
elif self.spin > 1:
s1, s2, L, R = 1, 11, 1, 1
elif name == "T":
s1, s2, L, R = 6, 12, 1, 1
elif name == "M":
s1, s2, L, R = 2, 3, 1, 1
elif name == "B":
s1, s2, L, R = 7, 11, 1, 1
else:
s1, s2, L, R = 0, 0, 0, 0
config_state = [s1 if ii % 2 == 0 else s2 for ii in range(self.n_sites)]
if name == "DW":
if self.spin == "integrated":
s1, s2, L, R = 0, 3, 0, 3
else:
s1, s2, L, R = 0, 4, 0, 2
config_state = [
s1 if ii < self.n_sites // 2 else s2 for ii in range(self.n_sites)
]
elif name == "DW2":
if self.spin == "integrated":
s1, s2, L, R = 0, 3, 0, 3
else:
s1, s2, L, R = 0, 4, 0, 2
config_state = [
s1 if (ii // 2) % 2 == 0 else s2 for ii in range(self.n_sites)
]
if self.has_obc[0]:
config_state[0] = L
config_state[-1] = R
# POLARIZED AND BARE VACUUM in 2D
else:
if not self.pure_theory:
if self.has_obc[0]:
if name == "V":
config_state = [0, 9, 0, 4, 4, 0, 9, 0]
elif name == "PV1":
config_state = [1, 12, 3, 5, 5, 2, 11, 1]
elif name == "PV2":
config_state = [1, 11, 1, 5, 5, 3, 10, 1]
else:
config_state = [0, 9, 0, 9, 9, 0, 9, 0]
else:
if name == "V":
config_state = [0 for _ in range(self.n_sites)]
config_state = config_state
return np.array(config_state)
[docs]
def SU2_Hamiltonian_couplings(self, g, m=None, theta=0.0):
"""Set SU(2) Hamiltonian couplings from physical parameters.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare mass parameter (used when matter fields are present).
theta : float, optional
Topological-angle parameter.
Returns
-------
None
Couplings are stored in ``self.coeffs``.
Notes
-----
In the current convention, the Hamiltonian is rescaled so that the
hopping term is dimensionless (as used in the project implementation and
in PRX Quantum 5, 040309).
The rescaling summary used in the code is:
- hopping: original coupling ``1/2`` -> ``2 * sqrt(2)``
- electric term: original ``g0^2 / 2`` -> ``8 g^2 / 3``
- magnetic term: rescaled convention factor applied in the model
- the symbol ``g`` in this implementation is used as the rescaled
coupling (effectively a ``g^2`` convention)
DFL-project convention (reference values often used in scripts):
- ``E = 8 * g / 3``
- ``B = -3 / g``
- ``t = 2 * sqrt(2)``
String-breaking convention (alternative reference choice):
- ``E = g``
- ``B = -1``
- ``t = 1``
"""
if self.dim == 1:
E = 8 * g / 3
B = 0
else:
E = g / 2
B = -1 / (2 * g)
# Dictionary with Hamiltonian COEFFICIENTS
self.coeffs = {
"g": g,
"E": E, # ELECTRIC FIELD coupling
"B": B, # MAGNETIC FIELD coupling
"theta": -theta * g, # THETA TERM coupling
}
if not self.pure_theory and m is not None:
# The correct hopping in original units should be 1/2
t = 2 * np.sqrt(2)
self.coeffs |= {
"t": -complex(0, t), # integrated-gauge hopping coefficient
"tx_even": -complex(0, t), # x HOPPING (EVEN SITES)
"tx_odd": -complex(0, t), # x HOPPING (ODD SITES)
"ty_even": -t, # y HOPPING (EVEN SITES)
"ty_odd": t, # y HOPPING (ODD SITES)
"tz_even": -t, # z HOPPING (EVEN SITES)
"tz_odd": t, # z HOPPING (ODD SITES)
"m": m,
"m_odd": -m, # EFFECTIVE MASS for ODD SITES
"m_even": m, # EFFECTIVE MASS for EVEN SITES
}
[docs]
def build_local_Hamiltonian(self, g, m, R0):
"""Build a local effective Hamiltonian around one 1D site.
Parameters
----------
g : float
Gauge coupling.
m : float
Bare mass parameter.
R0 : int
Central lattice site.
"""
logger.info("----------------------------------------------------")
logger.info("BUILDING local HAMILTONIAN around %s", R0)
# -------------------------------------------------------------------------------
if self.dim > 1:
raise ValueError(f"Local Hamiltonian valid only for D=1, got {self.dim}")
self.Hlocal = QMB_hamiltonian(self.lvals, size=self.sector_dim)
# -------------------------------------------------------------------------------
# Hamiltonian Coefficients
self.SU2_Hamiltonian_couplings(g, m)
# -------------------------------------------------------------------------------
hterms = {}
# ELECTRIC HAMILTONIAN
hterms["E2"] = LocalTerm(self.ops["E_square"], "E_square", **self.def_params)
# MASS TERM
hterms["N"] = LocalTerm(self.ops["N-1"], "N-1", **self.def_params)
# HOPPING
op_names_list = ["Qpx_dag", "Qmx"]
op_list = [self.ops[op] for op in op_names_list]
hterms["hop"] = TwoBodyTerm("x", op_list, op_names_list, **self.def_params)
# DAGGER HOPPING
op_names_list = ["Qpx", "Qmx_dag"]
op_list = [self.ops[op] for op in op_names_list]
hterms["hop_dag"] = TwoBodyTerm("x", op_list, op_names_list, **self.def_params)
hop_coeff = self.coeffs["tx_even"]
hop_coeff_half = 0.5 * hop_coeff
# ---------------------------------------------------------------------------
self.Hlocal.add_term(
hterms["E2"].get_Hamiltonian(self.coeffs["E"], self.get_mask([R0]))
)
# ---------------------------------------------------------------------------
self.Hlocal.add_term(
hterms["N"].get_Hamiltonian(
((-1) ** R0) * self.coeffs["m"], self.get_mask([R0])
)
)
# ---------------------------------------------------------------------------
# Add the hopping term (j,j+1)
self.Hlocal.add_term(
hterms["hop"].get_Hamiltonian(hop_coeff_half, mask=self.get_mask([R0]))
)
# ---------------------------------------------------------------------------
# Add the term (j-1,j)
self.Hlocal.add_term(
hterms["hop"].get_Hamiltonian(
hop_coeff_half,
mask=self.get_mask([(R0 - 1) % self.n_sites]),
)
)
# ---------------------------------------------------------------------------
# Add the hermitian conjugate
# Add the term (j,j+1)
self.Hlocal.add_term(
hterms["hop_dag"].get_Hamiltonian(-hop_coeff_half, mask=self.get_mask([R0]))
)
# ---------------------------------------------------------------------------
# Add the term (j-1,j)
self.Hlocal.add_term(
hterms["hop_dag"].get_Hamiltonian(
-hop_coeff_half,
mask=self.get_mask([(R0 - 1) % self.n_sites]),
)
)
[docs]
def get_mask(self, sites_list):
"""Build a boolean mask selecting the given sites."""
mask = np.zeros(self.lvals, dtype=bool)
for site in sites_list:
mask[site] = True
return mask
[docs]
def local_parity_labels(self, wrt_site):
"""
Local action of pure spatial inversion (left-right swap)
on the 6d dressed-site SU(2) basis.
Basis labels:
0: V, J=(0,0)
1: V, J=(1/2,1/2) (link singlet)
2: (1/2,0,1/2) (matter + right link)
3: (1/2,1/2,0) (matter + left link)
4: P, J=(0,0)
5: P, J=(1/2,1/2) (link singlet)
Returns
-------
tuple
``(loc_perm, loc_phase)`` where ``loc_perm`` contains the mapped
local basis indices and ``loc_phase`` contains the corresponding
phases (``+1`` or ``-1``).
"""
if self.dim > 1:
raise ValueError(f"Does not work in D={self.dim}>1")
if self.spin > 0.5:
raise ValueError(f"Does not work spin j={self.spin}>0.5")
logger.info("----------------------------------------------------")
if wrt_site:
# local map for "parity" (P_loc)
loc_perm = np.array([0, 1, 3, 2, 4, 5], dtype=np.int32)
loc_phase = np.array([1, -1, 1, -1, 1, -1], dtype=np.int32)
else:
# local map for "parity + particle-hole" (or C ∘ P_loc)
loc_perm = np.array([4, 5, 3, 2, 0, 1], dtype=np.int32)
loc_phase = np.array([1, -1, 1, 1, 1, -1], dtype=np.int32)
return loc_perm, loc_phase
[docs]
def get_parity_inversion_operator(self, wrt_site):
"""Construct the parity inversion operator in the current sector."""
# INVERSION SYMMETRY
if self.dim < 1:
raise NotImplementedError(f"Cannot apply PARITY in D={self.dim}>1")
# Acquire the phases and the permutation of the basis
if wrt_site:
inversion_wrt = f"site {self.n_sites//2-1}"
else:
inversion_wrt = f"bond ({self.n_sites//2-1}, {self.n_sites//2})"
loc_perm, loc_phase = self.local_parity_labels(wrt_site)
wrt_site = np.uint8(0 if wrt_site else 1)
logger.info("----------------------------------------------------")
logger.info("Parity symmetry: inversion wrt %s", inversion_wrt)
# Build the projector
r, c, d = build_parity_operator(
self.sector_configs, loc_perm, loc_phase, wrt_site
)
shape = (self.sector_dim, self.sector_dim)
logger.info("r %s c %s, d %s, %s", r.shape, c.shape, d.shape, shape)
self.parityOP = csr_matrix((d, (r, c)), shape=shape)
[docs]
def check_parity_inversion_operator(self):
"""Validate algebraic and dynamical properties of the parity operator."""
if self.momentum_basis is not None:
# Build the projector from the momentum sector to the global one
Pk = self._basis_Pk_as_csr()
# Check the commutator between the Hamiltonian H and parity P
err_PH = spnorm(
self.parityOP @ Pk @ self.H.Ham @ Pk.getH()
- Pk @ self.H.Ham @ Pk.getH() @ self.parityOP
)
else:
err_PH = spnorm(self.parityOP @ self.H.Ham - self.H.Ham @ self.parityOP)
logger.info("|PH-HP| = %s", err_PH)
# -----------------------------------------------------------
I = identity(self.sector_dim)
# P^2 = I
err_P2 = spnorm(self.parityOP @ self.parityOP - I)
logger.info("|P^2-I| = %s", err_P2)
# Hermitian: P^† = P
err_herm = spnorm(self.parityOP.getH() - self.parityOP)
logger.info("|P^†-P| = %s", err_herm)
# Unitary: P^† P = I
err_unit = spnorm(self.parityOP.getH() @ self.parityOP - I)
logger.info("|P^†P-I| = %s", err_unit)
# nnz per row
nnz_per_row = np.diff(self.parityOP.indptr)
# nnz per col
nnz_per_col = np.diff(self.parityOP.tocsc().indptr)
logger.info("row nnz counts: %s", np.unique(nnz_per_row))
logger.info("col nnz counts: %s", np.unique(nnz_per_col))
# values
unique_vals = np.unique(np.round(self.parityOP.data, 8))
logger.info("unique data values in P: %s", unique_vals)
if self.momentum_basis is None and self.sector_configs is not None:
for ii, cfg in enumerate(self.sector_configs):
s1 = self.get_qmb_state_from_configs([cfg])
Ss1 = QMB_state(s1, self.lvals, self.loc_dims)
Ps1 = QMB_state(self.parityOP @ s1, self.lvals, self.loc_dims)
exp_s1 = Ss1.expectation_value(self.H.Ham)
exp_Ps1 = Ps1.expectation_value(self.H.Ham)
Hpsi = self.H.Ham @ s1 # H|psi>
Ppsi = self.parityOP @ s1 # P|psi>
HPpsi = self.H.Ham @ Ppsi # HP|psi>
PHpsi = self.parityOP @ Hpsi # PH|psi>
diff = HPpsi - PHpsi
absnorm = np.linalg.norm(diff)
if np.abs(exp_Ps1 - exp_s1) > 1e-14:
logger.info("%s", exp_s1)
logger.info("%s", exp_Ps1)
raise ValueError(f"config {ii} {cfg} not symmetrized")
if np.abs(absnorm) > 1e-14:
logger.info("==============================================")
logger.info("State i")
Ss1.get_state_configurations(1e-3, self.sector_configs)
logger.info("State P|i>")
Ps1.get_state_configurations(1e-3, self.sector_configs)
A = QMB_state(HPpsi, self.lvals, self.loc_dims)
logger.info("State HP|i>")
A.get_state_configurations(1e-3, self.sector_configs)
B = QMB_state(PHpsi, self.lvals, self.loc_dims)
logger.info("State PH|i>")
B.get_state_configurations(1e-3, self.sector_configs)
logger.info("VDOT%s", np.vdot(PHpsi, HPpsi))
logger.info("Relative |HP-PH|=%s", absnorm)
raise ValueError(f"config {ii} {cfg} not symmetrized")
[docs]
def print_state_config(self, config, amplitude=None):
"""Log a readable per-site decomposition of an SU(2) basis configuration."""
logger.info("----------------------------------------------------")
msg = f"CONFIG {config}"
if amplitude is not None:
msg += f" |psi|^2={np.abs(amplitude)**2:.8f}"
logger.info(msg)
logger.info("----------------------------------------------------")
# Choose width (sign + digits).
max_abs = 1
max_abs = max(max_abs, int(abs(getattr(self, "spin", 1))))
max_abs = max(max_abs, int(abs(getattr(self, "background", 0))))
entry_width = len(str(max_abs)) + 2
logger.info("%s%s", " " * 19, self._format_header(entry_width))
for ii, cfg_idx in enumerate(config):
loc_basis_state = self.gauge_states_per_site[ii][cfg_idx]
state_str = "[" + " ".join(f"{val}" for val in loc_basis_state) + " ]"
logger.info("SITE %2d state %3d: %s", ii, cfg_idx, state_str)
def _local_state_labels(self) -> list[str]:
"""Return column labels used by :meth:`print_state_config`."""
labels: list[str] = []
# Background first (only if present in your effective gauge states)
if getattr(self, "background", 0) > 0:
labels.append("bg")
# Matter occupation (only if not pure theory)
if not getattr(self, "pure_theory", True):
labels.append("M")
# Links: -x,-y,-z,+x,+y,+z up to dim
dim = int(getattr(self, "dim", 0))
dirs = "xyz"[:dim]
labels += [f"-{d}" for d in dirs] + [f"+{d}" for d in dirs]
return labels
def _format_header(self, entry_width: int) -> str:
"""Format the header row for :meth:`print_state_config`."""
labels = self._local_state_labels()
header = "[" + " ".join(f"{lab:>{entry_width}s}" for lab in labels) + " ]"
return header
"""
def get_string_breaking_configs(self, finite_density=0):
Populate predefined string-breaking configurations for benchmark lattices.
Parameters
----------
finite_density : int, optional
Finite-density setting selecting one of the predefined datasets.
Returns
-------
None
Stores configurations in ``self.string_cfgs`` and related counters.
logger.info("finite density %s", finite_density)
if self.lvals == [5, 2]:
self.n_min_strings = 5
self.n_max_strings = 1
if finite_density == 0:
self.string_cfgs = {
"max0": np.array([6, 10, 2, 10, 1, 5, 3, 10, 3, 11], dtype=int),
"min0": np.array([7, 12, 3, 12, 1, 4, 0, 9, 0, 11], dtype=int),
"min1": np.array([7, 12, 3, 11, 0, 4, 0, 9, 1, 12], dtype=int),
"min2": np.array([7, 12, 2, 9, 0, 4, 0, 10, 2, 12], dtype=int),
"min3": np.array([7, 11, 0, 9, 0, 4, 1, 11, 2, 12], dtype=int),
"min4": np.array([6, 9, 0, 9, 0, 5, 2, 11, 2, 12], dtype=int),
}
elif finite_density == 2:
self.string_cfgs = {
"min0": np.array([7, 12, 12, 12, 1, 4, 0, 9, 0, 11], dtype=int),
"min1": np.array([7, 12, 12, 11, 0, 4, 0, 9, 1, 12], dtype=int),
"min2": np.array([7, 12, 11, 9, 0, 4, 0, 10, 2, 12], dtype=int),
"min3": np.array([7, 11, 9, 9, 0, 4, 1, 11, 2, 12], dtype=int),
"min4": np.array([6, 9, 9, 9, 0, 5, 2, 11, 2, 12], dtype=int),
}
elif finite_density == 4:
self.string_cfgs = {
"min0": np.array([7, 12, 12, 12, 5, 4, 0, 9, 0, 11], dtype=int),
"min1": np.array([7, 12, 12, 11, 4, 4, 0, 9, 1, 12], dtype=int),
"min2": np.array([7, 12, 11, 9, 4, 4, 0, 10, 2, 12], dtype=int),
"min3": np.array([7, 11, 9, 9, 4, 4, 1, 11, 2, 12], dtype=int),
"min4": np.array([6, 9, 9, 9, 4, 5, 2, 11, 2, 12], dtype=int),
}
elif finite_density == 6:
self.string_cfgs = {
"min0": np.array([12, 12, 12, 12, 5, 4, 0, 9, 0, 11], dtype=int),
"min1": np.array([12, 12, 12, 11, 4, 4, 0, 9, 1, 12], dtype=int),
"min2": np.array([12, 12, 11, 9, 4, 4, 0, 10, 2, 12], dtype=int),
"min3": np.array([12, 11, 9, 9, 4, 4, 1, 11, 2, 12], dtype=int),
"min4": np.array([11, 9, 9, 9, 4, 5, 2, 11, 2, 12], dtype=int),
}
elif self.lvals == [4, 3]:
self.n_min_strings = 10
self.n_max_strings = 24
self.string_cfgs = {
"min0": np.array([7, 12, 3, 5, 9, 0, 21, 1, 0, 9, 0, 11], dtype=int),
"min1": np.array([7, 12, 2, 4, 9, 0, 24, 2, 0, 9, 0, 11], dtype=int),
"min2": np.array([7, 12, 2, 4, 9, 0, 23, 0, 0, 9, 1, 12], dtype=int),
"min3": np.array([7, 11, 0, 4, 9, 3, 26, 2, 0, 9, 0, 11], dtype=int),
"min4": np.array([7, 11, 0, 4, 9, 3, 25, 0, 0, 9, 1, 12], dtype=int),
"min5": np.array([7, 11, 0, 4, 9, 2, 21, 0, 0, 10, 2, 12], dtype=int),
"min6": np.array([6, 9, 0, 4, 12, 5, 26, 2, 0, 9, 0, 11], dtype=int),
"min7": np.array([6, 9, 0, 4, 12, 5, 25, 0, 0, 9, 1, 12], dtype=int),
"min8": np.array([6, 9, 0, 4, 12, 4, 21, 0, 0, 10, 2, 12], dtype=int),
"min9": np.array([6, 9, 0, 4, 11, 0, 21, 0, 1, 11, 2, 12], dtype=int),
# --------------------------------------------------------------------
"max0": np.array([7, 12, 3, 5, 10, 5, 26, 3, 1, 11, 2, 12], dtype=int),
"max1": np.array([6, 10, 3, 5, 11, 3, 26, 3, 1, 11, 2, 12], dtype=int),
"max2": np.array([7, 11, 1, 5, 10, 6, 24, 3, 1, 11, 2, 12], dtype=int),
"max3": np.array([7, 12, 3, 5, 10, 4, 22, 3, 1, 12, 1, 12], dtype=int),
"max4": np.array([7, 12, 3, 5, 10, 5, 25, 1, 1, 11, 3, 11], dtype=int),
"max5": np.array([6, 10, 3, 5, 12, 7, 26, 3, 0, 10, 2, 12], dtype=int),
"max6": np.array([6, 10, 3, 5, 12, 8, 26, 3, 0, 10, 2, 12], dtype=int),
"max7": np.array([6, 10, 3, 5, 11, 2, 22, 3, 1, 12, 1, 12], dtype=int),
"max8": np.array([6, 10, 3, 5, 11, 3, 25, 1, 1, 11, 3, 11], dtype=int),
"max9": np.array([7, 11, 1, 5, 10, 6, 23, 1, 1, 11, 3, 11], dtype=int),
"max10": np.array([7, 12, 2, 4, 10, 5, 28, 2, 1, 11, 3, 11], dtype=int),
"max11": np.array([7, 12, 2, 4, 10, 5, 29, 2, 1, 11, 3, 11], dtype=int),
"max12": np.array([6, 10, 3, 5, 12, 7, 25, 1, 0, 10, 3, 11], dtype=int),
"max13": np.array([6, 10, 3, 5, 12, 8, 25, 1, 0, 10, 3, 11], dtype=int),
"max14": np.array([6, 9, 1, 5, 11, 1, 28, 3, 1, 12, 1, 12], dtype=int),
"max15": np.array([6, 9, 1, 5, 11, 1, 29, 3, 1, 12, 1, 12], dtype=int),
"max16": np.array([6, 10, 2, 4, 11, 3, 28, 2, 1, 11, 3, 11], dtype=int),
"max17": np.array([6, 10, 2, 4, 11, 3, 29, 2, 1, 11, 3, 11], dtype=int),
"max18": np.array([7, 11, 1, 5, 10, 7, 27, 1, 1, 12, 0, 11], dtype=int),
"max19": np.array([7, 11, 1, 5, 10, 8, 27, 1, 1, 12, 0, 11], dtype=int),
"max20": np.array([6, 10, 2, 4, 12, 7, 28, 2, 0, 10, 3, 11], dtype=int),
"max21": np.array([6, 10, 2, 4, 12, 8, 28, 2, 0, 10, 3, 11], dtype=int),
"max22": np.array([6, 10, 2, 4, 12, 7, 29, 2, 0, 10, 3, 11], dtype=int),
"max23": np.array([6, 10, 2, 4, 12, 8, 29, 2, 0, 10, 3, 11], dtype=int),
}
elif self.lvals == [3, 2]:
self.n_min_strings = 1
self.n_max_strings = 1
self.string_cfgs = {
"max0": np.array([25, 30, 2, 9, 9, 35], dtype=int),
"min0": np.array([27, 37, 2, 7, 0, 35], dtype=int),
}
elif self.lvals == [6, 2]:
self.n_min_strings = 6
self.n_max_strings = 0
if finite_density == 0:
self.string_cfgs = {
"min0": np.array([7, 12, 3, 12, 3, 5, 4, 0, 9, 0, 9, 6], dtype=int),
"min1": np.array(
[7, 12, 3, 12, 2, 4, 4, 0, 9, 0, 10, 7], dtype=int
),
"min2": np.array(
[7, 12, 3, 11, 0, 4, 4, 0, 9, 1, 11, 7], dtype=int
),
"min3": np.array(
[7, 12, 2, 9, 0, 4, 4, 0, 10, 2, 11, 7], dtype=int
),
"min4": np.array(
[7, 11, 0, 9, 0, 4, 4, 1, 11, 2, 11, 7], dtype=int
),
"min5": np.array([6, 9, 0, 9, 0, 4, 5, 2, 11, 2, 11, 7], dtype=int),
}
if finite_density == 2:
self.string_cfgs = {
"min0": np.array(
[7, 12, 12, 12, 3, 5, 4, 0, 9, 0, 9, 6], dtype=int
),
"min1": np.array(
[7, 12, 12, 12, 2, 4, 4, 0, 9, 0, 10, 7], dtype=int
),
"min2": np.array(
[7, 12, 12, 11, 0, 4, 4, 0, 9, 1, 11, 7], dtype=int
),
"min3": np.array(
[7, 12, 11, 9, 0, 4, 4, 0, 10, 2, 11, 7], dtype=int
),
"min4": np.array(
[7, 11, 9, 9, 0, 4, 4, 1, 11, 2, 11, 7], dtype=int
),
"min5": np.array([6, 9, 9, 9, 0, 4, 5, 2, 11, 2, 11, 7], dtype=int),
}
if finite_density == 4:
self.string_cfgs = {
"min0": np.array(
[7, 12, 12, 12, 12, 5, 4, 0, 9, 0, 9, 6], dtype=int
),
"min1": np.array(
[7, 12, 12, 12, 11, 4, 4, 0, 9, 0, 10, 7], dtype=int
),
"min2": np.array(
[7, 12, 12, 11, 9, 4, 4, 0, 9, 1, 11, 7], dtype=int
),
"min3": np.array(
[7, 12, 11, 9, 9, 4, 4, 0, 10, 2, 11, 7], dtype=int
),
"min4": np.array(
[7, 11, 9, 9, 9, 4, 4, 1, 11, 2, 11, 7], dtype=int
),
"min5": np.array([6, 9, 9, 9, 9, 4, 5, 2, 11, 2, 11, 7], dtype=int),
}
if finite_density == 6:
self.string_cfgs = {
"min0": np.array(
[7, 12, 12, 12, 12, 5, 4, 0, 9, 0, 9, 11], dtype=int
),
"min1": np.array(
[7, 12, 12, 12, 11, 4, 4, 0, 9, 0, 10, 12], dtype=int
),
"min2": np.array(
[7, 12, 12, 11, 9, 4, 4, 0, 9, 1, 11, 12], dtype=int
),
"min3": np.array(
[7, 12, 11, 9, 9, 4, 4, 0, 10, 2, 11, 12], dtype=int
),
"min4": np.array(
[7, 11, 9, 9, 9, 4, 4, 1, 11, 2, 11, 12], dtype=int
),
"min5": np.array(
[6, 9, 9, 9, 9, 4, 5, 2, 11, 2, 11, 12], dtype=int
),
}
if finite_density == 8:
self.string_cfgs = {
"min0": np.array(
[12, 12, 12, 12, 12, 5, 4, 0, 9, 0, 9, 11], dtype=int
),
"min1": np.array(
[12, 12, 12, 12, 11, 4, 4, 0, 9, 0, 10, 12], dtype=int
),
"min2": np.array(
[12, 12, 12, 11, 9, 4, 4, 0, 9, 1, 11, 12], dtype=int
),
"min3": np.array(
[12, 12, 11, 9, 9, 4, 4, 0, 10, 2, 11, 12], dtype=int
),
"min4": np.array(
[12, 11, 9, 9, 9, 4, 4, 1, 11, 2, 11, 12], dtype=int
),
"min5": np.array(
[11, 9, 9, 9, 9, 4, 5, 2, 11, 2, 11, 12], dtype=int
),
}
else:
msg = "String breaking in ED considered only for lvals=[5,2] or [4,3] or [3,2]"
raise ValueError(msg)
"""