Core

binary_to_MPS(n, L, sites) -> MPS

Return the computational-basis state |n> as an L-site MPS, where n is encoded in big-endian binary across the L qubit sites.

Example

sites = siteinds("Qubit", 4)
psi   = binary_to_MPS(5, 4, sites)   # |0101>

constant_mps(sites, value) -> MPS

Rank-1 MPS whose amplitude is value on every computational-basis state.

custom_mpo(mps, new_sites) -> MPO

Convert a 2N-site MPS produced by QTCI on an interleaved 2D quantics grid into an N-site MPO by contracting each pair of tensors (2i-1, 2i) and mapping old site indices to new_sites[i].

The first index of each pair becomes the bra (primed) site and the second becomes the ket (unprimed) site, consistent with ITensors MPO conventions.

custom_mps(qtt, sites) -> MPS

Replace the site indices of an MPS obtained from a 1D TCI tensor train with the physical sites of the target system.

eval_mps(A, n) -> Real

Evaluate a profile MPS at the 0-indexed basis coordinate n using the big-endian convention of binary_to_MPS. Do not use _eval_diag_mps for this purpose — that helper uses LSB-first ordering for the QFT momentum convention.

eval_mps_spatial(A::MPS; num_x, num_avg, x_start, x_end, x_groups,
                 box_half, Lx) -> (values, centers, groups)

Higher-level spatial sampler for a profile MPS such as an SCF occupation/density profile (res.rho_up). It mirrors get_ldos_spatial's num_x / num_avg / x_groups / box_half sampling-and-averaging API, but evaluates the MPS directly with eval_mps instead of running a KPM recursion — so it is cheap enough to sweep a very large system by sampling a grid of positions and averaging, rather than evaluating all 2^L sites.

For each sampled group of (1-indexed) site coordinates the returned value is the mean of eval_mps(A, x-1) over that group. With box_half > 0 each sampled position is expanded into a (2·box_half+1)² neighborhood on the 2D grid (periodic wrap), exactly like get_ldos_spatial; this needs the 2D layout, taken from Lx (defaults to L÷2, with Ly = L - Lx).

Keyword arguments

• num_x : number of sampled grid positions (default: all 2^L sites).
• num_avg : sub-positions averaged per grid point along the 1D index (stride).
• x_start, x_end : 1-indexed sampling window (default 1 … 2^L).
• x_groups : explicit groups — a vector of site indices (one per group) or a vector of vectors (each averaged). Overrides num_x/num_avg/x_start/x_end.
• box_half : 2D neighborhood half-width for averaging (0 = no box averaging).
• Lx : number of x qubits for the 2D layout (default L÷2).

Returns

• values : Vector{Float64}, the averaged MPS value per group.
• centers : Vector{Int}, the 1-indexed center site of each sampled group.
• groups : Vector{Vector{Int}}, the site indices averaged over per group.

For a 2D map, the center (ix, iy) of group g is ix = (centers[g]-1) % 2^Lx, iy = (centers[g]-1) ÷ 2^Lx.

fix_sites(mpo, sites) -> MPO

