"""Digital Flux Lattice (DFL) SU(2)-based model helpers."""
import logging
import numpy as np
from edlgt.modeling import (
LocalTerm,
NBodyTerm,
PlaquetteTerm,
TwoBodyTerm,
check_link_symmetry,
get_lattice_link_site_pairs,
staggered_mask,
)
from edlgt.operators import (
SU2_dressed_site_operators,
SU2_gauge_invariant_states,
SU2_gen_dressed_site_operators,
)
from edlgt.symmetries import (
get_state_configs,
get_symmetry_sector_generators,
symmetry_sector_configs,
)
from edlgt.tools.config_encoding import config_to_index_binarysearch
from .quantum_model import QuantumModel
logger = logging.getLogger(__name__)
__all__ = ["DFL_Model"]
[docs]
class DFL_Model(QuantumModel):
"""DFL model built from SU(2) dressed-site operators."""
def __init__(self, spin, pure_theory, ham_format, **kwargs):
"""Initialize the DFL model.
Parameters
----------
spin : float
Gauge-link spin representation.
pure_theory : bool
If ``True``, exclude matter fields.
background : int or list
Background-charge specification passed to the gauge-basis builder.
ham_format : str
Hamiltonian representation format.
**kwargs
Arguments forwarded to :class:`~edlgt.models.quantum_model.QuantumModel`.
"""
# Initialize base class with the common parameters
super().__init__(**kwargs)
# DFL uses SU(2) dressed-site operators that are intrinsically complex.
self.configure_dtype_mode(dtype_mode="complex", auto_mode="complex")
self.spin = spin
self.pure_theory = pure_theory
self.background = 0.5
self.staggered_basis = False
self.ham_format = ham_format
# -------------------------------------------------------------------------------
# Acquire gauge invariant basis and states
self.gauge_basis, self.gauge_states = SU2_gauge_invariant_states(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
)
# -------------------------------------------------------------------------------
# Acquire operators
if self.spin < 1:
ops = SU2_dressed_site_operators(
self.spin,
self.pure_theory,
lattice_dim=self.dim,
background=self.background,
)
else:
# Acquire operators
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)
# Rather than for SU2, here we do not select the symmetry sector
[docs]
def build_Hamiltonian(self, g, m=None, dtype_mode="auto"):
"""Assemble the standard DFL 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.
"""
self.configure_dtype_mode(dtype_mode=dtype_mode, auto_mode="complex")
logger.info("BUILDING HAMILTONIAN")
# Hamiltonian Coefficients
self.DFL_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)
if not np.isclose(self.coeffs["E"], 1e-10):
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
)
if not np.isclose(self.coeffs["B"], 1e-10):
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
)
if not np.isclose(self.coeffs["B"], 1e-10):
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
)
if not np.isclose(self.coeffs["B"], 1e-10):
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
h_terms["N_tot"] = LocalTerm(self.ops["N_tot"], "N_tot", **self.def_params)
for site in ["even", "odd"]:
if not np.isclose(self.coeffs["m"], 1e-10):
self.H.add_term(
h_terms["N_tot"].get_Hamiltonian(
self.coeffs[f"m_{site}"], staggered_mask(self.lvals, site)
)
)
# ---------------------------------------------------------------------
# HOPPING
for d in self.directions:
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["hopping"] = TwoBodyTerm(
d, op_list, op_names_list, **self.def_params
)
for site in ["even", "odd"]:
# Define the mask
mask = staggered_mask(self.lvals, site)
self.H.add_term(
h_terms["hopping"].get_Hamiltonian(
strength=self.coeffs[f"t{d}_{site}"],
add_dagger=True,
mask=mask,
)
)
self.H.build(format=self.ham_format)
[docs]
def build_gen_Hamiltonian(self, g, m=None, dtype_mode="auto"):
"""Assemble the generalized DFL 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.
"""
self.configure_dtype_mode(dtype_mode=dtype_mode, auto_mode="complex")
logger.info("BUILDING generalized HAMILTONIAN")
# Hamiltonian Coefficients
self.DFL_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 = ["xy", "xz", "yz"]
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"],
]
for ii, pdir in enumerate(plaquette_directions):
if (self.dim > 1 and ii == 0) or self.dim == 3:
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 DFL_Hamiltonian_couplings(self, g, m=None):
"""Set DFL Hamiltonian couplings from physical parameters.
Parameters
----------
g : float
Gauge coupling.
m : float, optional
Bare mass parameter.
Returns
-------
None
Stores couplings in ``self.coeffs``.
Notes
-----
The conventions used here follow the DFL normalization used in the
project code and may differ from other SU(2) Hamiltonian conventions.
"""
if self.dim == 1:
E = g / 2
B = 0
else:
E = g
B = -1 # -1/ (2 * g)
# Dictionary with Hamiltonian COEFFICIENTS
self.coeffs = {
"g": g,
"E": E, # ELECTRIC FIELD coupling
"B": B, # MAGNETIC FIELD coupling
}
if not self.pure_theory and m is not None:
t = 1
self.coeffs |= {
"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 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 get_background_charges_configs(self, logical_stag_basis):
"""Generate 1D background-charge-compatible reference configurations.
Parameters
----------
logical_stag_basis : int
Period used to alternate logical staggering sectors.
Returns
-------
tuple
``(A, BG)`` arrays containing basis-state labels and corresponding
background sectors.
"""
# NOTE! It works only in 1D
if len(self.lvals) > 1:
raise ValueError("SU2 background configs works only in 1D")
logical_stag_basis = validate_dfl_staggered_reference_lattice(
self.n_sites, logical_stag_basis
)
has_obc = self.has_obc[0]
def next_val(value, last):
"""Return the two possible next values based on the current value."""
if value in [0, 7]:
return [4, 11] if last else [0, 6]
elif value == 3:
# ``3`` is the left-edge state ``[1/2, 'V', 0, 1/2]`` in the
# reduced ``site_mx`` basis. It is not equivalent to the bulk
# vacuum-family states ``0`` and ``7`` because its outgoing
# right rishon carries ``j=1/2``. The first bulk site must
# therefore continue into the bulk states with incoming
# left-link ``j=1/2``.
return [5, 12] if last else [1, 7]
elif value in [1, 6]:
return [5, 12] if last else [1, 7]
elif value in [4, 12]:
return [0, 6] if last else [4, 11]
elif value in [5, 11]:
return [1, 7] if last else [5, 12]
else:
raise ValueError(f"got unexpected value {value}")
def next_right_edge_val(value):
"""Return the unique OBC right-edge state compatible with ``value``.
The last site lives in the reduced ``site_px`` basis, where the DFL
reference manifold only uses the pair-family states
- ``2`` with ``J_config = [0, 'P', 0, 0]``
- ``5`` with ``J_config = [1/2, 'P', 1/2, 0]``
Which one is allowed is fixed by the outgoing ``+x`` rishon of the
penultimate bulk state: bulk states with right-link ``j=0`` close
onto edge state ``2``, while those with right-link ``j=1/2`` close
onto edge state ``5``.
"""
if value in [0, 7, 4, 12]:
return 2
elif value in [1, 6, 5, 11]:
return 5
raise ValueError(f"got unexpected value {value}")
# In 1D the reference-state builder is an automaton:
# once the local state on site i is chosen, only two states on site i+1
# keep the local flux/matter pattern compatible with the DFL family.
#
# PBC starts from the four bulk seeds. OBC instead starts from the two
# left-boundary seeds of the reduced six-state edge basis and ends on a
# unique right-edge state fixed by the penultimate bulk configuration.
first_site_space = 2 if has_obc else 4
if has_obc:
loc_dims = np.array(
[first_site_space] + [2 for _ in range(self.n_sites - 2)] + [1],
dtype=int,
)
else:
loc_dims = np.array(
[first_site_space] + [2 for _ in range(self.n_sites - 1)], dtype=int
)
A = np.zeros((np.prod(loc_dims), self.n_sites), dtype=int)
local_configs = get_state_configs(loc_dims)
# Initial states and bg charge sector
if not has_obc:
init = np.array([0, 1, 6, 7], dtype=int)
else:
init = np.array([0, 3], dtype=int)
for tmp in range(len(A)):
local_config = local_configs[tmp]
# Save the config of the first site
A[tmp, 0] = init[local_config[0]]
for ii in range(1, self.n_sites, 1):
if has_obc and ii == self.n_sites - 1:
# The right boundary is not another bulk site: once the
# penultimate bulk state is fixed, Gauss law and the
# reduced edge basis determine a unique allowed state.
A[tmp, ii] = next_right_edge_val(A[tmp, ii - 1])
continue
# Every ``logical_stag_basis`` sites we switch to the companion
# local family. For ``logical_stag_basis=2`` this implements the
# familiar empty-empty / pair-pair alternation used in the DFL
# half-filled reference manifold.
last = ii % logical_stag_basis == 0
# Save the config of the ii-th site
A[tmp, ii] = next_val(A[tmp, ii - 1], last)[local_config[ii]]
# Filter rows if has_obc is False to make it compatible with PBC
if not has_obc:
# Close the automaton on the ring by matching the first and last
# bulk-state families.
mask1 = np.isin(A[:, 0], [0, 6]) & np.isin(A[:, -1], [4, 12])
mask2 = np.isin(A[:, 0], [1, 7]) & np.isin(A[:, -1], [5, 11])
mask = mask1 | mask2
else:
# OBC already closes onto the reduced right-edge basis explicitly in
# ``next_right_edge_val``. Keep a final sanity check so any future
# change to the automaton fails loudly instead of silently dropping
# the zero-background edge sector again.
mask = np.isin(A[:, -1], [2, 5])
A = A[mask]
if len(A) == 0:
return A, np.zeros((0, self.n_sites), dtype=int)
# Do not infer the background pattern from the binary automaton branch:
# at OBC the boundary basis is smaller, so the branch label no longer
# coincides with the physical background charge on the last site.
# Instead, read the background sector directly from the projected local
# ``bg`` operator that defines the actual string constraint used by the
# DFL sector-by-sector workflow.
if "bg" not in self.ops:
msg = "DFL background-charge configs require the local 'bg' operator"
raise ValueError(msg)
bg_diags = np.real(np.diagonal(self.ops["bg"], axis1=1, axis2=2))
bg_sector_value = float(np.max(np.abs(bg_diags)))
if np.isclose(bg_sector_value, 0.0):
msg = "Cannot extract background sectors: 'bg' operator is identically zero"
raise ValueError(msg)
BG = np.zeros_like(A, dtype=int)
for site_idx in range(self.n_sites):
ratios = bg_diags[site_idx, A[:, site_idx]] / bg_sector_value
BG[:, site_idx] = np.rint(ratios).astype(int)
# When the caller has already built the half-filled global+link sector
# (the normal DFL workflow does this first), discard any automaton branch
# that is not actually part of that many-body sector. This is what fixes
# the spurious OBC left-edge family generated by the old bookkeeping.
if self.sector_configs is not None:
valid_mask = np.zeros(len(A), dtype=bool)
for row_idx, config in enumerate(A):
valid_mask[row_idx] = (
config_to_index_binarysearch(config, self.sector_configs) >= 0
)
A = A[valid_mask]
BG = BG[valid_mask]
return A, BG
[docs]
def get_bgsector_groups(self, logical_stag_basis):
"""Enumerate the unique background-charge sectors and their statistical
weights from the staggered DFL reference manifold.
Parameters
----------
logical_stag_basis : int
Period passed through to :meth:`get_background_charges_configs`.
Returns
-------
tuple
``(bg_configs, bg_sectors, unique_bg_sectors, grouped_configs, weights)``
with one entry per sector in ``unique_bg_sectors`` and with
``weights`` summing to 1.
"""
# Raw (config, bg-vector) pairs from the automaton-based builder.
bg_configs, bg_sectors = self.get_background_charges_configs(logical_stag_basis)
# Deduplicate bg vectors and remember which raw configs map to each.
unique_bg_sectors, indices = np.unique(bg_sectors, axis=0, return_inverse=True)
# Bucket the configs by their sector index.
grouped_configs = [
bg_configs[indices == ii] for ii in range(len(unique_bg_sectors))
]
# Equal weight per config -> probability of each sector after normalisation.
weights = np.array([len(cfgs) for cfgs in grouped_configs], dtype=float)
weights /= weights.sum()
return bg_configs, bg_sectors, unique_bg_sectors, grouped_configs, weights
[docs]
def build_single_bgsector_configs(self, bg_sector):
"""Build the Hilbert-space configs of *one* bg-charge sector in a
single string-aware pass.
Combines the global+link constraints with the bg-string constraint up
front so the iterative builder prunes site-by-site and never
materialises the full DFL sector_configs (which OOMs at L=12).
Parameters
----------
bg_sector : ndarray
Background-charge vector identifying the sector.
Returns
-------
ndarray
Configurations belonging to the requested bg sector.
"""
# Standard DFL gauge constraints (global U(1) + link Gauss law).
global_ops, global_sectors, link_ops, link_sectors, pair_list = (
self.get_symmetry_inputs()
)
# Wrap the local bg operator as a per-site "string" generator.
bg_global_ops = get_symmetry_sector_generators(
[self.ops["bg"]], action="global"
)
# Scale factor that turns an integer bg vector into the actual diagonal
# value the constraint should match.
bg_sector_value = float(np.max(bg_global_ops))
return symmetry_sector_configs(
loc_dims=self.loc_dims,
glob_op_diags=global_ops,
glob_sectors=global_sectors,
sym_type_flag="U",
link_op_diags=link_ops,
link_sectors=link_sectors,
pair_list=pair_list,
string_op_diags=bg_global_ops,
string_sectors=bg_sector_value * np.atleast_2d(bg_sector).astype(float),
)
[docs]
def prepare_sector_split(self, logical_stag_basis):
"""Build the full global+link ``sector_configs`` once and enumerate the
bg sectors that live inside it.
Heavy: at L=12 the global+link sector has ~370 M entries and OOMs.
Use :meth:`build_single_bgsector_configs` per-sector when the
in-process loop is not an option.
Parameters
----------
logical_stag_basis : int
Period for the staggered reference manifold.
Returns
-------
tuple
``(global_ops, global_sectors, link_ops, link_sectors, pair_list,
bg_global_ops, bg_sector_value, unique_bg_sectors,
grouped_configs, weights)`` — every piece the in-process sector
loop needs.
Side Effects
------------
Sets ``self.sector_configs`` to the full global+link sector.
"""
# Standard DFL gauge constraints.
global_ops, global_sectors, link_ops, link_sectors, pair_list = (
self.get_symmetry_inputs()
)
# Local bg operator -> diagonals + max value (used later to scale the
# per-sector target as bg_sector_value * bg_sector).
bg_global_ops = get_symmetry_sector_generators(
[self.ops["bg"]], action="global"
)
bg_sector_value = float(np.max(bg_global_ops))
# Build the full global+link sector once. This is the expensive step.
self.sector_configs = symmetry_sector_configs(
loc_dims=self.loc_dims,
glob_op_diags=global_ops,
glob_sectors=global_sectors,
sym_type_flag="U",
link_op_diags=link_ops,
link_sectors=link_sectors,
pair_list=pair_list,
)
# Enumerate the bg sectors that the in-process loop will iterate over.
_, _, unique_bg_sectors, grouped_configs, weights = self.get_bgsector_groups(
logical_stag_basis
)
return (
global_ops,
global_sectors,
link_ops,
link_sectors,
pair_list,
bg_global_ops,
bg_sector_value,
unique_bg_sectors,
grouped_configs,
weights,
)
[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
-----
For the 1D truncated SU(2) dressed-site DFL model the off-diagonal
entries are
``C_ij = 1/4 <Qpx_dag(i) W(i+1) ... W(j-1) Qmx(j)>`` for ``i < j``
and the diagonal is ``C_ii = <N_tot(i)> / 2``.
See :meth:`~edlgt.models.SU2_Model.get_fermionic_string_correlator`
for a detailed derivation.
"""
if self.dim != 1 or self.pure_theory:
raise NotImplementedError(
"Fermionic string correlator is implemented only for 1D DFL with matter."
)
required_ops = {"N_tot", "Qpx_dag", "Qmx", "W"}
if not required_ops.issubset(self.ops):
raise NotImplementedError(
"Fermionic string correlator requires N_tot, Qpx_dag, Qmx, and W "
"in the operator set."
)
if not all(self.has_obc):
raise NotImplementedError(
"Fermionic string correlator is currently implemented only with OBC."
)
if not hasattr(self, "def_params"):
self.default_params()
target_state = self._get_measurement_state(
state=state, state_index=state_index, dynamics=dynamics
)
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):
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)]
op_names_list = ["Qpx_dag"] + ["W"] * (dist - 1) + ["Qmx"]
op_list = [self.ops[name] for name in op_names_list]
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 = False,
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("NG = %.16f", ng_value)
return ng_value, eigvals, corr
[docs]
def print_state_config(self, config):
"""Log a human-readable decomposition of a DFL basis configuration."""
logger.info("----------------------------------------------------")
logger.info("SINGLETS IN CONFIG %s", config)
logger.info("----------------------------------------------------")
for ii, cfg_idx in enumerate(config):
msg = f"site {ii} state {cfg_idx}"
lattice_label = self.lattice_labels[ii]
local_basis_state = self.gauge_states[lattice_label][cfg_idx]
local_basis_state.display_singlet(msg)
def validate_dfl_staggered_reference_lattice(
n_sites: int, logical_stag_basis: int
) -> int:
"""Validate the DFL reference-state staggering period.
The DFL initial-state automaton alternates between two complementary local
state families, and each family is held for ``logical_stag_basis`` sites.
Half-filled reference states therefore exist only when the chain length is
a multiple of ``2 * logical_stag_basis``.
"""
try:
stag_basis = int(logical_stag_basis)
except (TypeError, ValueError) as exc:
raise ValueError(
f"logical_stag_basis must be a positive integer: got {logical_stag_basis}"
) from exc
if stag_basis < 1 or stag_basis != logical_stag_basis:
raise ValueError(
f"logical_stag_basis must be a positive integer: got {logical_stag_basis}"
)
staggering_period = 2 * stag_basis
if n_sites % staggering_period != 0:
msg = (
"DFL half-filled reference states require n_sites to be a multiple "
f"of 2 * logical_stag_basis; got n_sites={n_sites}, "
f"logical_stag_basis={stag_basis}."
)
raise ValueError(msg)
return stag_basis