"""Bose-Hubbard model built on top of :class:`edlgt.models.quantum_model.QuantumModel`."""
import logging
import numpy as np
from edlgt.modeling import LocalTerm, TwoBodyTerm, zig_zag
from edlgt.operators import bose_operators
from .quantum_model import QuantumModel
logger = logging.getLogger(__name__)
__all__ = ["BoseHubbard_Model"]
[docs]
class BoseHubbard_Model(QuantumModel):
"""Bose-Hubbard lattice model with particle-number symmetry reduction."""
def __init__(self, n_max, sectors=None, **kwargs):
"""Initialize the Bose-Hubbard model.
Parameters
----------
n_max : int
Maximum onsite boson occupation.
sectors : list or None, optional
Global particle-number sector labels ``[N_total]``. If ``None``, no
symmetry reduction is applied.
**kwargs
Arguments forwarded to :class:`~edlgt.models.quantum_model.QuantumModel`.
"""
# Initialize base class with the common parameters
super().__init__(**kwargs)
self.n_max = n_max
self.coeffs = {}
# -------------------------------------------------------------------------------
# Acquire operators and project onto the local Hilbert space
ops = bose_operators(self.n_max)
self.project_operators(ops)
# -------------------------------------------------------------------------------
# GLOBAL SYMMETRIES
if sectors is not None:
global_ops = [self.ops["N"]]
global_sectors = sectors
else:
global_ops = None
global_sectors = None
# GET SYMMETRY SECTOR
self.get_abelian_symmetry_sector(
global_ops=global_ops,
global_sectors=global_sectors,
)
# DEFAULT PARAMS
self.default_params()
[docs]
def build_Hamiltonian(self, coeffs):
"""Assemble the Bose-Hubbard Hamiltonian.
Parameters
----------
coeffs : dict
Coupling dictionary with keys ``"t"`` (hopping amplitude)
and ``"U"`` (on-site interaction strength). Optional keys are
``"mu"`` for a uniform chemical potential and ``"local_potential"``
(or legacy ``"noise"``) for a site-resolved potential profile.
"""
logger.info("BUILDING HAMILTONIAN")
# Hamiltonian Coefficients
self.coeffs = coeffs
h_terms = {}
# -------------------------------------------------------------------------------
# NEAREST NEIGHBOR HOPPING
for axis in self.directions:
op_names_list = ["b_dagger", "b"]
op_list = [self.ops[op] for op in op_names_list]
h_terms[f"NN_{axis}"] = TwoBodyTerm(
axis=axis,
op_list=op_list,
op_names_list=op_names_list,
**self.def_params,
)
self.H.add_term(
h_terms[f"NN_{axis}"].get_Hamiltonian(
strength=-self.coeffs["t"], add_dagger=True
)
)
# -------------------------------------------------------------------------------
# ON-SITE INTERACTION
op_name = "N2"
h_terms[op_name] = LocalTerm(
operator=self.ops[op_name], op_name=op_name, **self.def_params
)
self.H.add_term(
h_terms[op_name].get_Hamiltonian(strength=0.5 * self.coeffs["U"])
)
op_name = "N"
h_terms[op_name] = LocalTerm(
operator=self.ops[op_name], op_name=op_name, **self.def_params
)
self.H.add_term(
h_terms[op_name].get_Hamiltonian(strength=-0.5 * self.coeffs["U"])
)
# -------------------------------------------------------------------------------
# OPTIONAL CHEMICAL POTENTIAL / SITE-RESOLVED POTENTIAL
chemical_potential = self.coeffs.get("mu", None)
if chemical_potential is not None:
self.H.add_term(
h_terms["N"].get_Hamiltonian(strength=-float(chemical_potential))
)
local_potential = self.coeffs.get("local_potential", self.coeffs.get("noise"))
if local_potential is not None:
local_potential = np.asarray(local_potential, dtype=float)
if local_potential.ndim == 0:
local_potential = np.full(self.n_sites, float(local_potential))
elif local_potential.shape != (self.n_sites,):
raise ValueError(
"local_potential must be scalar or have shape "
f"({self.n_sites},), got {local_potential.shape}"
)
for ii, potential in enumerate(local_potential):
if np.isclose(potential, 0.0):
continue
mask = np.zeros(self.n_sites, dtype=bool)
mask[ii] = True
self.H.add_term(
h_terms["N"].get_Hamiltonian(strength=potential, mask=mask)
)
# -------------------------------------------------------------------------------
self.H.build(self.ham_format)
[docs]
def print_state_config(self, config, amplitude=None):
"""Log a Bose-Hubbard occupation configuration site by site."""
if amplitude is not None:
logger.info(
"Config %s amplitude %.9f + %.9fi",
config,
np.real(amplitude),
np.imag(amplitude),
)
else:
logger.info("Config %s", config)
for site_index, occupation in enumerate(config):
logger.info(
"site %s coords %s occupation %s",
site_index,
tuple(zig_zag(self.lvals, site_index)),
occupation,
)