Replace the site indices of mpo (typically built from a TCI tensor train whose indices do not match the system's physical sites) with sites. Modifies mpo in-place and returns it.

fused_mpo(mps, new_sites) -> MPO

Convert an N-site MPS with dim-4 physical indices produced by QTCI on a :fused 2D quantics grid into an N-site MPO. Each dim-4 physical index encodes one (bra-bit, ket-bit) pair; a combiner splits it into new_sites[i]' (bra) and new_sites[i] (ket).

get_diagonal_mpo(L, sites, f; type=Float64, tol=1e-8) -> MPO

Build a diagonal MPO with entry f(x) at 1-indexed site x ∈ {1, ..., 2^L}. Wraps get_mpo with a 1→0 index shift for backward compatibility.

Example

pot = get_diagonal_mpo(L, sites, x -> 0.01 * x)

get_matrix(mpo, sites) -> Matrix{ComplexF64}
get_matrix(mpo, L, sites) -> Matrix{ComplexF64}   (L ignored)

Return the full D x D dense matrix of mpo, where D = prod(dim(s_i)). Works for any site types (Qubit, Layer, Spin, Nambu, ...). Feasible only for small systems (D <= 512).

get_mpo(L, sites, f; type=Float64, tol=1e-8, kwargs...) -> MPO

Compress a function into an MPO via QTCI using 0-indexed coordinates n = 0, ..., 2^L-1. Dispatches on arity:

• 1-argument f(n): builds a diagonal MPO (on-site potential).
• 2-argument f(i, j): builds a full interaction MPO by 2D QTCI, pairing the interleaved quantics legs into bra/ket via custom_mpo.

get_mps(L, sites, f; type=Float64, tol=1e-8) -> MPS

Compress a scalar profile f(n) on n = 0, ..., 2^L-1 into an MPS via QTCI.

hadamard_mpo(A, B, out_sites; maxdim=100, cutoff=1e-8) -> MPO

Site-wise Hadamard product of two MPOs: C[i,j] = A[i,j] * B[i,j]. out_sites may be any index set, including A/B's own sites; fresh working indices are created internally and the result is remapped to out_sites.

matrix_checker(mpo, sites, i, j) -> Number
matrix_checker(mpo, L, sites, i, j) -> Number   (L ignored)

Return the matrix element <i|mpo|j>. Works for any site types (Qubit, Layer, Spin, Nambu, ...). Intended for small-system validation.

mps2mpo(L, sites, density_mps) -> MPO

Convert a diagonal MPS (a function sampled on computational-basis states) into a diagonal MPO by calling Quantics._asdiagonal on each site tensor.

mps_to_diagonal_mpo(mps, sites) -> MPO

Convert an MPS into a diagonal MPO on sites by replacing each physical index with a bra-ket pair tied by a 3-leg delta. Used to convert the output of a 2D QTCI (encoded as a flat MPS) into a diagonal MPO on the interleaved (e.g. electron-hole) site space.

postpend_op(H_mpo, s, mat)        -> MPO   explicit matrix
postpend_op(H_mpo, s, op::Symbol) -> MPO   named op (Spin / Nambu index)
postpend_op(H_mpo, s, k, l)       -> MPO   sparse |k><l|  (Layer / any index)
postpend_op(H_mpo, s, k)          -> MPO   projector |k><k|

Append a single-site operator on index s to the end of H_mpo, extending it from L sites to L+1 sites. The returned MPO has site indices [original...; s].

Symmetric counterpart of prepend_op; dispatch rules are identical.

prepend_op(H_mpo, s, mat)       -> MPO   explicit matrix
prepend_op(H_mpo, s, op::Symbol)-> MPO   named op (Spin / Nambu index)
prepend_op(H_mpo, s, k, l)      -> MPO   sparse |k><l|  (Layer / any index)
prepend_op(H_mpo, s, k)         -> MPO   projector |k><k|

Prepend a single-site operator on index s to H_mpo, extending it from L sites to L+1 sites. The returned MPO has site indices [s; original...].

Dispatch rules

• Matrix form: mat[i,j] = <i|op|j> (1-indexed). Element type is preserved.
• Symbol form: named operator looked up by the type of s (tag "Spin" or "Nambu"). Defined in Supercond_tk.jl after the op dictionaries.
• Integer pair (k, l): places a single 1 at row k, col l in a dim(s) x dim(s) zero matrix. Covers layer hops and projectors for any dimension Layer index.
• Single integer k: shorthand for the projector |k><k|.

replace_sites(MPOin, newsites) -> MPO

Replace the physical (bra + ket) indices of each site in MPOin with the corresponding index from newsites, preserving prime levels.

rms_error(a, b) -> Float64

RMS distance between two MPS objects over all computational-basis states.

spatial_sampling_plan(L; Lx, grid, reduce, n_sub, num_x, num_y, num_avg,
                      x_start, x_end, xwin, ywin, x_groups, box_half, sublattice)
    -> (; centers, groups, resolve_sublattice, n_sub, stride_x, stride_y,
         grid, reduce, a, b)

Geometry-aware real-space sampling plan shared by every spatial sampler (eval_mps_spatial, get_ldos_spatial, get_ldos_spatial_gpu, get_scf_magnetization_gpu). It decides where to sample, how each output pixel reduces the cells under it (reduce), and — for multi-atom unit cells — whether to resolve or average the sublattice.

The three sampling procedures (reduce)

A spatial map of a 2^Lx × 2^Ly-unit-cell system at a coarse output resolution can reduce the cells beneath each pixel in three qualitatively different ways. The right choice depends on whether the quantity is smooth on the large scale (e.g. a Chern marker, an SCF density envelope) or a thin feature on a flat background (e.g. in-gap edge/domain-wall LDOS, width ξ ≪ system size).

1. :point (default) — point / box sampling. Lay out num_x[×num_y] sample positions and read the profile at each one. With box_half > 0 each pixel is the mean over a (2·box_half+1)² neighbourhood (smoothing). Cost ∝ (number of pixels) × (box cells).

   Aliasing caveat. The pixels probe only the cells they land on (± box_half). On a grid coarser than a feature's width this misses thin features that fall between pixels: a domain-wall LDOS channel of width ξ sampled at stride s ≫ ξ is caught only on the rare pixel within box_half of it. Making the box tile the plane (box_half ≈ s/2) closes the gaps but then evaluates essentially every cell — i.e. full-resolution cost. Use :point for smooth quantities or for a fully-resolved zoom (grid=true + a small window).

2. :block — block integration (gap-free coarse-graining). Partition the system into num_x × num_y equal blocks (num_x = 2^a, num_y = 2^b, powers of two) and report, per pixel, the sum over its whole block. This is computed by tracing out the low-order position bits (contracting the within-block bits of the profile MPS with [1,1] and keeping the a + b high-order block bits) — a partial contraction, not a per-cell sweep, so the cost is independent of block size and scales to Lx, Ly ≈ 14+.

   Because every cell belongs to exactly one block, a thin feature cannot fall between pixels — whichever blocks it threads light up, on an otherwise dark (gapped) background. This is the tool for imaging edge / domain-wall networks on a heavily downsampled map. Block centres are reported in centers; the per-axis block widths are stride_x = 2^(Lx-a), stride_y = 2^(Ly-b). If Lx is omitted, the same mode becomes a 1D partition into num_x = 2^a contiguous intervals with stride_x = 2^(L-a).

The fields reduce, a, b echo the chosen mode back to the caller; for :point they are (:point, 0, 0).

Sublattice resolve vs average

For a multi-atom unit cell (n_sub > 1) the plan also decides whether to resolve the sublattice (one output column per atom) or average it (one value per unit cell, atoms traced out), via sublattice:

• :auto (default) — resolve only at the atomic scale: consecutive samples are adjacent unit cells (:point with stride == 1 and box_half == 0). Otherwise (coarse grid, box_half > 0, or any :block map) average, since the intra-cell sublattice is below the sampling resolution.
• :resolve / :average force the choice. n_sub == 1 is always false.

Layout (:point mode)

groups/centers are 1-indexed unit-cell indices with n = ix + iy·2^Lx.

• grid=false (default) — centers on a 1D linear sweep of the row-major index (x_start/x_end, num_x points, num_avg sub-probes per block). Stride dx = window ÷ num_x. Lx is used only for the optional box_half neighbourhood.
• grid=true (needs Lx) — centers on a 2D xy grid of num_x × num_y points over the unit-cell window xwin=(ix0,ix1), ywin=(iy0,iy1) (0-indexed; default full system). Per-axis strides Nx_win÷num_x, Ny_win÷num_y.

x_groups overrides the :point layout entirely; the stride is then unknown, so :auto resolves (treats it as atomic) unless box_half > 0.

to_binary_vector(n, L) -> Vector{String}

Convert non-negative integer n to a length-L vector of "0"/"1" strings (big-endian), suitable as state labels for MPS(sites, state).

Example

to_binary_vector(5, 4)   # -> ["0", "1", "0", "1"]

circular_mod(i, L, L_chain, k) -> Float64

Circularly symmetric on-site modulation cos(k * r) where r is the distance from the centre of the 2D lattice.

compose_power(base, nn; side=:right, apply_kwargs=NamedTuple()) -> MPO

Compose base with itself nn times using exponentiation-by-squaring (O(log n) applies). Replaces the old arbitarty_offline helper which used O(n) sequential applies.

• side = :right: acc = apply(acc, base) at each set bit
• side = :left: acc = apply(base, acc) at each set bit

apply_kwargs (e.g. (; cutoff=1e-8, maxdim=200)) are forwarded to every apply call. nn = 0 returns the identity MPO; nn = 1 returns base unchanged.

hopping2MPO(f, N, sites; tol=1e-8, initial_positions=[], type=Float64,
            unfoldingscheme=:interleaved) -> MPO

Compress an arbitrary NxN hopping matrix H[i,j] = f(i,j) into an MPO using Quantics Tensor Cross Interpolation on a 2D quantics grid (N must be a power of 2).

unfoldingscheme controls the bit ordering of the 2D quantics grid:

• :interleaved (default): row and column bits alternate: rL cL ... r1 c1
• :fused: all row bits first, then column bits: rL ... r1 cL ... c1

initial_positions seeds the TCI pivots; useful when the matrix has known structure (e.g. near-diagonal for short-time propagators).

kineticNNN(L, sites, hopping, nn; apply_kwargs=NamedTuple()) -> MPO

Build a kinetic MPO for a 1D chain with a spatially varying hopping field encoded as the diagonal MPO hopping, and a neighbor reach controlled by nn.

Construction uses the shared shift-hopping helper: hopping * shift(nn) + shift(nn)' * dag(hopping).

boundary=:open is the default. Use boundary=:periodic (or bc=:pbc) for a cyclic shift on the quantics chain.

kinetic_1d_nn(L, sites; boundary=:open) -> MPO

Build the nearest-neighbour hopping MPO for a 1D chain of 2^L sites in the quantics binary representation. Hopping amplitude = 1; scale by multiplying the result. The default boundary=:open preserves the package's open-chain convention; boundary=:periodic adds the wrap-around bond.

kinetic_1d_nn_custom(L, sites, hopping; boundary=:open) -> MPO

Nearest-neighbour 1D kinetic MPO with a site-dependent hopping encoded as a diagonal MPO hopping. Useful for spatially varying hopping amplitudes (e.g. SSH model, quasicrystals).

qtci_matrix_to_MPO(A_fun, L, sites; tol=1e-8, type=Float64,
                   initial_positions=[]) -> MPO

Like hopping2MPO but works with a (2^L)x(2^L) matrix function and applies an extra truncation step with maxdim=20.

qtt_mpo(L, xvals, sites, func; tol_quantics=1e-8, maxbonddim_quantics=50) -> MPO

Compress a scalar function func(x) evaluated on the explicit integer grid xvals (typically 0:2^L-1) into a diagonal MPO via Quantics Tensor Cross Interpolation.

The result is diag(func(0), func(1), ..., func(2^L-1)) stored as an L-site MPO. Use this to encode spatially varying on-site potentials or hopping amplitudes as diagonal MPOs for use with kineticNNN and the 2D kinetic builders.

xvals = 0:2^L-1 for a 1D chain of 2^L sites xvals = 0:Nx*Ny-1 for a row-major flattened 2D grid

See also get_diagonal_mpo in utils.jl for a simpler 1-based-index wrapper.

quasicrystal_modulation_30deg(i, L, L_chain, k, p) -> Float64

On-site modulation for a p-fold quasicrystal pattern at wavevector k, centred on the middle of the 2D lattice.

TBHamiltonian

Central object representing a tight-binding Hamiltonian in the quantics MPO framework.

Fields

Core

• L : number of position qubit sites (log₂ of the physical system size)
• N : number of physical sites / unit cells (2^L)
• sites : ITensor site indices (position qubits + any auxiliary DOF indices)
• mpo : accumulated Hamiltonian as an ITensor MPO
• geometry : function i -> position_vector (1-indexed); nothing for implicit 1D

KPM spectral bounds

• scale : energy half-bandwidth; H/scale has spectrum in [-1, 1]. Set analytically at construction for standard geometries. Reset to 0.0 by _invalidate_cache! after any modification — _ensure_scale! then re-estimates it on demand via DMRG.
• center : spectral centre; 0.0 for particle-hole symmetric Hamiltonians.

Auxiliary DOF indices (nothing until the corresponding add_*! is called)

• spin_s : dim-2 spin index (prepended by add_spin!)
• nambu_s : dim-2 Nambu index (prepended by add_superconductivity!)
• layer_s : dim-k layer index (set by bilayer/multilayer constructors)
• sublattice_s : dim-k sublattice index (set by kagomé/Lieb/honeycomb/dice constructors)
• aux_side : :pre or :post — position of the outermost aux index in sites

Lazy caches (cleared by _invalidate_cache! whenever mpo changes)

• _tn_cache : MPO Chebyshev list; set by KPM_Tn(H, N; mode=:mpo)
• _tn_mps_cache : MPS Chebyshev state list; set by KPM_Tn(H, N; mode=:mps, psi0=…)
• _tn_Ncheb : order of the cached Chebyshev expansion
• _density_cache : cached density-matrix MPO

Stored interaction MPOs (set via add_interaction!)

• interaction_mpo : Hartree/CDW interaction kernel (fully scaled); used by get_scf(H, channel)
• fock_mpo : Fock/exchange interaction kernel (fully scaled)

Do not construct directly — use get_Hamiltonian.

_invalidate_cache!(H) -> H

Clear all cached intermediate results. Called automatically whenever mpo changes (e.g. add_hopping!, add_onsite!, add_zeeman!, add_superconductivity!).

Clears: _tn_cache, _tn_mps_cache, _tn_Ncheb, _density_cache. Resets scale and center to 0.0 so that _ensure_scale! re-estimates them via DMRG on the next KPM call.

_pos_sites(H) -> Vector{<:Index}

Return the L position-qubit indices, filtering out any auxiliary (spin, Nambu, layer) indices regardless of whether they sit at the front or back of H.sites.

add_hopping!(H, f; nn, sublat, sublat_from, sublat_to, maxdim, tol, ...) -> H

Add a hopping term to H.

f modes (selected automatically by type):

• f::Number — uniform nn-th-neighbour hopping (amplitude f, shell nn).
• f(i), 1-arg Function — spatially varying amplitude per site, QTCI-compressed.
• f(i,j), 2-arg Function — full N×N matrix, QTCI-compressed (nn ignored).

Sublattice keywords (H.sublattice_s must be set; call before add_spin! etc.)

• sublat=k — intra-sublattice: restrict the full kinetic term to sublattice k by postpending the projector |k⟩⟨k|. Works with all three f modes.

  add_hopping!(H, -0.1; sublat=1, nn=2)           # NNN within sublattice A

• sublat_from=k, sublat_to=l — inter-sublattice: adds f · K_u^nn ⊗ |l⟩⟨k| + conj(f) · K_d^nn ⊗ |k⟩⟨l| (Hermitian, f must be scalar). Use nn=0 for intra-cell bonds (no position shift; useful to correct individual bond amplitudes in an existing sublattice Hamiltonian):

  add_hopping!(H, -0.3; sublat_from=1, sublat_to=2, nn=1)   # A↔B NN inter-cell
  add_hopping!(H, δt;   sublat_from=2, sublat_to=3, nn=0)   # B↔C intra-cell only

Without sublattice keywords the function errors if any auxiliary DOF is already attached. Invalidates all caches.

add_interaction!(H, V; channel=:hartree, type=Float64, tol=1e-8) -> H

Store a pre-built interaction MPO in H for later use with get_scf(H, channel).

V can be:

• MPO — stored directly
• Number — scalar coupling; stored as V · Id on the position sites
• Function with 1 arg — on-site potential V(i), 0-indexed; stored as diagonal MPO
• Function with 2 args — interaction kernel V(i,j), 0-indexed; stored as full 2D MPO

channel:

• :hartree or :default → stored in H.interaction_mpo (used by :cdw and :magnetic SCF)
• :fock or :exchange → stored in H.fock_mpo

add_onsite!(H, f; layer=nothing, sublat=nothing, Lx=nothing,
            tol=1e-8, maxdim=nothing) -> H

Add a diagonal (on-site) potential to H.

f argument — same conventions as add_hopping_2D!

• f::Number — uniform constant; builds f · Id directly (no QTCI).
• f(n) — 1-arg function; n ∈ {0, …, N-1} is the 0-indexed unit-cell index.
• f(ix, iy) — 2-arg function; ix, iy are 0-indexed 2D coordinates (ix = n % Nx, iy = n ÷ Nx). Requires Lx= keyword so that Nx = 2^Lx.

sublat keyword (H.sublattice_s must be set)

• sublat=nothing (default) — applies f equally to all sublattices.

• sublat=k — restricts the potential to sublattice k only via |k⟩⟨k|.

  Canonical use-case — Semenoff mass on a honeycomb Hamiltonian:

  add_onsite!(H_hc, +M; sublat=1)   # +M on sublattice A
  add_onsite!(H_hc, -M; sublat=2)   # −M on sublattice B  → gap = 2M

  Works for any sublattice-carrying TBHamiltonian (kagome, Lieb, dice, …) and for spatially-varying f as well.

layer keyword (H.layer_s must be set)

• layer=nothing — apply the same onsite term to every layer.
• layer=k — apply only to layer k.
• layer=[...] — apply to the selected layers.

For layered sublattice Hamiltonians, the onsite term is first built on the position ⊗ sublattice monolayer space, then wrapped with the selected layer projectors.

Invalidates all caches.

add_soc!(H, λ; type=:rashba, direction=:z, tol=1e-8, maxdim=200) -> H

Add spin-orbit coupling to H. Calls add_spin! automatically if needed.

type:

• :rashba — nearest-neighbour Rashba SOC on the position chain: λ · (S_y ⊗ K_u − S_y ⊗ K_d) where K_u/K_d are the ±1 shift operators. λ must be a scalar. Breaks SU(2) spin symmetry while preserving time-reversal.
• :ising — diagonal Ising SOC λ(i) · S_z (equivalent to a position-dependent Zeeman along z; useful for Kane–Mele type models).
• :custom — arbitrary position-space MPO λ_mpo tensor-producted with the spin operator given by direction (:x, :y, or :z). λ may be a Number, a 1-arg Function λ(i), or a 2-arg Function λ(i,j) (the last compressed via TCI). For the result to be Hermitian, the position-space matrix must itself be Hermitian: λ(i,j) = conj(λ(j,i)). Diagonal and real-symmetric inputs satisfy this automatically.

Examples

add_soc!(H, 0.05)                             # Rashba λ=0.05
add_soc!(H, i -> 0.1*cos(2π*i/H.N); type=:ising)
add_soc!(H, (i,j)->...; type=:custom, direction=:y)

add_spin!(H; cutoff=1e-8, maxdim=200) -> H

Extend H to a spin-½ degenerate system by prepending a spin-½ index. The resulting Hamiltonian is I_spin ⊗ H (both spin sectors identical).

No-op if H is already spinful (H.spin_s !== nothing). Invalidates all caches.

add_superconductivity!(H, Δ; type=:swave, tol=1e-8, maxdim=200) -> H

Extend H to a Bogoliubov–de Gennes (BdG) Hamiltonian by prepending a Nambu (particle–hole) index.

The BdG structure is: HBdG = τz ⊗ Hkin + τ+ ⊗ Hpair + τ- ⊗ H_pair†

• Spinless + p-wave (type=:pwave, or auto-selected when spinless + :swave): H_pair = Δ·(K_forward − K_backward), the antisymmetric nearest-neighbour matrix required by Fermi statistics (Δ(i,j) = −Δ(j,i)). This is the Kitaev chain. Δ must be a Number.
• Spinful + s-wave (add_spin! called first, type=:swave): H_pair = (i·σ_y)_spin ⊗ Δ(r), the standard BCS singlet Cooper-pair operator. On-site (s-wave) pairing is allowed here because the antisymmetry is carried by the spin singlet factor i·σ_y. Δ can be a Number or 1-arg Function.
• Custom (type=:custom): arbitrary pairing matrix via 2-arg function Δ(i,j), compressed with TCI.

Note on spinless s-wave: on-site pairing is forbidden for spinless fermions (Δ(i,i) = 0 by antisymmetry). Calling with type=:swave on a spinless chain automatically redirects to :pwave (uniform Δ) or errors (spatially varying Δ).

type:

• :swave (default) — diagonal pairing for spinful chains; auto-redirects to :pwave for spinless chains when Δ isa Number
• :pwave — antisymmetric NN pairing Δ·(K_f − K_b) for spinless chains; Δ must be a Number
• :custom — arbitrary Δ(i,j); pass a 2-arg function

Errors if BdG has already been applied. Invalidates all caches.

Examples

add_superconductivity!(H_spinless, 0.1)              # auto p-wave (Kitaev chain)
add_superconductivity!(H_spinless, 0.1; type=:pwave) # explicit p-wave
add_superconductivity!(H_spinful,  0.1)              # singlet s-wave (spinful required)
add_superconductivity!(H_spinful,  i -> i < N÷2 ? 0.1 : 0.0)  # spatially varying s-wave
add_superconductivity!(H, (i,j) -> ...; type=:custom)          # general pairing

add_zeeman!(H, h; direction=:z, tol=1e-8, maxdim=200) -> H

Add a Zeeman coupling h · Sα to H. Calls add_spin! automatically if H is not yet spinful.

h can be:

• a Number — uniform field amplitude h₀
• a Function — spatially varying h(i), i ∈ {1, …, N} (1-indexed)

direction: :x, :y, or :z (default).

If add_superconductivity! was already called, the Zeeman term is wrapped in τ_z so it enters with opposite sign in the hole sector, as required in BdG.

Examples

add_zeeman!(H, 0.1)                         # uniform h = 0.1 along z
add_zeeman!(H, i -> 0.05 * sin(2π*i/H.N))  # oscillating field
add_zeeman!(H, 0.05; direction=:x)          # in-plane

central_index(geom, N) -> Int
central_index(H)       -> Int

Return the 1-indexed site index whose position is closest to the geometric centroid of the lattice. Accepts either a geometry function geom(i) and system size N, or a TBHamiltonian directly.

Errors if H.geometry is nothing.

get_Hamiltonian(geometry, params; L, [scale, tol, maxdim, kwargs...])
    -> TBHamiltonian

Build a TBHamiltonian from a named geometry and model parameters.

Supported geometry strings

For "kagome" and "lieb", L = Lx + Ly counts only the position qubits; the total atom count is 3 × 2^L. The sublattice index is stored in H.sublattice_s with H.aux_side = :post. H.geometry returns the full real-space position of each atom (1-indexed over all 3 × 2^L atoms).

Common keyword arguments

• L : number of qubit sites (system size = 2^L)
• scale : energy half-bandwidth for KPM normalisation (estimated if nothing)
• tol : QTCI tolerance (default 1e-8)
• maxdim : maximum MPO bond dimension after construction (default 15)

Examples

H = get_Hamiltonian("chain_1d", 1.0;    L=10)
H = get_Hamiltonian("square_2d", 1.0;   L=10, Lx=32)

rs = honeycomb_positions(10)
H  = get_Hamiltonian("haldane", (t2=0.2, phi=π/2, M=0.0); L=10, rs=rs)

H  = get_Hamiltonian("custom", (i,j) -> ...; L=10, scale=5.0, geometry=rs)

After construction, add further interaction terms with add_hopping! and add_onsite!.

honeycomb_positions(L; Lx=L÷2) -> Matrix{Float64}

Generate N = 2^L physical honeycomb positions consistent with the quantics row-major encoding n = ix + iy * 2^Lx, bond length = 1.

The lattice is an armchair ribbon: even rows have intra-row bonds (2k, 2k+1) and odd rows have intra-row bonds (2k+1, 2k+2), with all inter-row bonds (iy, ix) ↔ (iy+1, ix).

Returns an N × 2 matrix where row i (1-indexed) is the 2D position of quantics site i-1.

square_positions(L; Lx=L÷2) -> Matrix{Float64}

Physical positions for the 2^L-site square lattice in quantics row-major encoding n = ix + iy·2^Lx. Site i (1-indexed) maps to (ix, iy).

triangular_bravais_positions(L; Lx=L÷2) -> Matrix{Float64}

Physical positions for the 2^L-site Bravais triangular lattice in quantics row-major encoding n = ix + iy·2^Lx, bond length = 1. Bravais vectors a1=(1,0), a2=(1/2,√3/2): x = ix + iy/2, y = iy·√3/2.

triangular_positions(L; Lx=L÷2) -> Matrix{Float64}

Physical positions for the 2^L-site triangular lattice in quantics row-major encoding n = ix + iy·2^Lx, bond length = 1. Odd rows are offset by 0.5 in x: x = ix + 0.5·(iy % 2), y = iy·√3/2.

truncate!(H::TBHamiltonian; cutoff=1e-10, maxdim=nothing) -> H

Truncate the Hamiltonian MPO in-place using ITensorMPS.truncate!. Invalidates all caches (Chebyshev list, density matrix, scale/center).

Useful after a series of add_hopping! / add_onsite! calls that may have inflated the bond dimension.