diff --git a/Project.toml b/Project.toml
index b753c53..43121bb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -6,7 +6,6 @@ authors = ["Tiago Antao, Anouar Moustaj, Yitao Sun"]
[deps]
Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -20,7 +19,6 @@ TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
[compat]
Arpack = "0.5"
FFTW = "1"
-HDF5 = "0.17"
ITensorMPS = "0.2"
ITensors = "0.6"
NDTensors = "0.3"
diff --git a/README.md b/README.md
index aec5344..0c6299f 100644
--- a/README.md
+++ b/README.md
@@ -8,7 +8,7 @@
[](https://github.com/TensorBinding/TensorBinding.jl/actions/workflows/CI.yml?query=branch%3Amain)
[](https://codecov.io/gh/TensorBinding/TensorBinding.jl)
-**TensorBinding.jl** is a Julia package for constructing and studying tight-binding Hamiltonians as **Matrix Product Operators (MPOs)** in the *quantics binary* (QTT) representation. A system of *N = 2L* sites is encoded in *L* qubit sites, keeping bond dimensions small (typically ≤ 10) for physically relevant models. Arbitrary hopping matrices are compressed automatically via **Quantics Tensor Cross Interpolation (QTCI)**.
+**TensorBinding.jl** is a Julia package for constructing and studying tight-binding Hamiltonians as **Matrix Product Operators (MPOs)** in the *quantics binary* (QTT) representation. Ordinarily, *N = 2L* sites are encoded in *L* qubit sites; projected position spaces can embed a different physical site count in the same register. This keeps bond dimensions small (typically ≤ 10) for physically relevant models. Arbitrary hopping matrices are compressed automatically via **Quantics Tensor Cross Interpolation (QTCI)**.
---
@@ -33,6 +33,9 @@ See the [`examples/`](examples/) folder for notebooks covering the main workflow
**Hamiltonian construction**
- 1D: nearest-neighbour chain, SSH (uniform and sublattice-explicit), Aubry–André–Harper quasicrystal, uniform with on-site potential
+- Fibonacci onsite and hopping quasicrystals in a projected Zeckendorf basis, with open or physical-periodic boundaries and conumber ordering
+- Metallic-mean quasicrystals (`A → AᵐB`, `B → A`: silver mean, bronze mean, …) in a projected `(m+1)`-ary numeration basis on Qudit registers; `m = 1` reproduces the Fibonacci chain
+- k-bonacci quasicrystals (`aᵢ → a₁aᵢ₊₁`, `aₖ → a₁`: Tribonacci, Tetranacci, …) in a projected binary basis with no `k` consecutive ones; `k = 2` reproduces the Fibonacci chain
- 2D: square, triangular, honeycomb, kagomé, Lieb, and dice lattices — including sublattice-explicit models with an explicit unit-cell index
- Generic *n*th-nearest-neighbour hopping on any 2D geometry (`add_hopping_2D!`): uniform, direction-dependent, site-dependent, or fully position+direction-dependent amplitude functions
- Arbitrary hopping matrix `f(i,j)` compressed via QTCI (`hopping2MPO`)
@@ -47,6 +50,7 @@ See the [`examples/`](examples/) folder for notebooks covering the main workflow
- Chebyshev expansion of spectral functions, LDOS, Green's functions, and density matrices
- Kernels: Jackson (default), Lorentz, Fejér, Dirichlet, HODC
- Three complementary modes: MPO (full operator), diagonal/online (memory-efficient LDOS), MPS (reference-state propagation)
+- Exact trace DOS from an online three-MPO recursion, including projected position spaces
- Band structure *A(k,ω)* via QFT conjugation (`get_bands`); supports spin, BdG, layer, and sublattice projections via `aux_proj`
- Density matrix purification: McWeeny (cubic convergence) and SP2 (electron-number controlled)
diff --git a/docs/dev/REORGANISATION_TODO.md b/docs/dev/REORGANISATION_TODO.md
new file mode 100644
index 0000000..a7d0096
--- /dev/null
+++ b/docs/dev/REORGANISATION_TODO.md
@@ -0,0 +1,192 @@
+# TensorBinding reorganisation — todo
+
+Outcome of the code-organisation review of 2026-09-23 (four area sweeps over
+`lattice/` + `core/Hamiltonian.jl`, `solvers/`, `physics/`, `gpu/`, plus repo-wide
+metrics). Nothing here is implemented yet. Tiers are ordered so that each one can be
+done and merged on its own with the test suite as the guard; Tier 1 changes no
+behaviour, Tier 2 changes internals only, Tier 3 is user-visible.
+
+Line numbers refer to the working tree on that date and will drift.
+
+## Bugs (fix first, independently of the reorganisation)
+
+- [ ] `core/TBSystem.jl` `_build_haldane` calls `haldane_hoppingf`, which is not defined
+ anywhere; `get_Hamiltonian("haldane", …)` throws. Restore the function or drop the model.
+- [ ] `physics/RPA_tk.jl:935` `get_bubble_mpo_haydock` calls `_build_heff` with 3 arguments;
+ the definition (l.454) takes 4.
+- [ ] `physics/RPA_tk.jl` ~1065 and ~1249: `get_rpa_susceptibility_wynn` and
+ `get_magnon_susceptibility_wynn` assign `nq` inside `if chi_partial === nothing` inside
+ the frequency loop, so the second frequency hits an undefined variable. Hoist `nq`.
+- [ ] `physics/Topology_tk.jl` `get_C` accepts `Lambda` (ASCII alias) and never uses it; the
+ GPU twin honours it.
+- [ ] `core/TBSystem.jl` 17-argument `TBHamiltonian` compatibility constructor silently drops
+ `Lx`, `interaction_mpo`, `fock_mpo` and `position_space`; SCF (l.364, 485, 1136),
+ RPA (l.1141) and NH (l.99) copy Hamiltonians through it. Replace with a keyword copy
+ constructor (see Tier 2) and delete the positional ones.
+- [ ] `solvers/KPM_tk.jl` `get_density_quantics` uses an undefined global `sites`. Delete.
+- [ ] `solvers/Timeev_tk.jl` `compare_propagator_and_tdvp_heatmaps` calls `heatmap`/`plot`/
+ `display` although Plots is not a dependency. Move to `examples/`.
+- [ ] `gpu/GPU_tk.jl:871` second `_sample_state_amplitudes_gpu` call drops `pointavg`.
+- [ ] `gpu/GPU_tk.jl:267` `_onehot_gpu` only accepts `T<:Complex`; the advertised
+ `type=Float32/Float64` paths fail.
+- [ ] `physics/Purification_tk.jl:20` header example uses `method=:KPM`; code accepts `:kpm`.
+- [ ] `README.md:22` claims CUDA is an installed dependency; `Project.toml` has none.
+- [ ] `core/Utils.jl` `_exciton_block_groups` is reachable only through a branch that
+ `get_exciton_ldos_spatial_gpu` rejects earlier (`reduce=:block`). Delete both.
+
+## Tier 1 — mechanical, no behaviour change
+
+### Split the three grab-bag files
+- [ ] `solvers/KPM_tk.jl` (2067 lines) → `solvers/kpm/kernels.jl` (`_kpm_kernel`,
+ `_dos_weight_matrix`, HODC helpers, `_kpm_weight_matrix` from QFT), `recursion.jl`
+ (`KPM_Tn`, `KPM_Tn_mps`, `_run_kpm_mps!`), `cached.jl` (`get_ldos`, `get_ldos_spectrum`,
+ `*_from_Tn`, `*_from_mun`, Green's functions), `ldos.jl` (`get_ldos_online`,
+ `get_ldos_spatial`, split into `_ldos_spatial_mps`/`_ldos_spatial_mpo`), `dos.jl`
+ (`get_dos_stochastic`, `get_dos_trace`), `exciton.jl` (l.1648–1984).
+- [ ] `physics/RPA_tk.jl` (2148 lines) → `physics/rpa/Bubble.jl`, `Cheb2D.jl`, `Dyson.jl`;
+ MPO kron/interleave plumbing (l.10–287) → `core/Utils.jl`; Haydock recursion →
+ `solvers/Krylov_tk.jl`; `get_spect_k` → QFT conjugation file; delete l.288–377.
+- [ ] `physics/QFT_tk.jl` (1643 lines) → `Conjugation.jl` (l.96–205), `Bands.jl`
+ (l.638–948, 1223–1370), `KPath.jl` (l.426–635); exciton spectra (l.206–281, 949–1220)
+ → exciton folder; aux projection (l.1373–1505) → `core/AuxDOF.jl`.
+- [ ] `physics/NH_tk.jl` → `NH_model.jl` (struct, `hermitize`, `add_nh_*`) and `NH_KPM.jl`.
+- [ ] `lattice/2Dlattice_tk.jl` (1615 lines) → `Masks2D.jl`, `Hopping2D.jl`, `Presets.jl`
+ (QTCI `H*` builders incl. the 1D `HUniform`/`HSSH`/`HAAH`), `Sublattice.jl`
+ (kagome/lieb/dice/honeycomb), `Geometry.jl`; `MODEL_REGISTRY`/`build_hamiltonian` →
+ `core/ModelRegistry.jl`.
+- [ ] `gpu/GPU_tk.jl` (3647 lines) → `device.jl`, `primitives.jl`, `kpm.jl`, `bands.jl`,
+ `topology.jl`, `purification.jl`, `scf.jl`, `exciton.jl`, `nh.jl`, `timeev.jl`;
+ the conductivity-only Tucker/QFT/Hadamard block (~300 lines) → its example.
+
+### Move misplaced helpers next to their callers
+- [ ] One `core/AuxDOF.jl` owning spin/Nambu indices and op tables, `prepend_spin`/`prepend_nambu`,
+ Symbol overloads of `prepend_op`/`postpend_op` (from `Supercond_tk.jl`), `project_aux`,
+ `aux_site`, `_autoenable_proj` (from QFT), `_aux_setup`, `_ldos_make_psi0` (from KPM),
+ and the four `add_spin!`/`add_zeeman!`/`add_superconductivity!`/`add_soc!` mutators
+ (from TBSystem). Include it right after TBSystem.
+- [ ] `_estimate_spectral_bounds` → `solvers/DMRG_tk.jl`; include DMRG before KPM.
+- [ ] `_eval_diag_mps` → `core/Utils.jl` beside `eval_mps`; `mpsexciton` → Utils beside the
+ other product-state builders.
+- [ ] `qtt_mpo`, `compose_power`, `_row_break/_row_select/_col_select/_row_checker_mpo`,
+ `_site_projector_mpo`, `sigma_d/sigma_u` ops, layer prepend helpers → `core/Utils.jl`
+ (or `lattice/Masks2D.jl` for the masks).
+- [ ] BdG/pairing builders in `SCF_tk.jl` (l.298–528) → AuxDOF / Supercond.
+- [ ] `_project_spin_sector` (RPA) → AuxDOF as `project_sector(H, :spin, σ)`.
+- [ ] All geometry (`*_positions`, `_*_geometry`, `lattice_positions`, `_resolve_2d_geometry`,
+ junction geometry, `geometry_uc` closures) → `lattice/Geometry.jl` with one `(Lx, Ly)`
+ signature.
+- [ ] `_reconstruct_ldos_moment_columns` (GPU) → `solvers/kpm/kernels.jl`; move its test out of
+ `test/gpu_mps_ldos.jl`.
+
+### Delete dead and legacy code
+- [ ] Confirmed unreferenced everywhere (incl. notebooks and generated docs):
+ `build_cyclic_shift_mpo`, `_geom_n_sub`, `_nsublat`, `nsitelegs`, `_tb_spatial_groups_gpu`,
+ `get_nh_state_trajectory_gpu`, the `Delta_*` one-liners in SCF.
+- [ ] Unreferenced in src/test/tracked examples: `projop_2DSL`, `projop_1DSL`, `sample_diag`,
+ `project_spin`, `get_density_quantics`, `_get_exciton_ldos_cached` + exciton
+ `KPM_Tn(H, N, X)`, `ldos_exc_KPM_Tn`, `get_mus_raw`, `compute_dos_ldos_hodc`,
+ `kinetic_1d_nn_custom`, `qtci_matrix_to_MPO`, `quasicrystal_modulation_30deg`,
+ `circular_mod`, `interchain_hopping_*` (2nd_plus/minus, triangle, honeycomb) with their
+ skeleton/template helpers, `postpend_layer_projector/hopping`, `sdf_interval`,
+ `mps_kron`, `merge_mps_to_mpo`, `convert_mpo`, `_swap_mpo`, `apply_interleave_swaps`,
+ `get_Tnlists`, `get_bublle_expanded_from_Tn`, `build_bubble_mpo`,
+ `get_bubble_mpo_haydock`, `hopping_mpo_exciton`, `get_valley_projectors`,
+ `fock_exchange_builder`, `initial_guess_trivial_*_1D`, `nh_imag_onsite_mpo`,
+ `add_nh_imag_onsite!`, `add_nh_loss!`, `nh_reconstruct_spectral_mpo`,
+ `nh_spectral_function_allsite_mpo`, `spin_hamiltonian`, `bdg_hamiltonian` (re-inlined
+ in TBSystem), `_onehot_gpu_f32`, `nh_spectrum_grid_gpu`. Check each once more before
+ deleting; `examples/nontracked/APSOS/Modified_GPU_funcs.jl` carries forks of some.
+- [ ] Commented-out legacy: `QFT_tk.jl:1511–1643` (old `get_bands`, `get_spect_k*`),
+ `Purification_tk.jl:95–96`, unreachable code after early `return` in
+ `2Dlattice_tk.jl` (`generate_kin_u/d` l.33–63, six kinetic builders l.388–543).
+- [ ] Six positional "backward-compatible" `TBHamiltonian` constructors (TBSystem l.98–116,
+ 190–214) once Tier 2 keyword constructor exists.
+- [ ] Unconditional `println` in library code (~70 in src): `Hamiltonian.jl` 85–123,
+ `KPM_tk.jl` 14/30/31, `QFT_tk.jl` 1453–1470, `Topology_tk.jl` 499–539,
+ `TBSystem.jl` 1175, RPA legacy pipeline; switch to `@info … maxlog=1` or `verbose` gates.
+
+### Make the structure legible
+- [ ] Explicit `export` list (today only ITensors names are exported) so public vs private is visible.
+- [ ] One banner style (`# ====` vs `# ───` vs none); numbered sections that match contents
+ (2Dlattice runs 8, 8b, 8c, 8d, 8f; SCF header lists 8 sections, file has 9).
+- [ ] Rewrite the load-order comment in `TensorBinding.jl` as a real dependency graph; fix the
+ include order where a solver depends on a physics file (KPM ↔ QFT, TBSystem → Supercond,
+ Krylov → RPA, Bilayer → Twisted, SCF/RPA/Topology → Purification).
+- [ ] File names: drop the `_tk` suffix; rename `2Dlattice_tk.jl`; fix header comments that
+ cite files that do not exist (`utils.jl`, `2D_lattice.jl`, `twoparticle_tk.jl`, `krylov_tk.jl`).
+- [ ] Re-save `2Dlattice_tk.jl` as UTF-8 and restore the mojibake symbols (√, ·, ≠ appear as
+ `-`/`_`, e.g. `b=(1+-)/2` for the golden ratio).
+- [ ] Docstrings vs signatures: `get_ldos_spatial` omits 9 kwargs; `get_ldos_from_mun` omits
+ `eta`/`m_order`; Bilayer/Twisted claim `(MPO, sites)` returns but return `TBHamiltonian`;
+ Flake/TBSystem examples pass `Lx=16`/`32` where `Lx` is a qubit count; `get_Hamiltonian`
+ table lists 8 of 21 names; QFT table of contents (l.76–92) wrong in five places;
+ Topology header lists `berry_curvature_integrand`, which does not exist.
+- [ ] Tests: lattice builders, RPA, SCF, NH, Topology have no tests; add smoke tests before
+ splitting so the moves are guarded.
+
+## Tier 2 — shared kernels (internal behaviour only)
+
+- [ ] `_scaled_hamiltonian(H; cutoff)` = `(1/scale)·(H − center·physical_projector(H))`,
+ replacing ~20 inline copies (some use `MPO(sites,"Id")` and mishandle projected spaces:
+ `KPM_tk.jl` 1799, 1917, 1675; `QPI_tk.jl` 155).
+- [ ] `chebyshev_foreach(f!, H̃, T₀; maxdim, cutoff)` working for MPO and MPS on any device,
+ replacing ~22 hand-written three-term loops (6 KPM, 14 GPU, QFT, QPI) and 5 NH partial
+ recurrences; one truncation policy.
+- [ ] `_kpm_energy_grid(H, ωs; kernel, …) -> (ω_r, W, denom, valid)` replacing 14 copies of the
+ rescale/weights/valid block and 7 hand-written `π²·N·√(1−ω²)` normalisations.
+- [ ] `_chebyshev_sum(Tn, coeffs; …)` replacing 6 weighted-sum copies; HODC variants become a
+ coefficient choice.
+- [ ] One Jackson kernel (`_kpm_kernel`) with a `normalize` keyword; delete `_jackson_kernel`
+ (RPA) and `nh_jackson_weights` (NH).
+- [ ] `AuxProjection` struct (or `aux...` kwargs forwarded to `_aux_setup`) replacing the
+ 8-keyword block copied into ~10 signatures; one `_project_aux_sectors` replacing the
+ nambu→spin→layer→sublattice chain written 4× (KPM, QFT, GPU ×2) and the 4 sector
+ projectors (`project_aux`, `_project_aux_block`, `_project_spin_sector`, `contract_nh_block`).
+- [ ] `probe_state(H, x, σ…)` replacing the psi0 selection duplicated 3× in KPM.
+- [ ] Keyword `TBHamiltonian(; L, N, sites, mpo, …)` plus `similar(H; mpo=, sites=, …)` copy
+ constructor; delete the six positional overloads.
+- [ ] One model registry entry per model (builder → `TBHamiltonian`, dim, params, geometry,
+ scale) replacing `get_Hamiltonian`'s if-chain + `build_hamiltonian` + `_build_preset` +
+ `_build_sublattice` + `_preset_geometry` + `_estimate_scale`; `_param(params, :t, default)`
+ replacing the parsing ternaries; remove drifted `kw_defaults` from the registry.
+- [ ] One `masked_shift_hopping(Lx, Ly, sites, hop, q; src_mask)` replacing six near-identical
+ 2D kinetic builders; retire `generate_kin_u/d` in favour of `shift_mpo`.
+- [ ] `_sublattice_bond` + `_sublattice_setup` replacing ~12 repeated bond blocks in
+ kagome/lieb/honeycomb/dice; `_basis_positions` replacing 4 identical position loops;
+ `sum_mpos(terms; cutoff)`.
+- [ ] `get_density` as the only projector dispatcher (delete `_get_density_matrix` in RPA and
+ `_get_projector` in Topology); `_purified_pair` for the ρ± blocks in Purification.
+- [ ] RPA: `_cheb2d_setup` + `_tucker_bases` (5 copied prologues, 2 Tucker blocks); one Wynn
+ driver (3 copies); magnon functions as `mode=:magnetic`.
+- [ ] Timeev: `_rk4_step(rhs, …)` (2 copies), one `evolve_rk4_dm_*`, one trajectory loop;
+ remove the double normalisation after `tdvp(normalize=true)`.
+- [ ] GPU: thin wrappers over CPU kernels with a `to_device` hook (stochastic DOS, McWeeny/SP2,
+ Chern operator assembly, NH kernels, `_eval_block_mps`, `extract_diagonal_to_mps`,
+ `mps_to_diagonal_mpo`, `density_profile_from_dm`); one `_to_gpu(x, T)`; one
+ `_resolve_gpu_type` with a single warning threshold; `_gpu_log`.
+- [ ] Decide the one remaining sampling divergence: `get_ldos_spatial_mps_gpu` automatic plan
+ (balanced `fld` bins, `unique(round.(range))` samples) vs `spatial_sampling_plan` 1D branch.
+
+## Tier 3 — API consistency (user-visible)
+
+- [ ] `Ncheb` everywhere, positional (today `N`, `Ncheb`, `Nchebychev`, NH `n` meaning 2n).
+- [ ] `cutoff` for SVD truncation; `tci_tol` / `krylov_tol` / `scf_tol` for the others
+ (`tol` currently means four things).
+- [ ] `boundary` only (drop `bc`, `cyclic` aliases); `maxdim` defaults from one
+ `const KPM_DEFAULTS`; document the loose `tol=1e-8, maxdim=15` that `get_Hamiltonian`
+ hands to every builder.
+- [ ] `dtype` only (drop `type`); one `verbose::Int` level (drop `printinfo`).
+- [ ] Method symbols in one case (`:kpm`, not `:KPM`); `fermi` vs `ϵF`; `Λ` vs `Lambda`;
+ `omega` vs `ω_phys_vals`; exciton momenta `Q_*` only, one indexing convention.
+- [ ] Return NamedTuples instead of kwarg-dependent shapes (`get_bands` Matrix/NamedTuple,
+ `get_ldos_spatial_mps_gpu` four shapes, `get_ldos` MPS/MPO/Real/nothing, `thouless_pump`,
+ `nh_spectrum_grid`, the four SCF drivers); an `SCFResult` struct.
+- [ ] Split `mode` into `output=:operator|:diagonal` and `algorithm=:mpo|:mps`.
+- [ ] Naming: `chern_marker`/`winding_marker` (keep `get_C`/`get_W` as deprecated aliases),
+ `_hamiltonian` everywhere, lowercase `_mpo` (`hopping2MPO` → `hopping_mpo`),
+ `exciton_mpo` for `Exciton_Hamiltonian`, fix `get_bublle_expanded_from_Tn`.
+- [ ] Replace hidden mutable caches (`_tn_cache`, `_tn_mps_cache`, `_density_cache`,
+ `_ensure_scale!` side effects, solvers mutating user Hamiltonians) with an explicit
+ `KPMExpansion` object passed to the reconstruction functions.
+- [ ] CUDA as a package extension (`[weakdeps] CUDA`, `ext/TensorBindingCUDAExt/`), replacing
+ the `Base.loaded_modules` UUID lookup; fix the README dependency statement.
diff --git a/docs/make.jl b/docs/make.jl
index 4db1488..45389cd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -18,6 +18,7 @@ makedocs(;
"Home" => "index.md",
"API Reference" => [
"Core" => "api/core.md",
+ "Position Spaces" => "api/position_spaces.md",
"Lattice" => "api/lattice.md",
"Solvers" => "api/solvers.md",
"Physics" => "api/physics.md",
diff --git a/docs/src/api/position_spaces.md b/docs/src/api/position_spaces.md
new file mode 100644
index 0000000..6be7843
--- /dev/null
+++ b/docs/src/api/position_spaces.md
@@ -0,0 +1,171 @@
+```@meta
+CurrentModule = TensorBinding
+```
+
+# Position Spaces
+
+## Fibonacci quasicrystal
+
+Fibonacci chains use a projected Zeckendorf position space: `H.N = F_(L+2)`
+physical positions are embedded in an ambient `2^L` qubit register. The physical
+identity is therefore `physical_projector(H)`, not the ambient identity.
+
+```julia
+H = TensorBinding.fibonacci_hamiltonian(
+ 8; A=1.0, B=2.0,
+ model=:hopping,
+ boundary=:periodic,
+)
+
+# Equivalent generic constructor
+H = TensorBinding.get_Hamiltonian(
+ "fibonacci", (A=1.0, B=2.0);
+ L=8, model=:hopping, boundary=:periodic,
+)
+```
+
+`model=:onsite` uses `A` and `B` as onsite energies with uniform hopping `t`.
+`model=:hopping` uses them as bond amplitudes with uniform `onsite`. Periodic
+boundaries close the physical Fibonacci approximant; they do not wrap at the
+last ambient binary state.
+
+CPU KPM construction, cached/online LDOS, spatial LDOS, stochastic DOS,
+deterministic trace DOS, and KPM density construction are projector-aware.
+Binary-only modifiers, QFT/bands, GPU, purification, topology, SCF, exciton,
+and non-Hermitian APIs currently reject projected position spaces explicitly.
+
+```julia
+energies = range(-4, 4; length=400)
+
+dos = TensorBinding.get_dos_trace(H, 200, energies)
+ldos = TensorBinding.get_ldos_spatial(
+ H, 200, energies;
+ ordering=:conumber,
+ conumber_orientation=:standard,
+ conumber_centered=true,
+ conumber_alignment=:atomic,
+)
+conumbers = TensorBinding.site_axis(H; ordering=:conumber)
+```
+
+The two multipliers are reflections because
+`F_(n-1) = F_n - F_(n-2)` modulo `F_n`; changing orientation cannot repair a
+wrong cyclic phase. The raw modular formula alone does not choose where the periodic perpendicular-
+space interval is cut. The default `alignment=:atomic` chooses that cut to give
+three contiguous blocks of sizes `F_L | F_(L-1) | F_L`: molecular, atomic
+(`AA`), molecular. `centered=true` then labels this ordered axis around zero;
+it does not rotate it again. Set `alignment=:raw` to inspect unshifted residues,
+or use `orientation=:reversed` for the reflected perpendicular-space direction.
+
+Pointwise helpers avoid allocating a full permutation for very large systems:
+
+```julia
+c = TensorBinding.fibonacci_conumber(43, site)
+site_again = TensorBinding.fibonacci_site_from_conumber(43, c)
+kind = TensorBinding.fibonacci_site_environment(43, site)
+depth = TensorBinding.fibonacci_atomic_depth(43, site)
+```
+
+Successive atomic deflations map `L -> L-3`. To zoom without losing the induced
+phason/origin, retain the original conumber coordinates and slice the nested
+window instead of assigning the selected sites fresh indices `1:F_(L'+2)`:
+
+```julia
+zoom = TensorBinding.fibonacci_rg_partition(43; depth=13)
+@assert zoom.effective_L == 4
+@assert zoom.window_count == 8
+@assert (zoom.molecular_count, zoom.atomic_count, zoom.molecular_count) == (3, 2, 3)
+
+# Probe those original physical sites in their inherited conumber order.
+zoom_sites = [
+ TensorBinding.fibonacci_site_from_conumber(43, c; centered=false)
+ for c in zoom.window_ranks
+]
+ldos_zoom = TensorBinding.get_ldos_spatial(
+ H, 4000, energies;
+ ordering=:physical,
+ x_groups=[[site] for site in zoom_sites],
+)
+# Plot ldos_zoom against zoom.window_axis; do not conumber it a second time.
+```
+
+A periodic hopping ring with odd `H.N` is not bipartite and is therefore not
+required to have exact `E -> -E` chiral symmetry.
+
+```@autodocs
+Modules = [TensorBinding]
+Pages = ["position_spaces/Fibonacci.jl"]
+```
+
+## Metallic-mean quasicrystals
+
+The metallic-mean chain with parameter `m` is the fixed point of `A -> A^m B`,
+`B -> A` (`m = 1` Fibonacci, `m = 2` silver mean, `m = 3` bronze mean). Sites
+are labelled in the numeration system with basis `q_0 = 1`, `q_1 = m + 1`,
+`q_(l+1) = m q_l + q_(l-1)`, digits in `0:m`, and the rule that a digit `m` must
+be followed by `0`. `L` digits enumerate `H.N = q_L` physical sites inside an
+ambient `(m+1)^L` register of `Qudit` sites of dimension `m + 1`. The letter at
+site `n` is `B` exactly when the least significant digit of `n` is `m`, so the
+word and the validity indicator are both bond-dimension-2 automaton MPS and the
+Hamiltonian MPO `P (V + T K + h.c.) P` is exact at any `L`.
+
+```julia
+H = TensorBinding.metallic_mean_hamiltonian(
+ 2, 8; A=1.0, B=2.0, model=:hopping, boundary=:periodic,
+)
+
+# Equivalent generic constructor (m is required)
+H = TensorBinding.get_Hamiltonian(
+ "metallic_mean", (A=1.0, B=2.0); L=8, m=2,
+)
+```
+
+The projector-aware CPU KPM entry points work exactly as for Fibonacci
+(`KPM_Tn`, `get_ldos_online`, `get_ldos_spatial` with `ordering=:physical`,
+`get_dos_stochastic`, `get_dos_trace`), as does `get_ldos_spatial_mps_gpu`.
+Conumbering and the inherited-conumber sampling plans are currently
+Fibonacci-only, so `ordering=:conumber` throws for metallic means.
+`m = 1` reproduces the Fibonacci chain of `fibonacci_hamiltonian` on
+dimension-2 `Qudit` sites.
+
+```@autodocs
+Modules = [TensorBinding]
+Pages = ["position_spaces/MetallicMean.jl"]
+```
+
+## k-bonacci quasicrystals
+
+The k-bonacci chain on the alphabet `a_1, …, a_k` (written `A, B, C, …`) is the
+fixed point of `a_i -> a_1 a_(i+1)` for `i < k` and `a_k -> a_1` (`k = 2`
+Fibonacci, `k = 3` Tribonacci `A -> AB, B -> AC, C -> A`, `k = 4` Tetranacci).
+Sites are labelled by binary strings with no `k` consecutive ones, read with the
+weights `w_l = 2^l` for `l < k` and `w_l = w_(l-1) + … + w_(l-k)` otherwise
+(Zeckendorf for `k = 2`, the Tribonacci numbers `T_(l+3)` for `k = 3`). `L`
+digits enumerate `H.N = w_L` physical sites inside the ambient `2^L` `Qubit`
+register, exactly like the Fibonacci chain. The letter at site `n` is `a_(r+1)`
+where `r` is the number of trailing ones of `n`, so the word and the validity
+indicator are `k`-state automaton MPS (bond dimension `k`) and the Hamiltonian
+MPO `P (V + T K + h.c.) P` is exact at any `L`. The decrement `K` clears the
+least significant one and rewrites the tail with the pattern `1^(k-1) 0`.
+
+```julia
+# Tribonacci hopping chain with t_A/t_B = t_B/t_C = 0.8 and t_C = 1
+H = TensorBinding.kbonacci_hamiltonian(
+ 3, 10; values=(0.64, 0.8, 1.0), model=:hopping, boundary=:periodic,
+)
+
+# Equivalent generic constructor (k is required); letter keys or values=(…)
+H = TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8, C=1.0); L=10, k=3,
+)
+```
+
+The projector-aware CPU KPM entry points and `get_ldos_spatial_mps_gpu` work as
+for Fibonacci with `ordering=:physical`; conumbering is Fibonacci-only, so
+`ordering=:conumber` throws. `k = 2` reproduces `fibonacci_hamiltonian` exactly,
+on the same `Qubit` sites.
+
+```@autodocs
+Modules = [TensorBinding]
+Pages = ["position_spaces/KBonacci.jl"]
+```
diff --git a/examples/basics/getting_started.ipynb b/examples/basics/getting_started.ipynb
index 59290c2..6e3b31f 100644
--- a/examples/basics/getting_started.ipynb
+++ b/examples/basics/getting_started.ipynb
@@ -25,7 +25,7 @@
},
{
"cell_type": "code",
- "execution_count": 3,
+ "execution_count": 1,
"id": "gs_imports",
"metadata": {},
"outputs": [],
diff --git a/examples/misc/nonhermitian_loss_chain.ipynb b/examples/misc/nonhermitian_loss_chain.ipynb
index 512fae5..50fd23f 100644
--- a/examples/misc/nonhermitian_loss_chain.ipynb
+++ b/examples/misc/nonhermitian_loss_chain.ipynb
@@ -16,7 +16,7 @@
},
{
"cell_type": "code",
- "execution_count": 2,
+ "execution_count": 1,
"id": "794ec812",
"metadata": {},
"outputs": [],
diff --git a/src/TensorBinding.jl b/src/TensorBinding.jl
index bfb372f..4805df8 100644
--- a/src/TensorBinding.jl
+++ b/src/TensorBinding.jl
@@ -23,11 +23,22 @@ export MPO, MPS, OpSum, expect, inner, siteinds
# Load order matters:
# core/Utils.jl — binary/index helpers, diagonal MPO construction,
-# shift/Hadamard operators (no deps)
+# shift/Hadamard operators, and every sampling plan
+# (spatial, k-space, Fibonacci) shared by the CPU and
+# GPU solvers (no deps)
# core/Hamiltonian.jl — 1D/2D kinetic operator and QTCI MPO builders,
# preset model Hamiltonians (uses Utils)
-# core/TBSystem.jl — TBHamiltonian struct, get_Hamiltonian, add_*!
-# mutators (uses Utils, Hamiltonian)
+# core/TBSystem.jl — position-space policy types, TBHamiltonian struct,
+# get_Hamiltonian, add_*! mutators and the
+# position-space interface (uses Utils, Hamiltonian)
+# position_spaces/Fibonacci.jl — projected Fibonacci space, automata,
+# constructors, and conumbering (uses TBSystem)
+# position_spaces/MetallicMean.jl — projected metallic-mean spaces (A -> A^m B,
+# B -> A) on (m+1)-dimensional Qudit
+# registers (uses TBSystem, Utils)
+# position_spaces/KBonacci.jl — projected k-bonacci spaces (Tribonacci,
+# Tetranacci, …) on the binary register with
+# no k consecutive ones (uses TBSystem, Fibonacci)
# lattice/2Dlattice_tk.jl — 2D shift operators, lattice hoppings, geometry
# positions (uses Utils, Hamiltonian, TBSystem)
# lattice/NNNeighbor_tk.jl — generic nth-neighbor hopping accumulator
@@ -67,6 +78,9 @@ export MPO, MPS, OpSum, expect, inner, siteinds
include("core/Utils.jl")
include("core/Hamiltonian.jl")
include("core/TBSystem.jl")
+include("position_spaces/Fibonacci.jl")
+include("position_spaces/MetallicMean.jl")
+include("position_spaces/KBonacci.jl")
include("lattice/2Dlattice_tk.jl")
include("lattice/NNNeighbor_tk.jl")
include("lattice/Flake_tk.jl")
diff --git a/src/core/TBSystem.jl b/src/core/TBSystem.jl
index d74f027..d015333 100644
--- a/src/core/TBSystem.jl
+++ b/src/core/TBSystem.jl
@@ -5,6 +5,24 @@
# and the density matrix. All observable methods (get_DoS, get_density,
# get_Chern, get_bands …) dispatch on this struct.
+# ============================================================
+# Position-space policy types
+# ============================================================
+
+"""
+ AbstractPositionSpace
+
+Policy object describing how physical positions are embedded in the tensor-product
+register. `BinaryPositionSpace` is the ordinary `N = 2^L` quantics basis. Other
+position spaces (see `position_spaces/`) specialize `physical_projector`,
+`physical_site_state`, `site_axis`, and `site_permutation` after `TBHamiltonian`
+is defined below.
+"""
+abstract type AbstractPositionSpace end
+
+"""Ordinary binary position register containing all `2^L` basis states."""
+struct BinaryPositionSpace <: AbstractPositionSpace end
+
# ============================================================
# TBHamiltonian struct
# ============================================================
@@ -18,10 +36,11 @@ Fields
------
**Core**
- `L` : number of position qubit sites (log₂ of the physical system size)
-- `N` : number of physical sites / unit cells (2^L)
+- `N` : number of physical sites / unit cells (`2^L` for the binary basis)
- `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
+- `position_space`: policy describing the physical basis inside the tensor register
**KPM spectral bounds**
- `scale` : energy half-bandwidth; `H/scale` has spectrum in `[-1, 1]`.
@@ -73,8 +92,19 @@ mutable struct TBHamiltonian
interaction_mpo :: Union{Nothing, MPO}
fock_mpo :: Union{Nothing, MPO}
Lx :: Union{Nothing, Int} # x-qubit count for 2D (Ly = L - Lx); nothing for 1D
+ position_space :: AbstractPositionSpace
end
+# Backward-compatible full constructor (pre-position_space callers).
+TBHamiltonian(L, N, sites, mpo, geometry, geometry_uc, scale, center,
+ spin_s, nambu_s, layer_s, sublattice_s, aux_side,
+ _tn_cache, _tn_mps_cache, _tn_Ncheb, _density_cache,
+ interaction_mpo, fock_mpo, Lx) =
+ TBHamiltonian(L, N, sites, mpo, geometry, geometry_uc, scale, center,
+ spin_s, nambu_s, layer_s, sublattice_s, aux_side,
+ _tn_cache, _tn_mps_cache, _tn_Ncheb, _density_cache,
+ interaction_mpo, fock_mpo, Lx, BinaryPositionSpace())
+
# Backward-compatible 17-arg constructor (pre-interaction_mpo/pre-fock_mpo/pre-Lx callers);
# appends nothing, nothing, nothing.
TBHamiltonian(L, N, sites, mpo, geometry, geometry_uc, scale, center,
@@ -85,6 +115,78 @@ TBHamiltonian(L, N, sites, mpo, geometry, geometry_uc, scale, center,
_tn_cache, _tn_mps_cache, _tn_Ncheb, _density_cache,
nothing, nothing, nothing)
+# ============================================================
+# Position-space interface
+# ============================================================
+
+"""
+ ambient_dimension(H) -> Integer
+
+Dimension of the position tensor register before any physical-subspace projection.
+This is `2^H.L` for the quantics encodings supported by TensorBinding. Projected
+position spaces may return a `BigInt` when the ambient register exceeds `Int`.
+"""
+ambient_dimension(H::TBHamiltonian) = ambient_dimension(H.position_space, H)
+ambient_dimension(::BinaryPositionSpace, H::TBHamiltonian) = 2^H.L
+
+"""
+ physical_projector(H) -> MPO
+
+Identity operator on the physical position space. For ordinary binary systems this
+is the full identity; projected encodings specialize this method and return their
+valid-state projector.
+"""
+physical_projector(H::TBHamiltonian) = physical_projector(H.position_space, H)
+physical_projector(::BinaryPositionSpace, H::TBHamiltonian) = MPO(H.sites, "Id")
+
+"""
+ physical_site_state(H, x) -> MPS
+
+Product-state probe for 1-indexed physical position `x`. Auxiliary and two-particle
+spaces use their dedicated probe constructors.
+"""
+physical_site_state(H::TBHamiltonian, x::Integer) =
+ physical_site_state(H.position_space, H, x)
+
+function physical_site_state(::BinaryPositionSpace, H::TBHamiltonian, x::Integer)
+ 1 <= x <= H.N || throw(BoundsError(1:H.N, x))
+ length(H.sites) == H.L ||
+ error("physical_site_state currently requires a position-only TBHamiltonian.")
+ return binary_to_MPS(x - 1, H.L, H.sites)
+end
+
+"""Return the plotting axis for physical positions or an encoding-defined ordering."""
+function site_axis(H::TBHamiltonian; ordering::Symbol=:physical, kwargs...)
+ return site_axis(H.position_space, H; ordering, kwargs...)
+end
+
+function site_axis(::BinaryPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical ||
+ throw(ArgumentError("ordering=:$ordering is not available for BinaryPositionSpace"))
+ return collect(0:(H.N - 1))
+end
+
+"""Return the 1-based physical-site permutation associated with a plotting ordering."""
+function site_permutation(H::TBHamiltonian; ordering::Symbol=:physical, kwargs...)
+ return site_permutation(H.position_space, H; ordering, kwargs...)
+end
+
+function site_permutation(::BinaryPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical ||
+ throw(ArgumentError("ordering=:$ordering is not available for BinaryPositionSpace"))
+ return collect(1:H.N)
+end
+
+_is_binary_position_space(H::TBHamiltonian) = H.position_space isa BinaryPositionSpace
+
+function _require_binary_position_space(H::TBHamiltonian, api::AbstractString)
+ _is_binary_position_space(H) && return nothing
+ throw(ArgumentError("$api is not yet supported for $(typeof(H.position_space)); " *
+ "the first projected-space release supports CPU KPM DOS/LDOS only."))
+end
+
# Backward-compatible 16-arg constructor (pre-geometry_uc callers); inserts geometry_uc=nothing.
TBHamiltonian(L, N, sites, mpo, geometry, scale, center,
spin_s, nambu_s, layer_s, sublattice_s, aux_side,
@@ -146,9 +248,13 @@ Useful after a series of `add_hopping!` / `add_onsite!` calls that may
have inflated the bond dimension.
"""
function truncate!(H::TBHamiltonian; cutoff::Real = 1e-10, maxdim = nothing)
+ old_scale, old_center = H.scale, H.center
kwargs = maxdim === nothing ? (cutoff=cutoff,) : (cutoff=cutoff, maxdim=maxdim)
ITensorMPS.truncate!(H.mpo; kwargs...)
_invalidate_cache!(H)
+ if !_is_binary_position_space(H)
+ H.scale, H.center = old_scale, old_center
+ end
return H
end
@@ -170,6 +276,9 @@ Supported geometry strings
| `"square_2d"` | hopping amplitude `t::Number` | `Lx`, `Ly` (default `L÷2` each) |
| `"haldane"` | `(t2, phi, M)` NamedTuple | `rs` (N×2 Float64 position matrix, required) |
| `"custom"` | hopping function `f(i,j)` | `geometry`, `scale` (required), `type` |
+| `"fibonacci"` | `(A, B[, t, onsite])` NamedTuple | `model=:hopping/:onsite`, `boundary=:periodic/:open` |
+| `"metallic_mean"` | `(A, B[, t, onsite])` NamedTuple | `m` (required; `m=2` silver mean), `model`, `boundary` |
+| `"kbonacci"` | `(A, B, C, ...[, t, onsite])` or `(values=(a_1, ..., a_k)[, t, onsite])` NamedTuple | `k` (required; `k=3` Tribonacci), `model`, `boundary` |
| `"kagome"` | hopping amplitude `t::Number` | `Lx`, `Ly`; 3-atom unit cell, sublattice index postpended |
| `"lieb"` | hopping amplitude `t::Number` | `Lx`, `Ly`; 3-atom unit cell, sublattice index postpended |
@@ -195,6 +304,9 @@ 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)
+Hf = get_Hamiltonian("fibonacci", (A=1.0, B=2.0); L=8, model=:hopping)
+Hs = get_Hamiltonian("metallic_mean", (A=1.0, B=2.0); L=8, m=2) # silver mean
+Ht = get_Hamiltonian("kbonacci", (A=0.64, B=0.8, C=1.0); L=8, k=3) # Tribonacci
```
After construction, add further interaction terms with
@@ -207,6 +319,22 @@ function get_Hamiltonian(geometry::String, params;
maxdim=15,
ref_sites::Union{Nothing,Vector{<:Index}}=nothing,
kwargs...)
+ if geometry == "fibonacci"
+ ref_sites === nothing ||
+ throw(ArgumentError("ref_sites is not supported for FibonacciPositionSpace"))
+ return _build_fibonacci(params, L; scale, tol, maxdim, kwargs...)
+ end
+ if geometry == "metallic_mean"
+ ref_sites === nothing ||
+ throw(ArgumentError("ref_sites is not supported for MetallicMeanPositionSpace"))
+ return _build_metallic_mean(params, L; scale, tol, maxdim, kwargs...)
+ end
+ if geometry == "kbonacci"
+ ref_sites === nothing ||
+ throw(ArgumentError("ref_sites is not supported for KBonacciPositionSpace"))
+ return _build_kbonacci(params, L; scale, tol, maxdim, kwargs...)
+ end
+
sites = siteinds("Qubit", L)
N = 2^L
@@ -241,7 +369,7 @@ function get_Hamiltonian(geometry::String, params;
return _build_preset(geometry, params, L, N, sites; scale, tol, maxdim, ref_sites, kwargs...)
else
- known = ("chain_1d", "haldane", "custom",
+ known = ("chain_1d", "haldane", "custom", "fibonacci", "metallic_mean", "kbonacci",
"uniform", "ssh", "ssh_sublattice", "aah",
"square_2d", "hex_2d", "triangular_2d", "triangular_bravais",
"chern8", "chernhex", "qc2dsquare",
@@ -609,6 +737,7 @@ function add_hopping!(H::TBHamiltonian, f;
sublat = nothing,
sublat_from = nothing,
sublat_to = nothing)
+ _require_binary_position_space(H, "add_hopping!")
if !isnothing(H.Lx)
(!isnothing(sublat) || !isnothing(sublat_from) || !isnothing(sublat_to)) &&
error("add_hopping! sublat keywords are not supported for 2D Hamiltonians; use add_hopping_2D! directly.")
@@ -734,6 +863,7 @@ Invalidates all caches.
"""
function add_onsite!(H::TBHamiltonian, f; layer=nothing, sublat=nothing,
Lx=nothing, tol=1e-8, maxdim=nothing)
+ _require_binary_position_space(H, "add_onsite!")
if H.layer_s !== nothing
(H.spin_s === nothing && H.nambu_s === nothing) ||
error("Layered add_onsite! currently supports layer/position/sublattice Hamiltonians only.")
@@ -873,6 +1003,7 @@ function add_interaction!(H::TBHamiltonian, V;
type::Type = Float64,
tol::Real = 1e-8,
kwargs...)
+ _require_binary_position_space(H, "add_interaction!")
pos_s = _pos_sites(H)
mpo = if V isa MPO
V
@@ -911,6 +1042,7 @@ Invalidates all caches.
"""
function add_spin!(H::TBHamiltonian; cutoff::Real=1e-8, maxdim::Int=200,
position::Symbol=:pre)
+ _require_binary_position_space(H, "add_spin!")
H.spin_s === nothing || return H
spin_s = spin_index()
if position === :pre
@@ -960,6 +1092,7 @@ function add_zeeman!(H::TBHamiltonian, h;
tol::Real = 1e-8,
maxdim::Int = 200,
position::Union{Nothing,Symbol} = nothing)
+ _require_binary_position_space(H, "add_zeeman!")
direction in (:x, :y, :z) ||
error("direction must be :x, :y, or :z; got :$direction")
pos = something(position, H.aux_side)
@@ -1036,6 +1169,7 @@ function add_superconductivity!(H::TBHamiltonian, Δ;
tol::Real = 1e-8,
maxdim::Int = 200,
position::Union{Nothing,Symbol} = nothing)
+ _require_binary_position_space(H, "add_superconductivity!")
H.nambu_s === nothing ||
error("BdG already applied (H.nambu_s is set). Cannot apply twice.")
@@ -1151,6 +1285,7 @@ function add_soc!(H::TBHamiltonian, λ;
tol::Real = 1e-8,
maxdim::Int = 200,
position::Union{Nothing,Symbol} = nothing)
+ _require_binary_position_space(H, "add_soc!")
pos = something(position, H.aux_side)
add_spin!(H; cutoff=tol, maxdim=maxdim, position=pos)
pos_s = _pos_sites(H)
diff --git a/src/core/Utils.jl b/src/core/Utils.jl
index 57cbdf0..b1edbf3 100644
--- a/src/core/Utils.jl
+++ b/src/core/Utils.jl
@@ -521,6 +521,24 @@ function spatial_sampling_plan(L::Int;
ycenters = ny <= 1 ? [iy0] : round.(Int, range(iy0, iy1; length=ny))
centers = Int[ix + iy * Nx + 1 for iy in ycenters for ix in xcenters]
groups = [[c] for c in centers]
+ if num_avg > 1
+ # Spread `num_avg × num_avg` sub-samples ACROSS each coarse block
+ # (spacing stride÷num_avg), then average. On a coarse grid this is the
+ # correct block-average: it washes out fast on-site modulation (period ≪
+ # stride) while preserving the slow structure. Contrast `box_half`, a
+ # *contiguous* neighbourhood that spans only ±box_half cells — far less
+ # than one coarse stride — so it cannot average out sub-stride modulation
+ # and leaves per-pixel aliasing/speckle on a coarse grid.
+ sx = max(1, stride_x ÷ num_avg)
+ sy = max(1, stride_y ÷ num_avg)
+ groups = [
+ let uc0 = c - 1, ix_c = uc0 % Nx, iy_c = uc0 ÷ Nx
+ unique([mod(ix_c + a * sx, Nx) + mod(iy_c + b * sy, Ny) * Nx + 1
+ for b in 0:num_avg-1 for a in 0:num_avg-1])
+ end
+ for c in centers
+ ]
+ end
else
window = x_end - x_start + 1
nx = num_x <= 0 ? window : num_x
@@ -535,7 +553,9 @@ function spatial_sampling_plan(L::Int;
end
# ── 2D box averaging (periodic wrap) ───────────────────────────────────────
- if box_half > 0 && Lx !== nothing
+ # Skipped when num_avg>1 on a grid already spread each pixel into a sub-grid
+ # (the two averaging modes are mutually exclusive; num_avg takes precedence).
+ if box_half > 0 && num_avg <= 1 && Lx !== nothing
Nx = 2^Lx
Ny = 2^(L - Lx)
groups = [
@@ -625,6 +645,326 @@ function eval_mps_spatial(A::MPS;
return (values=values, centers=centers, groups=groups)
end
+"""
+ fibonacci_ldos_sampling_plan(L; depth=0, num_x=100, num_avg=1,
+ orientation=:standard, alignment=:atomic,
+ centered=true, origin=0)
+ fibonacci_ldos_sampling_plan(H::TBHamiltonian; kwargs...)
+
+Build a deterministic, bounded-size LDOS sampling plan for an `L`-qubit
+Fibonacci approximant, the projected-space counterpart of
+[`spatial_sampling_plan`](@ref). At `depth == 0` the plan covers the complete
+`F_(L+2)` conumber interval. Each additional depth selects the nested atomic
+renormalization window from [`fibonacci_rg_partition`](@ref), with
+`effective_L == L - 3depth`.
+
+The selected inherited interval is split into `min(num_x, window_count)`
+contiguous integer intervals whose widths differ by at most one. Up to
+`num_avg` equidistant conumbers (including both interval endpoints when there
+is more than one sample) are chosen in each interval and mapped directly to
+physical sites with [`fibonacci_site_from_conumber`](@ref). No full conumber
+permutation or other `F_(L+2)`-element array is constructed; storage is
+proportional to the requested output and sample counts.
+
+Returned fields useful to an LDOS/Slurm/HDF5 workflow include:
+
+- `groups`: physical-site vectors to pass as `x_groups` with
+ `ordering=:physical`;
+- `centers`: physical sites at the representative interval conumbers;
+- `conumber_axis`: those representative conumbers in the original `L`
+ coordinate system;
+- `intervals`, `interval_first`, and `interval_last`: represented inherited
+ conumber intervals;
+- `sample_conumbers`: the inherited conumbers corresponding to `groups`;
+- `sample_sites_flat`, `sample_conumbers_flat`, one-based `group_offsets`, and
+ zero-based `group_offsets_zero`: flat representations convenient for Julia
+ and Python/HDF5 consumers respectively (`flat[group_offsets[i]:
+ group_offsets[i+1]-1]` reconstructs group `i` in Julia);
+- `column_indices`: stable one-based output-column identifiers;
+- `depth`, `effective_L`, the original `L` and `N`, requested/actual sampling
+ counts, and all conumber conventions.
+
+`intervals` and every field containing `conumber` use the requested `centered`
+label convention. The corresponding `*_rank*` fields are always uncentered
+ranks in `0:N-1`. Thus a zoom always retains its original-`L` coordinates;
+the selected sites are never re-conumbered as an independent shorter chain.
+
+When `num_avg` is larger than an interval, that interval is sampled at every
+integer conumber and its `group_sizes` entry is smaller than `num_avg`.
+`num_x` in the result is the actual number of output columns, while
+`num_x_requested` records the input value.
+
+Nested (`depth > 0`) windows are defined only for the canonical atomic phase,
+so they require `alignment=:atomic` and `origin=0`. Reversing the orientation
+is supported because it maps the canonical atomic interval onto itself.
+
+The `TBHamiltonian` convenience method is defined in `position_spaces/Fibonacci.jl`.
+"""
+function fibonacci_ldos_sampling_plan(
+ L::Integer;
+ depth::Integer=0,
+ num_x::Integer=100,
+ num_avg::Integer=1,
+ orientation::Symbol=:standard,
+ alignment::Symbol=:atomic,
+ centered::Bool=true,
+ origin::Integer=0,
+)
+ num_x > 0 || throw(ArgumentError("num_x must be positive"))
+ num_avg > 0 || throw(ArgumentError("num_avg must be positive"))
+ orientation in (:standard, :reversed) ||
+ throw(ArgumentError("orientation must be :standard or :reversed"))
+ alignment in (:atomic, :raw) ||
+ throw(ArgumentError("alignment must be :atomic or :raw"))
+ if depth > 0 && (alignment !== :atomic || !iszero(origin))
+ throw(ArgumentError(
+ "depth > 0 requires alignment=:atomic and origin=0 so the selected " *
+ "window remains the canonical nested atomic renormalization window",
+ ))
+ end
+
+ # Work in uncentered ranks while partitioning. This keeps the RG embedding
+ # independent of how callers choose to label the inherited conumber axis.
+ partition = fibonacci_rg_partition(L; depth, centered=false)
+ N = fibonacci_site_count(L)
+ window_rank_first = first(partition.window_ranks)
+ window_rank_last = last(partition.window_ranks)
+ window_count = partition.window_count
+ ncolumns = min(Int(num_x), window_count)
+
+ # Tile the window exactly. Putting the remainder in the first intervals is
+ # deterministic and makes every width either floor(W/n) or ceil(W/n).
+ base_width, remainder = divrem(window_count, ncolumns)
+ rank_intervals = Vector{UnitRange{Int}}(undef, ncolumns)
+ cursor = window_rank_first
+ for column in 1:ncolumns
+ width = base_width + Int(column <= remainder)
+ rank_intervals[column] = cursor:(cursor + width - 1)
+ cursor += width
+ end
+ @assert cursor == window_rank_last + 1
+
+ shift = centered ? fld(N, 2) : 0
+ rank_to_conumber(rank::Int) = rank - shift
+ to_axis_interval(interval::UnitRange{Int}) =
+ rank_to_conumber(first(interval)):rank_to_conumber(last(interval))
+
+ intervals = [to_axis_interval(interval) for interval in rank_intervals]
+ interval_first = first.(intervals)
+ interval_last = last.(intervals)
+ interval_rank_first = first.(rank_intervals)
+ interval_rank_last = last.(rank_intervals)
+
+ # Integer samples are as uniformly spaced as possible. With one requested
+ # sample use the lower integer midpoint; with two or more include endpoints.
+ function equidistant_ranks(interval::UnitRange{Int})
+ width = length(interval)
+ count = min(Int(num_avg), width)
+ lo = first(interval)
+ count == 1 && return Int[lo + fld(width - 1, 2)]
+ return Int[lo + fld(k * (width - 1), count - 1)
+ for k in 0:(count - 1)]
+ end
+
+ sample_ranks = [equidistant_ranks(interval) for interval in rank_intervals]
+ sample_conumbers = [[rank_to_conumber(rank) for rank in ranks]
+ for ranks in sample_ranks]
+ center_ranks = Int[first(interval) + fld(length(interval) - 1, 2)
+ for interval in rank_intervals]
+ conumber_axis = rank_to_conumber.(center_ranks)
+
+ site_from_rank(rank::Int) = fibonacci_site_from_conumber(
+ L, rank_to_conumber(rank);
+ orientation, alignment, centered, origin,
+ )
+ groups = [[site_from_rank(rank) for rank in ranks] for ranks in sample_ranks]
+ centers = site_from_rank.(center_ranks)
+
+ group_sizes = length.(groups)
+ group_offsets = Vector{Int}(undef, ncolumns + 1)
+ group_offsets[1] = 1
+ for column in 1:ncolumns
+ group_offsets[column + 1] = group_offsets[column] + group_sizes[column]
+ end
+ group_offsets_zero = group_offsets .- 1
+ group_offsets_base = 1
+ total_samples = group_offsets[end] - 1
+ sample_sites_flat = Vector{Int}(undef, total_samples)
+ sample_conumbers_flat = Vector{Int}(undef, total_samples)
+ sample_ranks_flat = Vector{Int}(undef, total_samples)
+ for column in 1:ncolumns
+ destination = group_offsets[column]:(group_offsets[column + 1] - 1)
+ sample_sites_flat[destination] = groups[column]
+ sample_conumbers_flat[destination] = sample_conumbers[column]
+ sample_ranks_flat[destination] = sample_ranks[column]
+ end
+
+ window_first = rank_to_conumber(window_rank_first)
+ window_last = rank_to_conumber(window_rank_last)
+ metadata = (;
+ format="TensorBinding.fibonacci_ldos_sampling_plan",
+ format_version=1,
+ L=Int(L),
+ N,
+ depth=Int(depth),
+ effective_L=partition.effective_L,
+ window_count,
+ window_first,
+ window_last,
+ window_rank_first,
+ window_rank_last,
+ num_x=ncolumns,
+ num_x_requested=Int(num_x),
+ num_avg=Int(num_avg),
+ total_samples,
+ group_offsets_base,
+ orientation=String(orientation),
+ alignment=String(alignment),
+ centered,
+ origin=Int(origin),
+ )
+
+ return (;
+ groups,
+ centers,
+ conumber_axis,
+ intervals,
+ interval_first,
+ interval_last,
+ interval_rank_first,
+ interval_rank_last,
+ sample_conumbers,
+ sample_ranks,
+ sample_sites_flat,
+ sample_conumbers_flat,
+ sample_ranks_flat,
+ group_offsets,
+ group_offsets_zero,
+ group_offsets_base,
+ group_sizes,
+ column_indices=collect(1:ncolumns),
+ L=Int(L),
+ N,
+ depth=Int(depth),
+ effective_L=partition.effective_L,
+ window_count,
+ window_first,
+ window_last,
+ window_rank_first,
+ window_rank_last,
+ num_x=ncolumns,
+ num_x_requested=Int(num_x),
+ num_avg=Int(num_avg),
+ total_samples,
+ orientation,
+ alignment,
+ centered,
+ origin=Int(origin),
+ metadata,
+ )
+end
+
+# Enumerate all block members for exciton block-reduce (positional averaging).
+# For :block, spatial_sampling_plan gives singleton groups; this expands each to the
+# full set of probe positions inside the coarse block, enumerated from plan.stride_x/y.
+function _exciton_block_groups(plan, Lx::Union{Nothing,Int}, L::Int)
+ nblocks = length(plan.centers)
+ Wx = plan.stride_x
+ if Lx === nothing
+ return [[ixp * Wx + d + 1 for d in 0:Wx-1] for ixp in 0:nblocks-1]
+ end
+ a = plan.a
+ Wy = plan.stride_y
+ Nx = 2^Lx
+ return [let ixp = (iblock-1) % 2^a, iyp = (iblock-1) ÷ 2^a
+ [ixp*Wx + dx + (iyp*Wy + dy)*Nx + 1 for dy in 0:Wy-1 for dx in 0:Wx-1]
+ end
+ for iblock in 1:nblocks]
+end
+
+"""
+ ilinspace(xmin, xmax, num_x) -> Vector{Int}
+
+Return `num_x` as almost evenly spaced integers in `[xmin, xmax]`, inclusive,
+with a preference for the endpoints. Used to build the k-point center
+grid for band-structure sampling.
+"""
+function ilinspace(xmin, xmax, num_x::Int)
+ xvals = xmin:xmax
+ _N = length(xvals)
+ @assert 1 ≤ num_x ≤ _N
+ num_x == 1 && return [0]
+ step = (_N - 1) ÷ (num_x - 1)
+ return collect(xmin:step:(xmin+step*(num_x-1)))
+end
+
+"""
+ kspace_sampling_plan(L_pos, D; num_x, num_y=10, num_avg=1,
+ xmin=0, xmax=nothing, ymin=0, ymax=nothing,
+ k_groups_override=nothing) -> (; k_groups, num_x)
+
+Momentum-space sampling plan shared by `get_bands` and `get_bands_gpu`, the
+k-space counterpart of [`spatial_sampling_plan`](@ref). Momenta are 0-indexed
+QFT register labels `k in 0:2^L_pos-1`.
+
+- `k_groups_override` (for example from `kpath_2d`) is passed through untouched
+ and sets `num_x = length(k_groups)`.
+- `D == 1`: `num_x` centres from [`ilinspace`](@ref) over `[xmin, xmax]`
+ (`xmax` defaults to `2^L_pos - 1`); with `num_avg > 1` each centre is widened
+ to `num_avg` equidistant offsets within half a step on either side, clamped
+ to the register.
+- `D == 2`: `Lx = L_pos ÷ 2`; the first `min(num_x, 2^Lx)` points of the
+ `ilinspace` grids in `x` and `y` are zipped diagonally into row-major labels
+ `(y << Lx) | x`, again with optional `num_avg` widening. This is the legacy
+ diagonal cut through the 2D zone; for high-symmetry paths use `kpath_2d`.
+"""
+function kspace_sampling_plan(L_pos::Int, D::Int;
+ num_x::Int,
+ num_y::Int = 10,
+ num_avg::Int = 1,
+ xmin::Int = 0,
+ xmax = nothing,
+ ymin::Int = 0,
+ ymax = nothing,
+ k_groups_override = nothing)
+ if !isnothing(k_groups_override)
+ return (; k_groups = k_groups_override, num_x = length(k_groups_override))
+ end
+ N = 2^L_pos
+ if D == 1
+ _xmax = xmax === nothing ? N - 1 : Int(xmax)
+ xcenters = ilinspace(xmin, _xmax, num_x)
+ half_step = num_x > 1 ? (_xmax - xmin) / (2 * num_x) : 0
+ offsets = num_avg > 1 ? round.(Int, range(-half_step, half_step; length=num_avg)) : Int[0]
+ k_groups = [clamp.(xcenters[i] .+ offsets, 0, N - 1) for i in 1:num_x]
+ return (; k_groups, num_x)
+ elseif D == 2
+ Lx = div(L_pos, 2)
+ Nx_loc = 2^Lx
+ Ny_loc = 2^(L_pos - Lx)
+ nx = min(num_x, Nx_loc) # can't have more output pts than grid positions
+ _xmax = xmax === nothing ? Nx_loc - 1 : Int(xmax)
+ _ymax = ymax === nothing ? Ny_loc - 1 : Int(ymax)
+ xcenters = ilinspace(xmin, _xmax, Nx_loc)
+ ycenters = ilinspace(ymin, _ymax, Ny_loc)
+ half_step_x = nx > 1 ? (_xmax - xmin) / (2 * nx) : 0
+ half_step_y = num_y > 1 ? (_ymax - ymin) / (2 * num_y) : 0
+ x_offs = num_avg > 1 ? round.(Int, range(-half_step_x, half_step_x; length=num_avg)) : Int[0]
+ y_offs = num_avg > 1 ? round.(Int, range(-half_step_y, half_step_y; length=num_avg)) : Int[0]
+ k_groups = [
+ begin
+ xs = clamp.(xcenters[i] .+ x_offs, 0, Nx_loc - 1)
+ ys = clamp.(ycenters[i] .+ y_offs, 0, Ny_loc - 1)
+ [(y << Lx) | x for (x, y) in zip(xs, ys)] # diagonal zip in 2D k-space
+ end
+ for i in 1:nx
+ ]
+ return (; k_groups, num_x = nx)
+ else
+ error("kspace_sampling_plan: D must be 1 or 2")
+ end
+end
+
"""
rms_error(a, b) -> Float64
@@ -705,6 +1045,31 @@ function get_diagonal_mpo(L, sites, f; type=Float64, tol::Real=1e-8)
end
+"""
+ extract_diagonal_to_mps(M) -> MPS
+
+Extract the diagonal of an MPO `M` as an MPS by projecting each local bra/ket
+pair onto equal physical values. This is shared by KPM trace/LDOS, SCF, RPA,
+QFT, and purification routines.
+"""
+function extract_diagonal_to_mps(M::MPO)::MPS
+ N = length(M)
+ new_tensors = Vector{ITensor}(undef, N)
+ for i in 1:N
+ tensor = M[i]
+ bra, ket = siteinds(M, i)
+ diagonal_inds = uniqueinds(tensor, ket, bra)
+ result = ITensor(diagonal_inds..., ket)
+ for value in 1:dim(ket)
+ slice = tensor * onehot(ket => value) * onehot(bra => value)
+ result += slice * onehot(ket => value)
+ end
+ new_tensors[i] = result
+ end
+ return MPS(new_tensors)
+end
+
+
# ---------------------------------------------------------------------
# MPS -> diagonal MPO conversion
diff --git a/src/gpu/GPU_tk.jl b/src/gpu/GPU_tk.jl
index 561c6a8..78088ad 100644
--- a/src/gpu/GPU_tk.jl
+++ b/src/gpu/GPU_tk.jl
@@ -20,6 +20,9 @@
# get_ldos_spatial_gpu(H, Ncheb, ω; reduce=..., ...) — A(r,ω) real-space LDOS
# (:point or :block sampling,
# sublattice :average/:resolve)
+# get_ldos_spatial_mps_gpu(H, Ncheb, ω; ...) — A(r,ω), independent
+# GPU MPS recursions (including
+# projected position spaces)
# get_dos_stochastic_gpu(H, Ncheb, ω; ...) — stochastic-trace DOS
# get_nh_dos_grid_gpu(H, xlims, nx, ylims, ny, n; ...) — NH stochastic DOS
# get_nh_dos_points_gpu(H, z_points, n; ...) — NH stochastic DOS at selected z
@@ -1142,6 +1145,7 @@ function get_bands_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
type::Type{<:Number} = ComplexF32,
dtype::Union{Nothing,Type{<:Number}} = nothing)
+ _require_binary_position_space(H, "get_bands_gpu")
_check_gpu("get_bands_gpu")
gpu_type = dtype === nothing ? type : dtype
dtype !== nothing && dtype != type && type != ComplexF32 &&
@@ -1184,40 +1188,12 @@ function get_bands_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
kpath_setup(kpath_lattice, Lx_kp, Ly_kp, kpath; npts_per_segment=num_x)
end
- # ── k-groups (same logic as low-level CPU get_bands) ────────────────────
- Lx_pos = D == 2 ? div(L_pos, 2) : 0
- N_pos = 2^L_pos
- if !isnothing(k_groups_override)
- k_groups = k_groups_override
- num_x = length(k_groups)
- elseif D == 1
- _xmax = xmax === nothing ? N_pos - 1 : Int(xmax)
- xcenters = ilinspace(xmin, _xmax, num_x)
- half_step = num_x > 1 ? (_xmax - xmin) / (2 * num_x) : 0
- offsets = num_avg > 1 ? round.(Int, range(-half_step, half_step; length=num_avg)) : Int[0]
- k_groups = [clamp.(xcenters[i] .+ offsets, 0, N_pos - 1) for i in 1:num_x]
- elseif D == 2
- Nx_loc = 2^Lx_pos; Ny_loc = 2^(L_pos - Lx_pos)
- num_x = min(num_x, Nx_loc)
- _xmax = xmax === nothing ? Nx_loc - 1 : Int(xmax)
- _ymax = ymax === nothing ? Ny_loc - 1 : Int(ymax)
- xcenters = ilinspace(xmin, _xmax, Nx_loc)
- ycenters = ilinspace(ymin, _ymax, Ny_loc)
- hsx = num_x > 1 ? (_xmax - xmin) / (2 * num_x) : 0
- hsy = num_y > 1 ? (_ymax - ymin) / (2 * num_y) : 0
- x_offs = num_avg > 1 ? round.(Int, range(-hsx, hsx; length=num_avg)) : Int[0]
- y_offs = num_avg > 1 ? round.(Int, range(-hsy, hsy; length=num_avg)) : Int[0]
- k_groups = [
- begin
- xs = clamp.(xcenters[i] .+ x_offs, 0, Nx_loc - 1)
- ys = clamp.(ycenters[i] .+ y_offs, 0, Ny_loc - 1)
- [(y << Lx_pos) | x for (x, y) in zip(xs, ys)]
- end
- for i in 1:num_x
- ]
- else
- error("D must be 1 or 2")
- end
+ # ── k-groups (shared planner in core/Utils.jl, same as CPU get_bands) ────
+ Lx_pos = D == 2 ? div(L_pos, 2) : 0
+ kplan = kspace_sampling_plan(L_pos, D; num_x, num_y, num_avg,
+ xmin, xmax, ymin, ymax, k_groups_override)
+ k_groups = kplan.k_groups
+ num_x = kplan.num_x
Ak_w = zeros(Float64, Nω, num_x)
@@ -1447,6 +1423,8 @@ function get_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
type::Type{<:Number} = ComplexF32,
dtype::Union{Nothing,Type{<:Number}} = nothing)
+ _require_binary_position_space(H, "get_ldos_spatial_gpu")
+
_check_gpu("get_ldos_spatial_gpu")
gpu_type = dtype === nothing ? type : dtype
dtype !== nothing && dtype != type && type != ComplexF32 &&
@@ -1625,6 +1603,304 @@ function get_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
end
+"""
+ _reconstruct_ldos_moment_columns(moments, W, denom, valid)
+ -> Matrix{Float64}
+
+Reconstruct one LDOS column per column of raw Chebyshev `moments`. The weight
+matrix follows `_dos_weight_matrix`: `W[n, iω]` multiplies moment order `n-1`,
+and `denom[iω]` supplies the kernel-specific normalization. Invalid energies
+are returned as zero columns in energy space.
+"""
+function _reconstruct_ldos_moment_columns(
+ moments::AbstractMatrix{<:Real},
+ W::AbstractMatrix{<:Real},
+ denom::AbstractVector{<:Real},
+ valid::AbstractVector{Bool},
+)
+ Ncheb, ncols = size(moments)
+ size(W, 1) == Ncheb || throw(DimensionMismatch(
+ "moment rows ($(size(moments, 1))) must match weight rows ($(size(W, 1))).",
+ ))
+ Nω = size(W, 2)
+ length(denom) == Nω || throw(DimensionMismatch(
+ "denominator length ($(length(denom))) must match energy count ($Nω).",
+ ))
+ length(valid) == Nω || throw(DimensionMismatch(
+ "valid-mask length ($(length(valid))) must match energy count ($Nω).",
+ ))
+
+ result = zeros(Float64, Nω, ncols)
+ mul!(result, transpose(W), moments)
+ for iω in 1:Nω
+ if valid[iω]
+ view(result, iω, :) ./= denom[iω]
+ else
+ fill!(view(result, iω, :), 0.0)
+ end
+ end
+ return result
+end
+
+
+"""
+ get_ldos_spatial_mps_gpu(H, Ncheb, ω_phys_vals;
+ x_groups=nothing,
+ num_x=min(H.N, 100), num_avg=1,
+ x_start=1, x_end=H.N,
+ kernel=:jackson, lambda=4.0, eta=0.0, m_order=4,
+ maxdim=100, cutoff=1e-8,
+ type=ComplexF32, dtype=nothing,
+ verbose=false, printinfo=false,
+ return_maxlinkdim=false,
+ return_moments=false)
+ -> Matrix{Float64}
+
+GPU spatial LDOS from one independent MPS Chebyshev recursion per physical-site
+probe. Unlike [`get_ldos_spatial_gpu`](@ref), this path does not construct an MPO
+Chebyshev series and supports projected position spaces such as
+`FibonacciPositionSpace`.
+
+The rescaled operator is `H̃ = (H - H.center * P) / H.scale`, where
+`P = physical_projector(H)`. Probe `x` is constructed with
+`physical_site_state(H, x)`, so `x` is always a 1-based *physical* site rather
+than an ambient tensor-register index.
+
+`x_groups` can be a vector of positions (one output column per position) or a
+vector of position vectors. In the latter case, all probe LDOS values in a group
+are averaged into one output column. Without explicit groups, `num_x` intervals
+over `x_start:x_end` are sampled with `num_avg` approximately equidistant probes
+per interval. Automatic planning allocates only `O(num_x * num_avg)` probe
+indices, so callers can sample a huge projected space without enumerating it by
+choosing a modest `num_x` (or by supplying `x_groups`).
+The default is at most 100 output columns.
+
+`kernel=:hodc` uses HODC reconstruction (`eta`, `m_order`; `eta=0` uses
+`1/(Ncheb+1)`). Other supported kernels are `:jackson`, `:lorentz` (`lambda`),
+`:fejer`, and `:dirichlet`.
+
+Use `type=ComplexF32` (default) or a supported real/complex GPU tensor type;
+`dtype` is an alias. With `return_moments=true`, the group-averaged raw
+Chebyshev moments are also returned as a `Matrix{Float64}` of size
+`(Ncheb, length(x_groups))` (or `(Ncheb, num_x)` for automatic groups):
+
+`moments[n, j] = mean(x -> real(), group[j])`,
+
+where `Htilde = (H - H.center * P) / H.scale`. These moments contain no kernel
+weights or energy-dependent normalization, and can therefore be reconstructed
+later on a different energy grid or with a different KPM kernel.
+
+With `return_maxlinkdim=true`, `linkdims[j]` is the largest MPS bond dimension
+reached by any probe in group `j`. Return values are unambiguous for all keyword
+combinations:
+
+- neither keyword: `ldos`
+- `return_maxlinkdim=true`: `(ldos, linkdims)` (the existing API)
+- `return_moments=true`: `(ldos, moments)`
+- both keywords: `(ldos, moments, linkdims)`
+
+This entry point intentionally supports position-only, one-dimensional point or
+explicit-group sampling. Grid/window/box/block sampling, non-physical ordering,
+and auxiliary degrees of freedom are rejected with targeted errors. For those
+features use the MPO GPU path or the CPU `get_ldos_spatial` implementation.
+"""
+function get_ldos_spatial_mps_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
+ x_groups = nothing,
+ num_x::Int = min(H.N, 100),
+ num_avg::Int = 1,
+ x_start::Int = 1,
+ x_end::Int = H.N,
+ kernel::Symbol = :jackson,
+ lambda::Real = 4.0,
+ eta::Real = 0.0,
+ m_order::Int = 4,
+ maxdim::Int = 100,
+ cutoff::Real = 1e-8,
+ type::Type{<:Number} = ComplexF32,
+ dtype::Union{Nothing,Type{<:Number}} = nothing,
+ verbose::Bool = false,
+ printinfo::Bool = false,
+ return_maxlinkdim::Bool = false,
+ return_moments::Bool = false,
+ # Accepted only to provide clear compatibility errors.
+ num_y = nothing,
+ grid::Bool = false,
+ xwin = nothing,
+ ywin = nothing,
+ box_half::Int = 0,
+ reduce::Symbol = :point,
+ ordering::Symbol = :physical,
+ sublattice::Symbol = :auto,
+ nambu_proj::Bool = false,
+ proj_nambu = nothing,
+ spin_proj::Bool = false,
+ proj_s = nothing,
+ layer_proj::Bool = false,
+ proj_layer = nothing,
+ sublat_proj::Bool = false,
+ proj_sl = nothing)
+
+ Ncheb >= 2 || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: Ncheb must be at least 2."
+ ))
+ reduce === :point || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: only reduce=:point is supported; " *
+ "block reduction belongs to the MPO GPU path."
+ ))
+ if grid || num_y !== nothing || xwin !== nothing || ywin !== nothing || box_half != 0
+ throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: grid, num_y, windows, and box averaging " *
+ "are unsupported. Supply 1-based physical positions through x_groups."
+ ))
+ end
+ ordering === :physical || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: only ordering=:physical is supported. " *
+ "Map alternate coordinates to physical sites before passing x_groups."
+ ))
+ sublattice === :auto || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: sublattice resolution/averaging is unsupported."
+ ))
+
+ aux_requested = nambu_proj || spin_proj || layer_proj || sublat_proj ||
+ proj_nambu !== nothing || proj_s !== nothing ||
+ proj_layer !== nothing || proj_sl !== nothing
+ has_aux = !isnothing(H.nambu_s) || !isnothing(H.spin_s) ||
+ !isnothing(H.layer_s) || !isnothing(H.sublattice_s) ||
+ length(H.sites) != H.L
+ (aux_requested || has_aux) && throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: only position-only Hamiltonians are supported; " *
+ "auxiliary degrees of freedom and auxiliary projections are not available " *
+ "on this MPS GPU path."
+ ))
+
+ groups = if x_groups !== nothing
+ x_groups isa AbstractVector{<:AbstractVector} ?
+ [collect(Int, group) for group in x_groups] :
+ [[Int(x)] for x in x_groups]
+ else
+ num_x > 0 || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: num_x must be positive."
+ ))
+ num_avg > 0 || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: num_avg must be positive."
+ ))
+ 1 <= x_start <= x_end <= H.N || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: expected 1 <= x_start <= x_end <= H.N."
+ ))
+ window = x_end - x_start + 1
+ num_x <= window || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: num_x=$num_x exceeds the sampling " *
+ "window length $window."
+ ))
+ [let
+ lo = x_start + fld((i - 1) * window, num_x)
+ hi = x_start + fld(i * window, num_x) - 1
+ nsample = min(num_avg, hi - lo + 1)
+ nsample == 1 ? Int[lo] :
+ unique(round.(Int, range(lo, hi; length=nsample)))
+ end for i in 1:num_x]
+ end
+
+ isempty(groups) && throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: no spatial groups were selected."
+ ))
+ for group in groups
+ isempty(group) && throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: spatial groups must not be empty."
+ ))
+ all(x -> 1 <= x <= H.N, group) || throw(ArgumentError(
+ "get_ldos_spatial_mps_gpu: every position must lie in 1:H.N."
+ ))
+ end
+
+ _check_gpu("get_ldos_spatial_mps_gpu")
+ gpu_type = _resolve_gpu_type(
+ "get_ldos_spatial_mps_gpu", type, dtype, cutoff,
+ )
+ _ensure_scale!(H)
+
+ # P, rather than the ambient identity, is essential for projected position
+ # spaces: invalid register states must remain zero under the spectral shift.
+ P_cpu = physical_projector(H)
+ Ham_n_cpu = (1 / H.scale) * +(
+ H.mpo, (-H.center) * P_cpu; cutoff=Float64(cutoff),
+ )
+ Ham_n_gpu = _to_gpu_mpo(Ham_n_cpu, gpu_type)
+
+ ω_vals = (collect(ω_phys_vals) .- H.center) ./ H.scale
+ Nω = length(ω_vals)
+ W, denom = _dos_weight_matrix(
+ Ncheb, ω_vals; kernel=kernel, lambda=lambda, eta=eta, m_order=m_order,
+ )
+ valid = [abs(ω) < 1.0 for ω in ω_vals]
+ # Store kernel-independent, group-averaged moments. Besides making them
+ # available for offline reconstruction, this avoids applying all Nω
+ # energy weights separately for every probe in an averaged group.
+ moments = zeros(Float64, Ncheb, length(groups))
+ linkdims = zeros(Int, length(groups))
+
+ apply_kwargs = (cutoff=Float64(cutoff), maxdim=maxdim)
+ two = gpu_type(2)
+ negone = gpu_type(-1)
+ printinfo && println(
+ " [gpu] spatial MPS LDOS dtype=$gpu_type, groups=$(length(groups)), " *
+ "projected=$( !(H.position_space isa BinaryPositionSpace) )",
+ )
+
+ for (j, group) in enumerate(groups)
+ group_moments = view(moments, :, j)
+ group_weight = inv(Float64(length(group)))
+ group_maxlinkdim = 0
+
+ for x in group
+ psi0_gpu = _to_gpu_mps(physical_site_state(H, x), gpu_type)
+
+ function kpm_step!(phi, n)
+ mu = Float64(real(inner(psi0_gpu, phi)))
+ group_moments[n] += group_weight * mu
+ end
+
+ phi_km2 = psi0_gpu
+ kpm_step!(phi_km2, 1)
+ group_maxlinkdim = max(group_maxlinkdim, maxlinkdim(phi_km2))
+
+ phi_km1 = apply(Ham_n_gpu, phi_km2; apply_kwargs...)
+ kpm_step!(phi_km1, 2)
+ group_maxlinkdim = max(group_maxlinkdim, maxlinkdim(phi_km1))
+
+ for k in 3:Ncheb
+ phi_k = +(
+ two * apply(Ham_n_gpu, phi_km1; apply_kwargs...),
+ negone * phi_km2;
+ apply_kwargs...,
+ )
+ kpm_step!(phi_k, k)
+ group_maxlinkdim = max(group_maxlinkdim, maxlinkdim(phi_k))
+ phi_km2 = phi_km1
+ phi_km1 = phi_k
+ end
+
+ _gpu_gc!()
+ end
+
+ linkdims[j] = group_maxlinkdim
+ (verbose || printinfo) && (j % 5 == 0 || j == length(groups)) &&
+ println(
+ " [gpu] spatial MPS LDOS $j/$(length(groups)) " *
+ "(x=$(first(group)), n_avg=$(length(group))) " *
+ "maxlinkdim=$group_maxlinkdim",
+ )
+ end
+
+ result = _reconstruct_ldos_moment_columns(moments, W, denom, valid)
+
+ if return_moments
+ return return_maxlinkdim ? (result, moments, linkdims) : (result, moments)
+ end
+ return return_maxlinkdim ? (result, linkdims) : result
+end
+
+
"""
get_dos_stochastic_gpu(H, Ncheb, ω_phys_vals; kwargs...)
-> Vector{Float64} length Nω
@@ -1675,6 +1951,7 @@ function get_dos_stochastic_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
type::Type{<:Number} = ComplexF32,
dtype::Union{Nothing,Type{<:Number}} = nothing)
+ _require_binary_position_space(H, "get_dos_stochastic_gpu")
_check_gpu("get_dos_stochastic_gpu")
gpu_type = _resolve_gpu_type("get_dos_stochastic_gpu", type, dtype, cutoff)
_ensure_scale!(H)
@@ -2450,23 +2727,8 @@ function get_nh_dos_grid_diag_trace_gpu(H::TBHamiltonian, xlims, nx::Int, ylims,
end
-# Enumerate all block members for exciton block-reduce (positional averaging).
-# For :block, spatial_sampling_plan gives singleton groups; this expands each to the
-# full set of probe positions inside the coarse block, enumerated from plan.stride_x/y.
-function _exciton_block_groups(plan, Lx::Union{Nothing,Int}, L::Int)
- nblocks = length(plan.centers)
- Wx = plan.stride_x
- if Lx === nothing
- return [[ixp * Wx + d + 1 for d in 0:Wx-1] for ixp in 0:nblocks-1]
- end
- a = plan.a
- Wy = plan.stride_y
- Nx = 2^Lx
- return [let ixp = (iblock-1) % 2^a, iyp = (iblock-1) ÷ 2^a
- [ixp*Wx + dx + (iyp*Wy + dy)*Nx + 1 for dy in 0:Wy-1 for dx in 0:Wx-1]
- end
- for iblock in 1:nblocks]
-end
+# Block-member enumeration for exciton block-reduce (`_exciton_block_groups`)
+# lives in core/Utils.jl next to spatial_sampling_plan.
"""
get_exciton_ldos_spatial_gpu(H, Ncheb, ω_phys_vals;
@@ -2497,6 +2759,11 @@ Use `type=ComplexF32` (default, faster) or `type=ComplexF64` (safer at tight cut
or on large systems where F32 eigendecomposition can produce NaN). `dtype` is accepted
as an alias for `type` for consistency with other GPU entry points.
+`return_maxlinkdim=true` returns `(result, linkdims)` instead of just `result`, where
+`linkdims::Vector{Int}` is the reached MPS bond dimension per output column (the χ the
+Chebyshev recursion hit under the given `maxdim`/`cutoff`). Useful for cutoff/tolerance
+studies where χ is the observable.
+
!!! note "Block averaging not supported"
`reduce=:block` is **not available** for the exciton LDOS. In the MPO-based LDOS
functions (`get_ldos_spatial_gpu`), block averaging is a cheap O(1) partial trace
@@ -2526,7 +2793,8 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
type::Type{<:Number} = ComplexF32,
dtype::Union{Nothing,Type{<:Number}} = nothing,
verbose::Bool = false,
- printinfo::Bool = false)
+ printinfo::Bool = false,
+ return_maxlinkdim::Bool = false)
_check_gpu("get_exciton_ldos_spatial_gpu")
gpu_type = dtype === nothing ? type : dtype
@@ -2595,6 +2863,8 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
printinfo && println(" [gpu] exciton ldos dtype=$gpu_type")
+ linkdims = zeros(Int, nX) # reached MPS bond dim per output column (see return_maxlinkdim)
+
for (j, group) in enumerate(groups)
last_linkdim = 0
@@ -2634,11 +2904,12 @@ function get_exciton_ldos_spatial_gpu(H::TBHamiltonian, Ncheb::Int, ω_phys_vals
for iω in 1:Nω
result[iω, j] /= length(group)
end
+ linkdims[j] = last_linkdim
(verbose || printinfo) && (j % 5 == 0 || j == nX) &&
println(" [gpu] exciton ldos $j/$nX (X=$(Xs[j]), n_avg=$(length(group))) maxlinkdim=$last_linkdim")
end
- return result
+ return return_maxlinkdim ? (result, linkdims) : result
end
@@ -2805,6 +3076,7 @@ function get_C_gpu(H::TBHamiltonian, xfunc=nothing, yfunc=nothing;
dtype::Type{<:Complex} = ComplexF32,
printinfo::Bool = false)
+ _require_binary_position_space(H, "get_C_gpu")
_check_gpu("get_C_gpu")
gpu_type = _resolve_gpu_type("get_C_gpu", dtype, nothing, cutoff)
Λ_val = Lambda !== nothing ? Float64(Lambda) : Float64(Λ)
diff --git a/src/physics/NH_tk.jl b/src/physics/NH_tk.jl
index 0147f76..44a241b 100644
--- a/src/physics/NH_tk.jl
+++ b/src/physics/NH_tk.jl
@@ -74,6 +74,7 @@ function hermitized_hamiltonian(H::TBHamiltonian;
scale::Real = 0.0,
convention::Symbol = :z_minus_H,
block_placement::Symbol = :post)
+ _require_binary_position_space(H, "hermitized_hamiltonian")
convention in (:H_minus_z, :z_minus_H) ||
error("convention must be :H_minus_z or :z_minus_H; got :$convention")
I_H = MPO(H.sites, "Id")
@@ -323,6 +324,7 @@ function add_nh_onsite!(H::TBHamiltonian, v;
tol::Real = 1e-8,
maxdim::Int = 200,
type = ComplexF64)
+ _require_binary_position_space(H, "add_nh_onsite!")
pos_s = _nh_position_sites_only(H)
term = _nh_diagonal_mpo(H.L, pos_s, v; Lx=Lx, type=type)
H.mpo = +(H.mpo, term; cutoff=tol, maxdim=maxdim)
@@ -385,6 +387,7 @@ function add_loss!(H::TBHamiltonian, f;
maxdim::Int = 200,
type = Float64,
space::Symbol = :full)
+ _require_binary_position_space(H, "add_loss!")
term = ComplexF64(coefficient) * loss_profile_mpo(H, f; Lx=Lx, type=type, space=space)
H.mpo = +(H.mpo, term; cutoff=tol, maxdim=maxdim)
ITensorMPS.truncate!(H.mpo; cutoff=tol, maxdim=maxdim)
@@ -455,6 +458,7 @@ function nh_nonreciprocal_hopping_mpo(H::TBHamiltonian, t_forward, t_backward;
tol::Real = 1e-8,
maxdim::Int = 200,
type = ComplexF64)
+ _require_binary_position_space(H, "nh_nonreciprocal_hopping_mpo")
pos_s = _nh_position_sites_only(H)
Hf = _nh_directional_hop(pos_s, H.N, t_forward, nn, :forward;
L=H.L, tol=tol, maxdim=maxdim, type=type)
diff --git a/src/physics/Purification_tk.jl b/src/physics/Purification_tk.jl
index a64239e..3395e97 100644
--- a/src/physics/Purification_tk.jl
+++ b/src/physics/Purification_tk.jl
@@ -182,6 +182,7 @@ end
function purification_initial_guess(H::TBHamiltonian; ϵF::Real=0.0,
maxdim::Int=40, cutoff::Float64=1e-8)
+ _require_binary_position_space(H, "purification_initial_guess")
_ensure_scale!(H)
Id = MPO(H.sites, "Id")
coeff_I = 0.5 + (ϵF + H.center) / (2 * H.scale)
@@ -208,6 +209,7 @@ function mcweeny_purify(H::TBHamiltonian;
cutoff::Float64 = 1e-8,
tol::Float64 = 1e-5,
verbose::Bool = false)
+ _require_binary_position_space(H, "mcweeny_purify")
ρ0 = purification_initial_guess(H; ϵF=ϵF, maxdim=maxdim, cutoff=cutoff)
ρ = mcweeny_purify(ρ0; maxiters=maxiters, maxdim=maxdim, cutoff=cutoff,
tol=tol, verbose=verbose)
@@ -230,6 +232,7 @@ function sp2_purify(H::TBHamiltonian;
cutoff::Float64 = 1e-8,
tol::Float64 = 1e-5,
verbose::Bool = false)
+ _require_binary_position_space(H, "sp2_purify")
ρ0 = purification_initial_guess(H; maxdim=maxdim, cutoff=cutoff)
ρ = sp2_purify(ρ0, Nel; maxiters=maxiters, maxdim=maxdim, cutoff=cutoff,
tol=tol, verbose=verbose)
@@ -284,6 +287,8 @@ function get_density(H::TBHamiltonian;
tol::Float64 = 1e-5,
verbose::Bool = false)
+ method === :kpm || _require_binary_position_space(H, "get_density(method=:$method)")
+
if H._density_cache !== nothing
verbose && println("get_density: returning cached density matrix")
return H._density_cache
@@ -408,6 +413,7 @@ function get_ldos_drho(H::TBHamiltonian, ω::Real;
maxiters::Int = 30,
tol::Float64 = 1e-5,
verbose::Bool = false)
+ _require_binary_position_space(H, "get_ldos_drho")
mode in (:mpo, :mps) ||
error("get_ldos_drho: mode must be :mpo or :mps, got :$mode")
_ensure_scale!(H)
@@ -455,6 +461,7 @@ function get_dos_drho(H::TBHamiltonian, ω::Real;
maxiters::Int = 30,
tol::Float64 = 1e-5,
verbose::Bool = false)
+ _require_binary_position_space(H, "get_dos_drho")
_ensure_scale!(H)
ρ0_p = purification_initial_guess(H; ϵF = ω + dmu, maxdim=maxdim, cutoff=cutoff)
diff --git a/src/physics/QFT_tk.jl b/src/physics/QFT_tk.jl
index acc21cb..306372b 100644
--- a/src/physics/QFT_tk.jl
+++ b/src/physics/QFT_tk.jl
@@ -80,8 +80,8 @@
# 1a. Single-particle QFT conjugate_by_qft
# 1b. Exciton QFT conjugate_by_qft_exciton
# 2. Legacy sublattice projectors projop_2DSL, projop_1DSL
-# 3. Internal utilities ilinspace, _eval_diag_mps, sample_diag,
-# _kpm_weight_matrix
+# 3. Internal utilities _eval_diag_mps, sample_diag, _kpm_weight_matrix
+# (ilinspace / kspace_sampling_plan: core/Utils.jl)
# (exciton MPS probes mpsexciton/Q/QTrace/KQ now live in TwoParticle_tk.jl)
# 3b. High-symmetry k-path kpath_2d, hsk_honeycomb/square/triangular,
# kpath_setup, _hs_label, _hsk
@@ -351,21 +351,8 @@ end
# _kpm_weight_matrix — precomputed Chebyshev-KPM weights W[n, iω]
# ============================================================
-"""
- ilinspace(xmin, xmax, num_x) -> Vector{Int}
-
-Return `num_x` as almost evenly spaced integers in `[xmin, xmax]`, inclusive,
-with a preference for the endpoints. Used to build the k-point center
-grid for band-structure sampling.
-"""
-function ilinspace(xmin, xmax, num_x::Int)
- xvals = xmin:xmax
- _N = length(xvals)
- @assert 1 ≤ num_x ≤ _N
- num_x == 1 && return [0]
- step = (_N - 1) ÷ (num_x - 1)
- return collect(xmin:step:(xmin+step*(num_x-1)))
-end
+# `ilinspace` and `kspace_sampling_plan` (k-point centre placement and grouping
+# shared with get_bands_gpu) live in core/Utils.jl with the other sampling plans.
"""
@@ -822,41 +809,12 @@ function get_bands(H_mpo::MPO, scale::Real, center::Real, sites,
# compute it unconditionally so it is always in scope when D==2.
Lx = D == 2 ? div(L_pos, 2) : 0
- # ── Build k-point groups ──────────────────────────────────────────────────
+ # ── Build k-point groups (shared planner in core/Utils.jl) ────────────────
# k_groups_override (from kpath_2d) bypasses the grid sampling entirely.
- if !isnothing(k_groups_override)
- k_groups = k_groups_override
- num_x = length(k_groups)
- elseif D == 1
- _xmax = xmax === nothing ? N - 1 : Int(xmax)
- xcenters = ilinspace(xmin, _xmax, num_x)
- half_step = num_x > 1 ? (_xmax - xmin) / (2 * num_x) : 0
- offsets = num_avg > 1 ? round.(Int, range(-half_step, half_step; length=num_avg)) : Int[0]
- k_groups = [clamp.(xcenters[i] .+ offsets, 0, N - 1) for i in 1:num_x]
- elseif D == 2
- Lx = div(L_pos, 2) # also computed above; repeated here keeps the branch self-contained
- Nx_loc = 2^Lx
- Ny_loc = 2^(L_pos - Lx)
- num_x = min(num_x, Nx_loc) # can't have more output pts than grid positions
- _xmax = xmax === nothing ? Nx_loc - 1 : Int(xmax)
- _ymax = ymax === nothing ? Ny_loc - 1 : Int(ymax)
- xcenters = ilinspace(xmin, _xmax, Nx_loc)
- ycenters = ilinspace(ymin, _ymax, Ny_loc)
- half_step_x = num_x > 1 ? (_xmax - xmin) / (2 * num_x) : 0
- half_step_y = num_y > 1 ? (_ymax - ymin) / (2 * num_y) : 0
- x_offs = num_avg > 1 ? round.(Int, range(-half_step_x, half_step_x; length=num_avg)) : Int[0]
- y_offs = num_avg > 1 ? round.(Int, range(-half_step_y, half_step_y; length=num_avg)) : Int[0]
- k_groups = [
- begin
- xs = clamp.(xcenters[i] .+ x_offs, 0, Nx_loc - 1)
- ys = clamp.(ycenters[i] .+ y_offs, 0, Ny_loc - 1)
- [(y << Lx) | x for (x, y) in zip(xs, ys)] # diagonal zip in 2D k-space
- end
- for i in 1:num_x
- ]
- else
- error("D must be 1 or 2")
- end
+ kplan = kspace_sampling_plan(L_pos, D; num_x, num_y, num_avg,
+ xmin, xmax, ymin, ymax, k_groups_override)
+ k_groups = kplan.k_groups
+ num_x = kplan.num_x
Ak_w = zeros(Float64, Nω, num_x)
@@ -1036,6 +994,7 @@ function get_exciton_bands(H_QFT::MPO, H::TBHamiltonian, Ncheb::Int, omega_phys_
cutoff::Real = 1e-8,
verbose::Bool = false,
printinfo::Bool = false)
+ _require_binary_position_space(H, "get_exciton_bands")
_ensure_scale!(H)
length(H.sites) == 2 * H.L ||
error("get_exciton_bands: H is not an exciton Hamiltonian (expected length(H.sites) == 2*H.L).")
@@ -1060,9 +1019,7 @@ function get_exciton_bands(H_QFT::MPO, H::TBHamiltonian, Ncheb::Int, omega_phys_
q_end_eff = q_end === nothing ? k_end : Int(q_end)
groups = if group_arg !== nothing
- group_arg isa AbstractVector{<:AbstractVector} ?
- [collect(Int, grp) for grp in group_arg] :
- [[Int(q)] for q in group_arg]
+ spatial_sampling_plan(H.L; x_groups=group_arg).groups
elseif list_arg !== nothing
[[Int(q)] for q in list_arg]
else
@@ -1073,12 +1030,10 @@ function get_exciton_bands(H_QFT::MPO, H::TBHamiltonian, Ncheb::Int, omega_phys_
window = q_end_eff - q_start_eff + 1
num_q_eff <= window ||
error("get_exciton_bands: num_q=$num_q_eff exceeds sampling window length $window.")
- dq = div(window, num_q_eff)
- dq_sub = max(1, div(dq, num_avg))
- [[q_start_eff + (i - 1) * dq + k * dq_sub
- for k in 0:num_avg-1
- if q_start_eff + (i - 1) * dq + k * dq_sub <= q_end_eff]
- for i in 1:num_q_eff]
+ # 1D point layout of the shared planner (core/Utils.jl): stride
+ # window ÷ num_q with num_avg sub-probes per coarse cell.
+ spatial_sampling_plan(H.L; num_x=num_q_eff, num_avg,
+ x_start=q_start_eff, x_end=q_end_eff).groups
end
isempty(groups) && error("get_exciton_bands: no momentum groups were selected.")
@@ -1173,6 +1128,7 @@ function get_exciton_continuum(H_QFT::MPO, H::TBHamiltonian, Ncheb::Int, omega_p
cutoff::Real = 1e-8,
verbose::Bool = false,
printinfo::Bool = false)
+ _require_binary_position_space(H, "get_exciton_continuum")
_ensure_scale!(H)
length(H.sites) == 2 * H.L ||
error("get_exciton_continuum: H is not an exciton Hamiltonian (expected length(H.sites) == 2*H.L).")
@@ -1336,6 +1292,7 @@ function get_bands(H::TBHamiltonian, Ncheb::Int, D::Int, ω_phys_vals;
cutoff::Real = 1e-10,
printinfo::Bool = false)
+ _require_binary_position_space(H, "get_bands")
_ensure_scale!(H)
nambu_proj, spin_proj, layer_proj, sublat_proj =
_autoenable_proj(H, nambu_proj, spin_proj, layer_proj, sublat_proj)
@@ -1624,6 +1581,7 @@ end
function get_bands(H::TBHamiltonian, ω_phys_vals;
aux_proj = nothing, tol=1e-9, maxdim::Int=100)
+ _require_binary_position_space(H, "get_bands")
H._tn_cache === nothing &&
error("No Chebyshev cache found. Call KPM_Tn(H, Ncheb; ...) first.")
pos_sites = _pos_sites(H)
diff --git a/src/physics/RPA_tk.jl b/src/physics/RPA_tk.jl
index 03fce9c..2a33427 100644
--- a/src/physics/RPA_tk.jl
+++ b/src/physics/RPA_tk.jl
@@ -194,34 +194,6 @@ function interleave_mpo(target_mpo, phys_sites, n)
return new_mpo
end
-# ============================================================
-# Diagonal extraction
-# ============================================================
-
-"""
- extract_diagonal_to_mps(M) -> MPS
-
-Extract the diagonal of an MPO `M` as an MPS by projecting each site
-tensor onto the subspace where bra and ket indices are equal.
-"""
-function extract_diagonal_to_mps(M::MPO)::MPS
- N = length(M)
- new_tensors = Vector{ITensor}(undef, N)
- for i in 1:N
- t = M[i]
- s2, s1 = siteinds(M, i) # s2 = bra, s1 = ket
- dim_s = dim(s1)
- v_inds = uniqueinds(t, s1, s2)
- res = ITensor(v_inds..., s1)
- for v in 1:dim_s
- slice = t * onehot(s1 => v) * onehot(s2 => v)
- res += slice * onehot(s1 => v)
- end
- new_tensors[i] = res
- end
- return MPS(new_tensors)
-end
-
# ============================================================
# MPO/MPS merging utilities
# ============================================================
diff --git a/src/physics/SCF_tk.jl b/src/physics/SCF_tk.jl
index 85ec809..82de4a8 100644
--- a/src/physics/SCF_tk.jl
+++ b/src/physics/SCF_tk.jl
@@ -1205,6 +1205,7 @@ require an explicit coupling constant.
function get_scf(H0::TBHamiltonian, channel::Symbol;
interaction::Symbol = :dense,
kwargs...)
+ _require_binary_position_space(H0, "get_scf")
ch = _canonical_channel(channel)
ch === :swave &&
error("get_scf(H0, :swave) requires an explicit coupling: use get_scf(H0, g, :swave).")
@@ -1256,6 +1257,7 @@ function get_scf(H0::TBHamiltonian, U, channel::Symbol;
stop_on_increase::Bool = false,
verbose::Bool = true,
builder_kwargs...)
+ _require_binary_position_space(H0, "get_scf")
ch = _canonical_channel(channel)
dmethod = _canonical_density_method(method === nothing ? density_method : method)
diff --git a/src/physics/Topology_tk.jl b/src/physics/Topology_tk.jl
index 24cc6dc..734c3c9 100644
--- a/src/physics/Topology_tk.jl
+++ b/src/physics/Topology_tk.jl
@@ -164,6 +164,7 @@ function get_W(H::TBHamiltonian, xfunc=nothing;
quenched::Bool = true,
l = nothing,
Λ::Real = 10)
+ _require_binary_position_space(H, "get_W")
H.sublattice_s === nothing || dim(H.sublattice_s) == 2 ||
error("get_W requires a 2-component sublattice index (dim=2); got dim=$(dim(H.sublattice_s)).")
H.sublattice_s !== nothing ||
@@ -635,6 +636,7 @@ function get_C(H::TBHamiltonian, xfunc=nothing, yfunc=nothing;
Nel = nothing,
quenched::Bool = true,
sequential::Bool = false)
+ _require_binary_position_space(H, "get_C")
if xfunc === nothing || yfunc === nothing
geom = H.geometry_uc !== nothing ? H.geometry_uc :
H.geometry !== nothing ? H.geometry :
@@ -673,6 +675,7 @@ index. The returned MPO shares the same site indices as `H.mpo`.
function get_valley_operator(H::TBHamiltonian;
maxdim::Int = 500,
cutoff::Float64 = 1e-8)
+ _require_binary_position_space(H, "get_valley_operator")
H.Lx !== nothing ||
error("get_valley_operator requires a 2D Hamiltonian (H.Lx must be set).")
H.sublattice_s !== nothing ||
@@ -769,6 +772,7 @@ function get_valley_C(H::TBHamiltonian,
Nel = nothing,
quenched::Bool = true,
sequential::Bool = false)
+ _require_binary_position_space(H, "get_valley_C")
valley in (:K, :K_prime) ||
error("valley must be :K or :K_prime, got :$valley")
diff --git a/src/physics/TwoParticle_tk.jl b/src/physics/TwoParticle_tk.jl
index fed02c7..799d0ed 100644
--- a/src/physics/TwoParticle_tk.jl
+++ b/src/physics/TwoParticle_tk.jl
@@ -79,6 +79,8 @@ function exciton_hamiltonian(H_c::TBHamiltonian, H_v::TBHamiltonian, Ufunc;
tol = 1e-8,
cutoff = 1e-8,
maxdim = 200)
+ _require_binary_position_space(H_c, "exciton_hamiltonian")
+ _require_binary_position_space(H_v, "exciton_hamiltonian")
H_exc_mpo = Exciton_Hamiltonian(H_c, H_v, Ufunc;
on_site = on_site,
tol_quantics = tol_quantics,
diff --git a/src/position_spaces/Fibonacci.jl b/src/position_spaces/Fibonacci.jl
new file mode 100644
index 0000000..f186926
--- /dev/null
+++ b/src/position_spaces/Fibonacci.jl
@@ -0,0 +1,505 @@
+# Fibonacci.jl — Fibonacci position space and Hamiltonian constructors
+#
+# Fibonacci chains use Zeckendorf strings (no adjacent ones) inside an ambient
+# 2^L qubit register. The validity projector is therefore the physical identity
+# for every projected-space solver operation.
+
+"""
+ FibonacciPositionSpace(projector)
+
+Zeckendorf-encoded Fibonacci position space. `projector` is the identity on the
+`F_(L+2)` valid strings embedded in the ambient `2^L` qubit register.
+"""
+struct FibonacciPositionSpace <: AbstractPositionSpace
+ projector::MPO
+end
+
+ambient_dimension(::FibonacciPositionSpace, H::TBHamiltonian) = big(2)^H.L
+
+"""Return the `n`th Fibonacci number with `F_0=0` and `F_1=1`."""
+function fibonacci_number(n::Integer)
+ n >= 0 || throw(ArgumentError("n must be non-negative"))
+ a, b = big(0), big(1)
+ for _ in 1:n
+ a, b = b, a + b
+ end
+ return a
+end
+
+"""Number `F_(L+2)` of valid length-`L` Zeckendorf strings."""
+fibonacci_site_count(L::Integer) = Int(fibonacci_number(L + 2))
+
+"""
+ fibonacci_zeckendorf_digits(n, L) -> Vector{Int}
+
+Length-`L`, most-significant-first Zeckendorf representation of the physical
+site label `n in 0:F_(L+2)-1`.
+"""
+function fibonacci_zeckendorf_digits(n::Integer, L::Integer)
+ L >= 1 || throw(ArgumentError("L must be positive"))
+ N = fibonacci_site_count(L)
+ 0 <= n < N || throw(ArgumentError("site label must satisfy 0 <= n < $N"))
+ digits = zeros(Int, L)
+ remainder = big(n)
+ for (position, k) in enumerate((L + 1):-1:2)
+ weight = fibonacci_number(k)
+ if weight <= remainder
+ digits[position] = 1
+ remainder -= weight
+ end
+ end
+ iszero(remainder) || error("Zeckendorf conversion failed for n=$n, L=$L")
+ return digits
+end
+
+function _fibonacci_automaton_mps(sites; A=0.0, B=1.0)
+ L = length(sites)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ T = promote_type(Float64, typeof(A), typeof(B))
+ links = [Index(2, "FibAutomaton,Link,l=$i") for i in 1:(L - 1)]
+ word = MPS(sites)
+
+ first = ITensor(T, sites[1], links[1])
+ first[sites[1] => 1, links[1] => 1] = one(T)
+ first[sites[1] => 2, links[1] => 2] = one(T)
+ word[1] = first
+
+ for i in 2:(L - 1)
+ bulk = ITensor(T, links[i - 1], sites[i], links[i])
+ bulk[links[i - 1] => 1, sites[i] => 1, links[i] => 1] = one(T)
+ bulk[links[i - 1] => 2, sites[i] => 1, links[i] => 1] = one(T)
+ bulk[links[i - 1] => 1, sites[i] => 2, links[i] => 2] = one(T)
+ word[i] = bulk
+ end
+
+ last = ITensor(T, links[end], sites[end])
+ last[links[end] => 1, sites[end] => 1] = A
+ last[links[end] => 2, sites[end] => 1] = A
+ last[links[end] => 1, sites[end] => 2] = B
+ word[end] = last
+ return word
+end
+
+ITensors.op(::OpName"FibLower", ::SiteType"Qubit") = [0 1; 0 0]
+ITensors.op(::OpName"FibRaise", ::SiteType"Qubit") = [0 0; 1 0]
+ITensors.op(::OpName"FibP0", ::SiteType"Qubit") = [1 0; 0 0]
+
+"""
+ fibonacci_decrement_mpo(sites; boundary=:open) -> MPO
+
+Physical decrement `K|n> = |n-1>` in the Zeckendorf basis. With periodic
+boundaries the only added automaton transition is `|0> -> |F_(L+2)-1>`;
+wrapping never occurs at the ambient binary state `2^L-1`.
+"""
+function fibonacci_decrement_mpo(sites; boundary::Symbol=:open)
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ L = length(sites)
+ shifts = OpSum()
+ for i in 1:L
+ term = OpSum()
+ term += 1.0, "FibLower", i
+ for j in (i + 1):L
+ term *= (isodd(j - i) ? "FibRaise" : "FibP0", j)
+ end
+ shifts += term
+ end
+
+ if boundary === :periodic
+ last_digits = fibonacci_zeckendorf_digits(fibonacci_site_count(L) - 1, L)
+ wrap = OpSum()
+ first_op = isone(last_digits[1]) ? "FibRaise" : "FibP0"
+ wrap += 1.0, first_op, 1
+ for i in 2:L
+ op = isone(last_digits[i]) ? "FibRaise" : "FibP0"
+ wrap *= (op, i)
+ end
+ shifts += wrap
+ end
+ return MPO(shifts, sites)
+end
+
+_fibonacci_mpo_adjoint(A::MPO) = swapprime(dag(A), 0, 1)
+
+function physical_projector(space::FibonacciPositionSpace, H::TBHamiltonian)
+ length(H.sites) == H.L ||
+ error("FibonacciPositionSpace currently supports position-only Hamiltonians")
+ return copy(space.projector)
+end
+
+function physical_site_state(::FibonacciPositionSpace, H::TBHamiltonian, x::Integer)
+ 1 <= x <= H.N || throw(BoundsError(1:H.N, x))
+ length(H.sites) == H.L ||
+ error("FibonacciPositionSpace currently supports position-only Hamiltonians")
+ return MPS(H.sites, string.(fibonacci_zeckendorf_digits(x - 1, H.L)))
+end
+
+"""
+ fibonacci_bond_symbol(L, bond) -> Symbol
+
+Return `:A` or `:B` for the 1-indexed bond beginning at `bond` in the
+canonical `F_(L+2)`-bond periodic approximant. Bond `N` joins site `N` to
+site `1` when periodic boundaries are used.
+"""
+function fibonacci_bond_symbol(L::Integer, bond::Integer)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ N = fibonacci_site_count(L)
+ 1 <= bond <= N || throw(BoundsError(1:N, bond))
+ return iszero(fibonacci_zeckendorf_digits(bond - 1, L)[end]) ? :A : :B
+end
+
+"""
+ fibonacci_site_environment(L, site; boundary=:periodic) -> Symbol
+
+Classify a site from its adjacent bonds. `:atomic` means `AA`, while
+`:molecular_AB` and `:molecular_BA` retain the orientation of the molecular
+site. Open-chain endpoints return `:boundary`.
+"""
+function fibonacci_site_environment(L::Integer, site::Integer;
+ boundary::Symbol=:periodic)
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ N = fibonacci_site_count(L)
+ 1 <= site <= N || throw(BoundsError(1:N, site))
+ boundary === :open && site in (1, N) && return :boundary
+ left = fibonacci_bond_symbol(L, site == 1 ? N : site - 1)
+ right = fibonacci_bond_symbol(L, site)
+ left === :A && right === :A && return :atomic
+ left === :A && right === :B && return :molecular_AB
+ left === :B && right === :A && return :molecular_BA
+ error("invalid Fibonacci bond environment $left$right at site $site")
+end
+
+function _fibonacci_conumber_rank(L::Integer, site::Integer;
+ orientation::Symbol=:standard,
+ alignment::Symbol=:atomic,
+ origin::Integer=0)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ orientation in (:standard, :reversed) ||
+ throw(ArgumentError("orientation must be :standard or :reversed"))
+ alignment in (:atomic, :raw) ||
+ throw(ArgumentError("alignment must be :atomic or :raw"))
+ N = fibonacci_site_count(L)
+ 1 <= site <= N || throw(BoundsError(1:N, site))
+ multiplier = Int(fibonacci_number(L)) # F_(n-2), with N=F_n
+ x = site - 1
+ rank = Int(mod((big(x) + origin) * multiplier, N))
+
+ # For the canonical bond word used by the MPO, the AA acceptance window is
+ # [0,F_(L-1)-1] for even L and [1,F_(L-1)] for odd L in raw standard
+ # conumbers. Move it between the two F_L molecular windows. This is a
+ # cyclic cut of perpendicular space, not a change of conumber multiplier.
+ if alignment === :atomic
+ rank = Int(mod(big(rank) + multiplier - Int(isodd(L)), N))
+ end
+ orientation === :reversed && (rank = N - 1 - rank)
+ return rank
+end
+
+"""
+ fibonacci_conumber(L, site; orientation=:standard, alignment=:atomic,
+ centered=true, origin=0) -> Int
+
+Conumber of a single 1-indexed physical site, evaluated without constructing a
+Hamiltonian or allocating an `F_(L+2)`-element permutation.
+
+`alignment=:atomic` makes the conumber cut compatible with TensorBinding's
+canonical Fibonacci bond phase: the `AA` sites form one central block, between
+the `AB` and `BA` molecular blocks. `alignment=:raw` exposes the unshifted
+modular residue. `centered=true` labels the ordered ranks around zero; it does
+not perform an additional cyclic permutation.
+"""
+function fibonacci_conumber(L::Integer, site::Integer;
+ orientation::Symbol=:standard,
+ alignment::Symbol=:atomic,
+ centered::Bool=true,
+ origin::Integer=0)
+ N = fibonacci_site_count(L)
+ rank = _fibonacci_conumber_rank(
+ L, site; orientation, alignment, origin,
+ )
+ return centered ? rank - fld(N, 2) : rank
+end
+
+"""Inverse of [`fibonacci_conumber`](@ref), returning a 1-indexed site."""
+function fibonacci_site_from_conumber(L::Integer, conumber::Integer;
+ orientation::Symbol=:standard,
+ alignment::Symbol=:atomic,
+ centered::Bool=true,
+ origin::Integer=0)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ orientation in (:standard, :reversed) ||
+ throw(ArgumentError("orientation must be :standard or :reversed"))
+ alignment in (:atomic, :raw) ||
+ throw(ArgumentError("alignment must be :atomic or :raw"))
+ N = fibonacci_site_count(L)
+ rank = centered ? conumber + fld(N, 2) : conumber
+ 0 <= rank < N || throw(BoundsError(0:(N - 1), rank))
+ orientation === :reversed && (rank = N - 1 - rank)
+ multiplier = Int(fibonacci_number(L))
+ if alignment === :atomic
+ rank = mod(rank - multiplier + Int(isodd(L)), N)
+ end
+ x = Int(mod(big(rank) * invmod(multiplier, N) - origin, N))
+ return x + 1
+end
+
+"""
+ fibonacci_rg_partition(L; depth=0, centered=true)
+
+Return the molecular–atomic–molecular conumber intervals after `depth`
+successive atomic deflations. Each deflation maps `L -> L-3`. The returned
+ranges are embedded in the original conumber ordering, so a zoom should slice
+these ranges directly rather than re-conumbering the selected sites as a new
+canonical chain.
+"""
+function fibonacci_rg_partition(L::Integer; depth::Integer=0,
+ centered::Bool=true)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ depth >= 0 || throw(ArgumentError("depth must be non-negative"))
+ effective_L = Int(L)
+ window_first = 0
+ for _ in 1:depth
+ effective_L >= 5 ||
+ throw(ArgumentError("depth=$depth deflates L=$L below the supported L=2 approximant"))
+ window_first += Int(fibonacci_number(effective_L))
+ effective_L -= 3
+ end
+
+ molecular_count = Int(fibonacci_number(effective_L))
+ atomic_count = Int(fibonacci_number(effective_L - 1))
+ window_count = fibonacci_site_count(effective_L)
+ @assert 2molecular_count + atomic_count == window_count
+
+ left = window_first:(window_first + molecular_count - 1)
+ atomic = (last(left) + 1):(last(left) + atomic_count)
+ right = (last(atomic) + 1):(window_first + window_count - 1)
+ window = window_first:(window_first + window_count - 1)
+ shift = centered ? fld(fibonacci_site_count(L), 2) : 0
+ shift_range(r) = (first(r) - shift):(last(r) - shift)
+ return (;
+ depth, effective_L, window_count, molecular_count, atomic_count,
+ window_ranks=window, left_molecular_ranks=left,
+ atomic_ranks=atomic, right_molecular_ranks=right,
+ window_axis=shift_range(window),
+ left_molecular_axis=shift_range(left),
+ atomic_axis=shift_range(atomic),
+ right_molecular_axis=shift_range(right),
+ )
+end
+
+"""
+ fibonacci_atomic_depth(L, site; kwargs...) -> Int
+
+Number of consecutive atomic deflations containing `site`. This directly
+tests whether a site remains inside the nested central atomic windows.
+"""
+function fibonacci_atomic_depth(L::Integer, site::Integer;
+ orientation::Symbol=:standard,
+ origin::Integer=0)
+ rank = _fibonacci_conumber_rank(
+ L, site; orientation, alignment=:atomic, origin,
+ )
+ effective_L = Int(L)
+ window_first = 0
+ depth = 0
+ while effective_L >= 2
+ molecular_count = Int(fibonacci_number(effective_L))
+ atomic_count = Int(fibonacci_number(effective_L - 1))
+ atomic_first = window_first + molecular_count
+ atomic_last = atomic_first + atomic_count - 1
+ atomic_first <= rank <= atomic_last || break
+ depth += 1
+ window_first = atomic_first
+ effective_L -= 3
+ end
+ return depth
+end
+
+function _fibonacci_conumbering(H::TBHamiltonian;
+ orientation::Symbol=:standard,
+ centered::Bool=true,
+ origin::Integer=0,
+ alignment::Symbol=:atomic)
+ H.position_space isa FibonacciPositionSpace ||
+ throw(ArgumentError("conumbering requires FibonacciPositionSpace"))
+ N = H.N
+ ranks = [_fibonacci_conumber_rank(
+ H.L, site; orientation, alignment, origin,
+ ) for site in 1:N]
+ labels = centered ? ranks .- fld(N, 2) : ranks
+ permutation = sortperm(labels)
+ axis = labels[permutation]
+ @assert length(unique(axis)) == N
+ multiplier = orientation === :standard ?
+ Int(fibonacci_number(H.L)) : Int(fibonacci_number(H.L + 1))
+ return (; axis, permutation, ranks, labels, multiplier,
+ orientation, centered, origin, alignment)
+end
+
+function site_axis(::FibonacciPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical,
+ orientation::Symbol=:standard,
+ centered::Bool=true,
+ origin::Integer=0,
+ alignment::Symbol=:atomic,
+ kwargs...)
+ ordering === :physical && return collect(0:(H.N - 1))
+ ordering === :conumber ||
+ throw(ArgumentError("ordering must be :physical or :conumber"))
+ return _fibonacci_conumbering(
+ H; orientation, centered, origin, alignment,
+ ).axis
+end
+
+function site_permutation(::FibonacciPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical,
+ orientation::Symbol=:standard,
+ centered::Bool=true,
+ origin::Integer=0,
+ alignment::Symbol=:atomic,
+ kwargs...)
+ ordering === :physical && return collect(1:H.N)
+ ordering === :conumber ||
+ throw(ArgumentError("ordering must be :physical or :conumber"))
+ return _fibonacci_conumbering(
+ H; orientation, centered, origin, alignment,
+ ).permutation
+end
+
+# Convenience overload of the sampling utility defined in core/Utils.jl; it
+# lives here because its signature needs TBHamiltonian and FibonacciPositionSpace.
+function fibonacci_ldos_sampling_plan(H::TBHamiltonian; kwargs...)
+ H.position_space isa FibonacciPositionSpace ||
+ throw(ArgumentError("fibonacci LDOS sampling requires FibonacciPositionSpace"))
+ expected_N = fibonacci_site_count(H.L)
+ H.N == expected_N ||
+ throw(ArgumentError("Hamiltonian has N=$(H.N), expected F_(L+2)=$expected_N for L=$(H.L)"))
+ return fibonacci_ldos_sampling_plan(H.L; kwargs...)
+end
+
+"""
+ fibonacci_hamiltonian(L; A, B, model=:hopping, t=1.0, onsite=0.0,
+ boundary=:periodic, scale=nothing, padding=1.05,
+ cutoff=1e-12, maxdim=200) -> TBHamiltonian
+
+Construct a Fibonacci chain in the projected Zeckendorf position space.
+
+- `model=:onsite`: `A` and `B` are onsite energies and `t` is uniform hopping.
+- `model=:hopping`: `A` and `B` are bond amplitudes and `onsite` is uniform.
+
+The default periodic boundary closes the physical `F_(L+2)`-site approximant.
+For odd `F_(L+2)`, a periodic hopping chain is an odd cycle and therefore is not
+exactly chiral even when `onsite=0`.
+"""
+function fibonacci_hamiltonian(
+ L::Integer; A, B,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+ scale=nothing,
+ padding::Real=1.05,
+ cutoff::Real=1e-12,
+ maxdim::Integer=200,
+)
+ L >= 2 || throw(ArgumentError("Fibonacci chains require L >= 2"))
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ padding > 1 || throw(ArgumentError("padding must be greater than 1"))
+ if model === :onsite
+ isreal(A) && isreal(B) ||
+ throw(ArgumentError("onsite Fibonacci values A and B must be real"))
+ else
+ isreal(onsite) ||
+ throw(ArgumentError("the uniform onsite energy must be real"))
+ end
+
+ sites = siteinds("Qubit", L; conserve_qns=false)
+ word_mps = _fibonacci_automaton_mps(sites; A, B)
+ valid_mps = _fibonacci_automaton_mps(sites; A=1.0, B=1.0)
+ word = mps_to_diagonal_mpo(word_mps, sites)
+ P = mps_to_diagonal_mpo(valid_mps, sites)
+ K = fibonacci_decrement_mpo(sites; boundary)
+
+ V, TK = if model === :onsite
+ word, t * K
+ else
+ onsite * P, apply(word, K; cutoff=cutoff, maxdim=maxdim)
+ end
+ hopping = +(TK, _fibonacci_mpo_adjoint(TK); cutoff=cutoff, maxdim=maxdim)
+ Hraw = +(V, hopping; cutoff=cutoff, maxdim=maxdim)
+ mpo = apply(P, apply(Hraw, P; cutoff=cutoff, maxdim=maxdim);
+ cutoff=cutoff, maxdim=maxdim)
+ ITensorMPS.truncate!(mpo; cutoff=cutoff, maxdim=maxdim)
+
+ center, halfwidth = if model === :onsite
+ lo, hi = extrema((Float64(real(A)), Float64(real(B))))
+ ((lo + hi) / 2, (hi - lo) / 2 + 2abs(t))
+ else
+ (Float64(real(onsite)), 2max(abs(A), abs(B)))
+ end
+ scale_value = isnothing(scale) ? padding * Float64(halfwidth) : Float64(scale)
+ scale_value > 0 || throw(ArgumentError("KPM scale must be positive"))
+
+ N = fibonacci_site_count(L)
+ H = TBHamiltonian(L, N, sites, mpo, _chain_geometry(),
+ scale_value, Float64(center),
+ nothing, nothing, nothing, nothing, 0, nothing)
+ H.position_space = FibonacciPositionSpace(P)
+ return H
+end
+
+function _build_fibonacci(params, L::Integer;
+ scale=nothing, tol=1e-12, maxdim=200, kwargs...)
+ p = if params isa NamedTuple
+ Dict{Symbol,Any}(pairs(params))
+ elseif params isa AbstractDict
+ Dict{Symbol,Any}(Symbol(k) => v for (k, v) in pairs(params))
+ else
+ throw(ArgumentError("fibonacci parameters must be a NamedTuple or dictionary containing A and B"))
+ end
+ haskey(p, :A) && haskey(p, :B) ||
+ throw(ArgumentError("fibonacci parameters must contain A and B"))
+ allowed = Set((:A, :B, :t, :onsite))
+ unknown = setdiff(Set(keys(p)), allowed)
+ isempty(unknown) || throw(ArgumentError("unknown fibonacci parameters: $(collect(unknown))"))
+ return fibonacci_hamiltonian(
+ L; A=p[:A], B=p[:B],
+ t=get(p, :t, 1.0), onsite=get(p, :onsite, 0.0),
+ scale=scale, cutoff=tol, maxdim=maxdim, kwargs...,
+ )
+end
+
+# Dense small-system oracle used only by the test suite.
+function _dense_fibonacci_hamiltonian(
+ L::Integer; A, B,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+)
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ N = fibonacci_site_count(L)
+ word = [iszero(fibonacci_zeckendorf_digits(n, L)[end]) ? A : B
+ for n in 0:(N - 1)]
+ diagonal = model === :onsite ? word : fill(onsite, N)
+ bonds = model === :onsite ? fill(t, N - 1) : word[1:(N - 1)]
+ H = zeros(ComplexF64, N, N)
+ H[diagind(H)] .= diagonal
+ for n in 1:(N - 1)
+ H[n, n + 1] = bonds[n]
+ H[n + 1, n] = conj(bonds[n])
+ end
+ if boundary === :periodic
+ wrap = model === :onsite ? t : word[end]
+ H[N, 1] = wrap
+ H[1, N] = conj(wrap)
+ end
+ return H
+end
diff --git a/src/position_spaces/KBonacci.jl b/src/position_spaces/KBonacci.jl
new file mode 100644
index 0000000..3aa0d5d
--- /dev/null
+++ b/src/position_spaces/KBonacci.jl
@@ -0,0 +1,380 @@
+# KBonacci.jl — k-bonacci quasicrystal position spaces
+#
+# The k-bonacci word on the alphabet a_1, …, a_k (written A, B, C, …) is the
+# fixed point of the substitution a_i -> a_1 a_(i+1) for i < k and a_k -> a_1:
+# k = 2 is Fibonacci (A -> AB, B -> A), k = 3 Tribonacci (A -> AB, B -> AC,
+# C -> A), k = 4 Tetranacci, and so on. Sites are labelled by their binary
+# expansion in the k-bonacci numeration system with weights
+# w_l = 2^l for l < k, w_l = w_(l-1) + w_(l-2) + … + w_(l-k) for l >= k,
+# and the admissibility rule "no k consecutive ones" (k = 2: Zeckendorf).
+# w_l counts the admissible length-l strings, so L digits enumerate exactly
+# w_L physical sites inside the ambient 2^L Qubit register (w_L = F_(L+2) for
+# k = 2 and the Tribonacci number T_(L+3) for k = 3). The letter at site n is
+# a_(r+1), where r in 0:k-1 is the number of trailing ones of n, so the word and
+# the validity indicator are k-state automaton MPS (bond dimension k) and the
+# Hamiltonian MPO is exact at any L. For k = 2 the construction coincides with
+# Fibonacci.jl on the same Qubit sites; the qubit digit operators FibLower,
+# FibRaise and FibP0 defined there are reused here.
+
+"""
+ KBonacciPositionSpace(k, projector)
+
+Projected position space of the k-bonacci chain of order `k`. `projector` is
+the identity on the `w_L` admissible binary strings (no `k` consecutive ones)
+embedded in the ambient `2^L` Qubit register.
+"""
+struct KBonacciPositionSpace <: AbstractPositionSpace
+ k::Int
+ projector::MPO
+end
+
+ambient_dimension(::KBonacciPositionSpace, H::TBHamiltonian) = big(2)^H.L
+
+"""
+ kbonacci_number(k, n) -> BigInt
+
+Return `w_n`, the number of length-`n` binary strings without `k` consecutive
+ones: `w_n = 2^n` for `n < k` and `w_n = w_(n-1) + … + w_(n-k)` otherwise. In the
+usual seeding of the k-bonacci sequence (`k - 1` zeros followed by a one) this is
+its `(n + k)`-th term: `F_(n+2)` for `k = 2`, the Tribonacci number `T_(n+3)` for
+`k = 3` (`1, 2, 4, 7, 13, 24, …`), the Tetranacci number for `k = 4`
+(`1, 2, 4, 8, 15, 29, …`).
+"""
+function kbonacci_number(k::Integer, n::Integer)
+ k >= 2 || throw(ArgumentError("k-bonacci order k must be at least 2"))
+ n >= 0 || throw(ArgumentError("n must be non-negative"))
+ n < k && return big(2)^n
+ window = [big(2)^l for l in 0:(k - 1)] # w_0, …, w_(k-1)
+ for _ in k:n
+ push!(window, sum(window))
+ popfirst!(window)
+ end
+ return window[end]
+end
+
+"""Number `w_L` of admissible length-`L` binary strings of the k-bonacci chain."""
+kbonacci_site_count(k::Integer, L::Integer) = Int(kbonacci_number(k, L))
+
+"""
+ kbonacci_digits(k, n, L) -> Vector{Int}
+
+Length-`L`, most-significant-first binary expansion of the physical site label
+`n in 0:w_L-1` in the k-bonacci numeration system, with weights
+`w_(L-1), …, w_0`, obtained greedily. The result never contains `k`
+consecutive ones. For `k = 2` this is [`fibonacci_zeckendorf_digits`](@ref).
+"""
+function kbonacci_digits(k::Integer, n::Integer, L::Integer)
+ L >= 1 || throw(ArgumentError("L must be positive"))
+ N = kbonacci_site_count(k, L)
+ 0 <= n < N || throw(ArgumentError("site label must satisfy 0 <= n < $N"))
+ digits = zeros(Int, L)
+ remainder = big(n)
+ for position in 1:L
+ weight = kbonacci_number(k, L - position)
+ if remainder >= weight
+ digits[position] = 1
+ remainder -= weight
+ end
+ end
+ iszero(remainder) || error("k-bonacci conversion failed for k=$k, n=$n, L=$L")
+ return digits
+end
+
+# Number of trailing ones of a digit string; this selects the letter.
+function _kbonacci_trailing_ones(digits::AbstractVector{<:Integer})
+ r = 0
+ for d in Iterators.reverse(digits)
+ d == 1 || break
+ r += 1
+ end
+ return r
+end
+
+"""Letter symbols `[:A, :B, …]` of the k-bonacci alphabet (`k <= 26`)."""
+function kbonacci_letters(k::Integer)
+ 2 <= k <= 26 || throw(ArgumentError("letter symbols are defined for 2 <= k <= 26"))
+ return [Symbol('A' + i) for i in 0:(k - 1)]
+end
+
+# k-state automaton MPS on Qubit sites. Link state s (1-based) means "the string
+# read so far ends in s-1 ones" (state 1 after a 0 or before reading anything).
+# A digit 0 resets to state 1, a digit 1 advances s -> s+1 and is forbidden from
+# state k, which annihilates every string with k consecutive ones. The amplitude
+# of an admissible string is values[r+1] with r its number of trailing ones.
+function _kbonacci_automaton_mps(k::Integer, sites; values)
+ L = length(sites)
+ k >= 2 || throw(ArgumentError("k-bonacci order k must be at least 2"))
+ L >= 2 || throw(ArgumentError("k-bonacci chains require L >= 2"))
+ length(values) == k ||
+ throw(ArgumentError("values must hold one amplitude per letter (k = $k), got $(length(values))"))
+ all(s -> dim(s) == 2, sites) ||
+ throw(ArgumentError("k-bonacci sites must be binary (Qubit) sites"))
+ T = promote_type(Float64, map(typeof, Tuple(values))...)
+ links = [Index(k, "KBAutomaton,Link,l=$i") for i in 1:(L - 1)]
+ word = MPS(sites)
+
+ first = ITensor(T, sites[1], links[1])
+ first[sites[1] => 1, links[1] => 1] = one(T)
+ first[sites[1] => 2, links[1] => 2] = one(T)
+ word[1] = first
+
+ for i in 2:(L - 1)
+ bulk = ITensor(T, links[i - 1], sites[i], links[i])
+ for s in 1:k
+ bulk[links[i - 1] => s, sites[i] => 1, links[i] => 1] = one(T)
+ end
+ for s in 1:(k - 1)
+ bulk[links[i - 1] => s, sites[i] => 2, links[i] => s + 1] = one(T)
+ end
+ word[i] = bulk
+ end
+
+ last = ITensor(T, links[end], sites[end])
+ for s in 1:k
+ last[links[end] => s, sites[end] => 1] = T(values[1])
+ end
+ for s in 1:(k - 1)
+ last[links[end] => s, sites[end] => 2] = T(values[s + 1])
+ end
+ word[end] = last
+ return word
+end
+
+"""
+ kbonacci_decrement_mpo(k, sites; boundary=:open) -> MPO
+
+Physical decrement `K|n> = |n-1>` in the k-bonacci numeration system. The
+`i`-th term fires when the least significant one sits at position `i`: that
+digit is cleared and the (all-zero) tail is rewritten with the repeating
+pattern `1, …, 1, 0` (`k - 1` ones then a zero), which encodes `w_j - 1` on `j`
+digits. With periodic boundaries the only added transition is
+`|0> -> |w_L - 1>`; the ambient register state `2^L - 1` is never wrapped.
+Uses the qubit operators `FibLower`, `FibRaise`, `FibP0` of Fibonacci.jl.
+"""
+function kbonacci_decrement_mpo(k::Integer, sites; boundary::Symbol=:open)
+ k >= 2 || throw(ArgumentError("k-bonacci order k must be at least 2"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ L = length(sites)
+ shifts = OpSum()
+ for i in 1:L
+ term = OpSum()
+ term += 1.0, "FibLower", i
+ for j in (i + 1):L
+ term *= ((j - i) % k == 0 ? "FibP0" : "FibRaise", j)
+ end
+ shifts += term
+ end
+
+ if boundary === :periodic
+ last_digits = kbonacci_digits(k, kbonacci_site_count(k, L) - 1, L)
+ op_for(digit) = isone(digit) ? "FibRaise" : "FibP0"
+ wrap = OpSum()
+ wrap += 1.0, op_for(last_digits[1]), 1
+ for i in 2:L
+ wrap *= (op_for(last_digits[i]), i)
+ end
+ shifts += wrap
+ end
+ return MPO(shifts, sites)
+end
+
+function physical_projector(space::KBonacciPositionSpace, H::TBHamiltonian)
+ length(H.sites) == H.L ||
+ error("KBonacciPositionSpace currently supports position-only Hamiltonians")
+ return copy(space.projector)
+end
+
+function physical_site_state(space::KBonacciPositionSpace, H::TBHamiltonian, x::Integer)
+ 1 <= x <= H.N || throw(BoundsError(1:H.N, x))
+ length(H.sites) == H.L ||
+ error("KBonacciPositionSpace currently supports position-only Hamiltonians")
+ return MPS(H.sites, string.(kbonacci_digits(space.k, x - 1, H.L)))
+end
+
+function site_axis(::KBonacciPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical || throw(ArgumentError(
+ "ordering=:$ordering is not available for KBonacciPositionSpace; " *
+ "only :physical is defined (conumbering is currently Fibonacci-only)"))
+ return collect(0:(H.N - 1))
+end
+
+function site_permutation(::KBonacciPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical || throw(ArgumentError(
+ "ordering=:$ordering is not available for KBonacciPositionSpace; " *
+ "only :physical is defined (conumbering is currently Fibonacci-only)"))
+ return collect(1:H.N)
+end
+
+"""
+ kbonacci_bond_symbol(k, L, bond) -> Symbol
+
+Return the letter (`:A`, `:B`, `:C`, …) of the 1-indexed bond beginning at
+`bond` in the canonical `w_L`-bond periodic approximant of the k-bonacci chain:
+the letter of physical site `bond - 1`, selected by its number of trailing ones.
+Bond `N` joins site `N` to site `1` when periodic boundaries are used.
+"""
+function kbonacci_bond_symbol(k::Integer, L::Integer, bond::Integer)
+ N = kbonacci_site_count(k, L)
+ 1 <= bond <= N || throw(BoundsError(1:N, bond))
+ r = _kbonacci_trailing_ones(kbonacci_digits(k, bond - 1, L))
+ return kbonacci_letters(k)[r + 1]
+end
+
+"""
+ kbonacci_hamiltonian(k, L; values, model=:hopping, t=1.0, onsite=0.0,
+ boundary=:periodic, scale=nothing, padding=1.05,
+ cutoff=1e-12, maxdim=200) -> TBHamiltonian
+
+Construct the k-bonacci chain of order `k` (`a_i -> a_1 a_(i+1)`, `a_k -> a_1`)
+in its projected numeration position space on `L` Qubit sites. `values` holds
+one amplitude per letter `a_1, …, a_k` (i.e. `A, B, C, …`) and the chain has
+`H.N = w_L` physical sites.
+
+- `model=:onsite`: `values` are onsite energies and `t` is uniform hopping.
+- `model=:hopping`: `values` are bond amplitudes and `onsite` is uniform.
+
+The Hamiltonian is assembled as `P (V + T K + h.c.) P`, where `P` is the
+validity projector, `T`/`V` the diagonal word MPO, and `K` the decrement
+[`kbonacci_decrement_mpo`](@ref). `k = 2` reproduces the Fibonacci chain of
+[`fibonacci_hamiltonian`](@ref) with `values = (A, B)`; `k = 3` is the
+Tribonacci chain. The default periodic boundary closes the physical `w_L`-site
+approximant.
+"""
+function kbonacci_hamiltonian(
+ k::Integer, L::Integer; values,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+ scale=nothing,
+ padding::Real=1.05,
+ cutoff::Real=1e-12,
+ maxdim::Integer=200,
+)
+ k >= 2 || throw(ArgumentError("k-bonacci order k must be at least 2"))
+ L >= 2 || throw(ArgumentError("k-bonacci chains require L >= 2"))
+ length(values) == k ||
+ throw(ArgumentError("values must hold exactly k = $k letter amplitudes, got $(length(values))"))
+ all(v -> v isa Number, values) ||
+ throw(ArgumentError("values must be numbers"))
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ padding > 1 || throw(ArgumentError("padding must be greater than 1"))
+ amplitudes = collect(values)
+ if model === :onsite
+ all(isreal, amplitudes) ||
+ throw(ArgumentError("onsite k-bonacci values must be real"))
+ else
+ isreal(onsite) ||
+ throw(ArgumentError("the uniform onsite energy must be real"))
+ end
+
+ sites = siteinds("Qubit", L; conserve_qns=false)
+ word_mps = _kbonacci_automaton_mps(k, sites; values=amplitudes)
+ valid_mps = _kbonacci_automaton_mps(k, sites; values=ones(k))
+ word = mps_to_diagonal_mpo(word_mps, sites)
+ P = mps_to_diagonal_mpo(valid_mps, sites)
+ K = kbonacci_decrement_mpo(k, sites; boundary)
+
+ V, TK = if model === :onsite
+ word, t * K
+ else
+ onsite * P, apply(word, K; cutoff=cutoff, maxdim=maxdim)
+ end
+ hopping = +(TK, shift_adjoint_mpo(TK); cutoff=cutoff, maxdim=maxdim)
+ Hraw = +(V, hopping; cutoff=cutoff, maxdim=maxdim)
+ mpo = apply(P, apply(Hraw, P; cutoff=cutoff, maxdim=maxdim);
+ cutoff=cutoff, maxdim=maxdim)
+ ITensorMPS.truncate!(mpo; cutoff=cutoff, maxdim=maxdim)
+
+ center, halfwidth = if model === :onsite
+ lo, hi = extrema(Float64.(real.(amplitudes)))
+ ((lo + hi) / 2, (hi - lo) / 2 + 2abs(t))
+ else
+ (Float64(real(onsite)), 2maximum(abs.(amplitudes)))
+ end
+ scale_value = isnothing(scale) ? padding * Float64(halfwidth) : Float64(scale)
+ scale_value > 0 || throw(ArgumentError("KPM scale must be positive"))
+
+ N = kbonacci_site_count(k, L)
+ H = TBHamiltonian(L, N, sites, mpo, _chain_geometry(),
+ scale_value, Float64(center),
+ nothing, nothing, nothing, nothing, 0, nothing)
+ H.position_space = KBonacciPositionSpace(Int(k), P)
+ return H
+end
+
+function _build_kbonacci(params, L::Integer;
+ k=nothing, scale=nothing, tol=1e-12, maxdim=200,
+ kwargs...)
+ k === nothing && throw(ArgumentError(
+ "get_Hamiltonian(\"kbonacci\", …) requires the keyword k " *
+ "(k=2 Fibonacci, k=3 Tribonacci, k=4 Tetranacci, …)"))
+ k >= 2 || throw(ArgumentError("k-bonacci order k must be at least 2"))
+ p = if params isa NamedTuple
+ Dict{Symbol,Any}(pairs(params))
+ elseif params isa AbstractDict
+ Dict{Symbol,Any}(Symbol(key) => v for (key, v) in pairs(params))
+ else
+ throw(ArgumentError("k-bonacci parameters must be a NamedTuple or dictionary " *
+ "containing values=(a_1, …, a_k) or the letter keys A, B, C, …"))
+ end
+ letters = k <= 26 ? kbonacci_letters(k) : Symbol[]
+ values = if haskey(p, :values)
+ any(letter -> haskey(p, letter), letters) &&
+ throw(ArgumentError("give either values=(…) or the letter keys $(join(letters, ", ")), not both"))
+ p[:values]
+ elseif !isempty(letters) && all(letter -> haskey(p, letter), letters)
+ [p[letter] for letter in letters]
+ else
+ throw(ArgumentError("k-bonacci parameters must contain values=(a_1, …, a_k)" *
+ (isempty(letters) ? "" : " or all of the letter keys $(join(letters, ", "))")))
+ end
+ allowed = Set((:values, :t, :onsite, letters...))
+ unknown = setdiff(Set(keys(p)), allowed)
+ isempty(unknown) || throw(ArgumentError("unknown k-bonacci parameters: $(collect(unknown))"))
+ return kbonacci_hamiltonian(
+ k, L; values,
+ t=get(p, :t, 1.0), onsite=get(p, :onsite, 0.0),
+ scale=scale, cutoff=tol, maxdim=maxdim, kwargs...,
+ )
+end
+
+# Dense small-system oracle used only by the test suite.
+function _dense_kbonacci_hamiltonian(
+ k::Integer, L::Integer; values,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+)
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ length(values) == k ||
+ throw(ArgumentError("values must hold exactly k = $k letter amplitudes"))
+ N = kbonacci_site_count(k, L)
+ amplitudes = collect(values)
+ word = [amplitudes[_kbonacci_trailing_ones(kbonacci_digits(k, n, L)) + 1]
+ for n in 0:(N - 1)]
+ diagonal = model === :onsite ? word : fill(onsite, N)
+ bonds = model === :onsite ? fill(t, N - 1) : word[1:(N - 1)]
+ H = zeros(ComplexF64, N, N)
+ H[diagind(H)] .= diagonal
+ for n in 1:(N - 1)
+ H[n, n + 1] = bonds[n]
+ H[n + 1, n] = conj(bonds[n])
+ end
+ if boundary === :periodic
+ wrap = model === :onsite ? t : word[end]
+ H[N, 1] = wrap
+ H[1, N] = conj(wrap)
+ end
+ return H
+end
diff --git a/src/position_spaces/MetallicMean.jl b/src/position_spaces/MetallicMean.jl
new file mode 100644
index 0000000..790fc6a
--- /dev/null
+++ b/src/position_spaces/MetallicMean.jl
@@ -0,0 +1,350 @@
+# MetallicMean.jl — metallic-mean quasicrystal position spaces
+#
+# The metallic-mean word with parameter m is the fixed point of the substitution
+# A -> A^m B, B -> A (m = 1 Fibonacci, m = 2 silver mean / Pell, m = 3 bronze mean).
+# Sites are labelled by their expansion in the numeration system with basis
+# q_0 = 1, q_1 = m + 1, q_(l+1) = m q_l + q_(l-1),
+# digits in 0:m and the admissibility rule "a digit m must be followed by 0".
+# L digits enumerate exactly q_L physical sites inside the ambient (m+1)^L
+# register of Qudit sites. The letter at site n is B iff the least significant
+# digit of n is m, so both the word and the validity indicator are
+# bond-dimension-2 automaton MPS and the Hamiltonian MPO is exact at any L.
+# For m = 1 this reproduces the Zeckendorf construction of Fibonacci.jl on
+# dimension-2 Qudit sites.
+
+"""
+ MetallicMeanPositionSpace(m, projector)
+
+Projected position space of the metallic-mean chain with parameter `m`.
+`projector` is the identity on the `q_L` admissible digit strings embedded in
+the ambient `(m+1)^L` Qudit register.
+"""
+struct MetallicMeanPositionSpace <: AbstractPositionSpace
+ m::Int
+ projector::MPO
+end
+
+ambient_dimension(space::MetallicMeanPositionSpace, H::TBHamiltonian) =
+ big(space.m + 1)^H.L
+
+"""
+ metallic_mean_number(m, n) -> BigInt
+
+Return `q_n` for the metallic mean with parameter `m`, where `q_0 = 1`,
+`q_1 = m + 1`, and `q_(n+1) = m q_n + q_(n-1)`. For `m = 1` this is `F_(n+2)`;
+for `m = 2` it is the sequence `1, 3, 7, 17, 41, …` of the silver-mean chain.
+"""
+function metallic_mean_number(m::Integer, n::Integer)
+ m >= 1 || throw(ArgumentError("metallic-mean parameter m must be at least 1"))
+ n >= 0 || throw(ArgumentError("n must be non-negative"))
+ a, b = big(1), big(m + 1)
+ for _ in 1:n
+ a, b = b, m * b + a
+ end
+ return a
+end
+
+"""Number `q_L` of admissible length-`L` digit strings for the metallic mean `m`."""
+metallic_mean_site_count(m::Integer, L::Integer) = Int(metallic_mean_number(m, L))
+
+"""
+ metallic_mean_digits(m, n, L) -> Vector{Int}
+
+Length-`L`, most-significant-first expansion of the physical site label
+`n in 0:q_L-1` in the metallic-mean numeration system: digits in `0:m` with
+weights `q_(L-1), …, q_0`, obtained greedily. The result never contains a
+digit `m` followed by a nonzero digit.
+"""
+function metallic_mean_digits(m::Integer, n::Integer, L::Integer)
+ L >= 1 || throw(ArgumentError("L must be positive"))
+ N = metallic_mean_site_count(m, L)
+ 0 <= n < N || throw(ArgumentError("site label must satisfy 0 <= n < $N"))
+ digits = zeros(Int, L)
+ remainder = big(n)
+ for position in 1:L
+ weight = metallic_mean_number(m, L - position)
+ digit, remainder = divrem(remainder, weight)
+ digits[position] = Int(digit)
+ end
+ iszero(remainder) || error("metallic-mean conversion failed for m=$m, n=$n, L=$L")
+ return digits
+end
+
+# Bond-dimension-2 automaton MPS on (m+1)-dimensional sites. Link state 1 means
+# "the last digit read was in 0:m-1 (or nothing was read yet)", link state 2
+# means "the last digit read was m". The amplitude is zero on inadmissible
+# strings, A on admissible strings whose last digit is below m, and B on
+# admissible strings ending in m.
+function _metallic_mean_automaton_mps(m::Integer, sites; A=0.0, B=1.0)
+ L = length(sites)
+ L >= 2 || throw(ArgumentError("metallic-mean chains require L >= 2"))
+ d = m + 1
+ all(s -> dim(s) == d, sites) ||
+ throw(ArgumentError("metallic-mean sites must have local dimension m+1 = $d"))
+ T = promote_type(Float64, typeof(A), typeof(B))
+ links = [Index(2, "MMAutomaton,Link,l=$i") for i in 1:(L - 1)]
+ word = MPS(sites)
+
+ first = ITensor(T, sites[1], links[1])
+ for σ in 0:(m - 1)
+ first[sites[1] => σ + 1, links[1] => 1] = one(T)
+ end
+ first[sites[1] => d, links[1] => 2] = one(T)
+ word[1] = first
+
+ for i in 2:(L - 1)
+ bulk = ITensor(T, links[i - 1], sites[i], links[i])
+ for σ in 0:(m - 1)
+ bulk[links[i - 1] => 1, sites[i] => σ + 1, links[i] => 1] = one(T)
+ end
+ bulk[links[i - 1] => 1, sites[i] => d, links[i] => 2] = one(T)
+ bulk[links[i - 1] => 2, sites[i] => 1, links[i] => 1] = one(T) # after m only 0
+ word[i] = bulk
+ end
+
+ last = ITensor(T, links[end], sites[end])
+ for σ in 0:(m - 1)
+ last[links[end] => 1, sites[end] => σ + 1] = A
+ end
+ last[links[end] => 1, sites[end] => d] = B
+ last[links[end] => 2, sites[end] => 1] = A
+ word[end] = last
+ return word
+end
+
+# Local digit operators on a Qudit of dimension d = m + 1 (1-based matrix
+# index = digit + 1). ITensors dispatches Qudit operators on the site dimension.
+# MMLower = Σ_{k=1}^{m} |k-1><0| projector on the digit 0
+# MMRaiseTop = |m><0| turns a 0 into the top digit m
+function ITensors.op(::OpName"MMLower", ::SiteType"Qudit", d::Int)
+ mat = zeros(Float64, d, d)
+ for k in 1:(d - 1)
+ mat[k, k + 1] = 1.0
+ end
+ return mat
+end
+function ITensors.op(::OpName"MMP0", ::SiteType"Qudit", d::Int)
+ mat = zeros(Float64, d, d)
+ mat[1, 1] = 1.0
+ return mat
+end
+function ITensors.op(::OpName"MMRaiseTop", ::SiteType"Qudit", d::Int)
+ mat = zeros(Float64, d, d)
+ mat[d, 1] = 1.0
+ return mat
+end
+
+"""
+ metallic_mean_decrement_mpo(m, sites; boundary=:open) -> MPO
+
+Physical decrement `K|n> = |n-1>` in the metallic-mean numeration system. The
+`i`-th term fires when the least significant nonzero digit sits at position `i`:
+that digit is lowered by one and the (all-zero) tail becomes `m, 0, m, 0, …`,
+which encodes `q_k - 1`. With periodic boundaries the only added transition is
+`|0> -> |q_L - 1>`; wrapping never occurs at an inadmissible register state.
+"""
+function metallic_mean_decrement_mpo(m::Integer, sites; boundary::Symbol=:open)
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ L = length(sites)
+ shifts = OpSum()
+ for i in 1:L
+ term = OpSum()
+ term += 1.0, "MMLower", i
+ for j in (i + 1):L
+ term *= (isodd(j - i) ? "MMRaiseTop" : "MMP0", j)
+ end
+ shifts += term
+ end
+
+ if boundary === :periodic
+ last_digits = metallic_mean_digits(m, metallic_mean_site_count(m, L) - 1, L)
+ op_for(digit) = digit == 0 ? "MMP0" :
+ digit == m ? "MMRaiseTop" :
+ error("unexpected digit $digit in the expansion of q_L - 1")
+ wrap = OpSum()
+ wrap += 1.0, op_for(last_digits[1]), 1
+ for i in 2:L
+ wrap *= (op_for(last_digits[i]), i)
+ end
+ shifts += wrap
+ end
+ return MPO(shifts, sites)
+end
+
+function physical_projector(space::MetallicMeanPositionSpace, H::TBHamiltonian)
+ length(H.sites) == H.L ||
+ error("MetallicMeanPositionSpace currently supports position-only Hamiltonians")
+ return copy(space.projector)
+end
+
+function physical_site_state(space::MetallicMeanPositionSpace, H::TBHamiltonian, x::Integer)
+ 1 <= x <= H.N || throw(BoundsError(1:H.N, x))
+ length(H.sites) == H.L ||
+ error("MetallicMeanPositionSpace currently supports position-only Hamiltonians")
+ return _product_state_mps(H.sites, metallic_mean_digits(space.m, x - 1, H.L) .+ 1)
+end
+
+function site_axis(::MetallicMeanPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical || throw(ArgumentError(
+ "ordering=:$ordering is not available for MetallicMeanPositionSpace; " *
+ "only :physical is defined (conumbering is currently Fibonacci-only)"))
+ return collect(0:(H.N - 1))
+end
+
+function site_permutation(::MetallicMeanPositionSpace, H::TBHamiltonian;
+ ordering::Symbol=:physical, kwargs...)
+ ordering === :physical || throw(ArgumentError(
+ "ordering=:$ordering is not available for MetallicMeanPositionSpace; " *
+ "only :physical is defined (conumbering is currently Fibonacci-only)"))
+ return collect(1:H.N)
+end
+
+"""
+ metallic_mean_bond_symbol(m, L, bond) -> Symbol
+
+Return `:A` or `:B` for the 1-indexed bond beginning at `bond` in the canonical
+`q_L`-bond periodic approximant of the metallic-mean chain. Bond `N` joins site
+`N` to site `1` when periodic boundaries are used.
+"""
+function metallic_mean_bond_symbol(m::Integer, L::Integer, bond::Integer)
+ N = metallic_mean_site_count(m, L)
+ 1 <= bond <= N || throw(BoundsError(1:N, bond))
+ return metallic_mean_digits(m, bond - 1, L)[end] == m ? :B : :A
+end
+
+"""
+ metallic_mean_hamiltonian(m, L; A, B, model=:hopping, t=1.0, onsite=0.0,
+ boundary=:periodic, scale=nothing, padding=1.05,
+ cutoff=1e-12, maxdim=200) -> TBHamiltonian
+
+Construct the metallic-mean chain with parameter `m` (`A -> A^m B`, `B -> A`)
+in its projected numeration position space on `L` Qudit sites of dimension
+`m + 1`. The chain has `H.N = q_L` physical sites.
+
+- `model=:onsite`: `A` and `B` are onsite energies and `t` is uniform hopping.
+- `model=:hopping`: `A` and `B` are bond amplitudes and `onsite` is uniform.
+
+The Hamiltonian is assembled as `P (V + T K + h.c.) P`, where `P` is the
+validity projector, `T`/`V` the diagonal word MPO, and `K` the decrement
+[`metallic_mean_decrement_mpo`](@ref). `m = 1` reproduces the Fibonacci chain of
+[`fibonacci_hamiltonian`](@ref) on dimension-2 Qudit sites. The default periodic
+boundary closes the physical `q_L`-site approximant.
+"""
+function metallic_mean_hamiltonian(
+ m::Integer, L::Integer; A, B,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+ scale=nothing,
+ padding::Real=1.05,
+ cutoff::Real=1e-12,
+ maxdim::Integer=200,
+)
+ m >= 1 || throw(ArgumentError("metallic-mean parameter m must be at least 1"))
+ L >= 2 || throw(ArgumentError("metallic-mean chains require L >= 2"))
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ padding > 1 || throw(ArgumentError("padding must be greater than 1"))
+ if model === :onsite
+ isreal(A) && isreal(B) ||
+ throw(ArgumentError("onsite metallic-mean values A and B must be real"))
+ else
+ isreal(onsite) ||
+ throw(ArgumentError("the uniform onsite energy must be real"))
+ end
+
+ sites = siteinds("Qudit", L; dim=m + 1)
+ word_mps = _metallic_mean_automaton_mps(m, sites; A, B)
+ valid_mps = _metallic_mean_automaton_mps(m, sites; A=1.0, B=1.0)
+ word = mps_to_diagonal_mpo(word_mps, sites)
+ P = mps_to_diagonal_mpo(valid_mps, sites)
+ K = metallic_mean_decrement_mpo(m, sites; boundary)
+
+ V, TK = if model === :onsite
+ word, t * K
+ else
+ onsite * P, apply(word, K; cutoff=cutoff, maxdim=maxdim)
+ end
+ hopping = +(TK, shift_adjoint_mpo(TK); cutoff=cutoff, maxdim=maxdim)
+ Hraw = +(V, hopping; cutoff=cutoff, maxdim=maxdim)
+ mpo = apply(P, apply(Hraw, P; cutoff=cutoff, maxdim=maxdim);
+ cutoff=cutoff, maxdim=maxdim)
+ ITensorMPS.truncate!(mpo; cutoff=cutoff, maxdim=maxdim)
+
+ center, halfwidth = if model === :onsite
+ lo, hi = extrema((Float64(real(A)), Float64(real(B))))
+ ((lo + hi) / 2, (hi - lo) / 2 + 2abs(t))
+ else
+ (Float64(real(onsite)), 2max(abs(A), abs(B)))
+ end
+ scale_value = isnothing(scale) ? padding * Float64(halfwidth) : Float64(scale)
+ scale_value > 0 || throw(ArgumentError("KPM scale must be positive"))
+
+ N = metallic_mean_site_count(m, L)
+ H = TBHamiltonian(L, N, sites, mpo, _chain_geometry(),
+ scale_value, Float64(center),
+ nothing, nothing, nothing, nothing, 0, nothing)
+ H.position_space = MetallicMeanPositionSpace(Int(m), P)
+ return H
+end
+
+function _build_metallic_mean(params, L::Integer;
+ m=nothing, scale=nothing, tol=1e-12, maxdim=200,
+ kwargs...)
+ m === nothing && throw(ArgumentError(
+ "get_Hamiltonian(\"metallic_mean\", …) requires the keyword m " *
+ "(m=1 Fibonacci, m=2 silver mean, m=3 bronze mean, …)"))
+ p = if params isa NamedTuple
+ Dict{Symbol,Any}(pairs(params))
+ elseif params isa AbstractDict
+ Dict{Symbol,Any}(Symbol(k) => v for (k, v) in pairs(params))
+ else
+ throw(ArgumentError("metallic-mean parameters must be a NamedTuple or dictionary containing A and B"))
+ end
+ haskey(p, :A) && haskey(p, :B) ||
+ throw(ArgumentError("metallic-mean parameters must contain A and B"))
+ allowed = Set((:A, :B, :t, :onsite))
+ unknown = setdiff(Set(keys(p)), allowed)
+ isempty(unknown) || throw(ArgumentError("unknown metallic-mean parameters: $(collect(unknown))"))
+ return metallic_mean_hamiltonian(
+ m, L; A=p[:A], B=p[:B],
+ t=get(p, :t, 1.0), onsite=get(p, :onsite, 0.0),
+ scale=scale, cutoff=tol, maxdim=maxdim, kwargs...,
+ )
+end
+
+# Dense small-system oracle used only by the test suite.
+function _dense_metallic_mean_hamiltonian(
+ m::Integer, L::Integer; A, B,
+ model::Symbol=:hopping,
+ t::Number=1.0,
+ onsite::Number=0.0,
+ boundary::Symbol=:periodic,
+)
+ model in (:onsite, :hopping) ||
+ throw(ArgumentError("model must be :onsite or :hopping"))
+ boundary in (:open, :periodic) ||
+ throw(ArgumentError("boundary must be :open or :periodic"))
+ N = metallic_mean_site_count(m, L)
+ word = [metallic_mean_digits(m, n, L)[end] == m ? B : A for n in 0:(N - 1)]
+ diagonal = model === :onsite ? word : fill(onsite, N)
+ bonds = model === :onsite ? fill(t, N - 1) : word[1:(N - 1)]
+ H = zeros(ComplexF64, N, N)
+ H[diagind(H)] .= diagonal
+ for n in 1:(N - 1)
+ H[n, n + 1] = bonds[n]
+ H[n + 1, n] = conj(bonds[n])
+ end
+ if boundary === :periodic
+ wrap = model === :onsite ? t : word[end]
+ H[N, 1] = wrap
+ H[1, N] = conj(wrap)
+ end
+ return H
+end
diff --git a/src/solvers/KPM_tk.jl b/src/solvers/KPM_tk.jl
index bb577ca..a12bbcb 100644
--- a/src/solvers/KPM_tk.jl
+++ b/src/solvers/KPM_tk.jl
@@ -78,6 +78,7 @@ Returns `(Tn_list, scale, center)`. To convert a physical energy ω:
function KPM_Tn(H_mpo::MPO, N::Int, sites;
scale::Union{Real, Nothing} = nothing,
center::Real = 0.0,
+ identity_mpo::Union{MPO,Nothing} = nothing,
maxdim::Int = 40,
dmrg_nsweeps::Int = 5,
dmrg_maxdim = [10, 20, 40],
@@ -94,7 +95,7 @@ function KPM_Tn(H_mpo::MPO, N::Int, sites;
end
# ── Scaled Hamiltonian: (H − center·I) / scale ────────────────────────
- I_mpo = MPO(sites, "Id")
+ I_mpo = isnothing(identity_mpo) ? MPO(sites, "Id") : copy(identity_mpo)
Ham_n = (1 / scale) * +(H_mpo, (-center) * I_mpo; cutoff = cutoff)
# ── Chebyshev recursion T_0 = I, T_1 = H_scaled, T_k = 2H·T_{k-1} − T_{k-2}
@@ -205,6 +206,7 @@ function KPM_Tn(H::TBHamiltonian, Ncheb::Int;
Tn, _, _ = KPM_Tn(H.mpo, Ncheb, H.sites;
scale = H.scale,
center = H.center,
+ identity_mpo = physical_projector(H),
maxdim = maxdim,
cutoff = cutoff,
verbose = verbose)
@@ -214,6 +216,7 @@ function KPM_Tn(H::TBHamiltonian, Ncheb::Int;
Tn, _, _ = KPM_Tn_mps(H.mpo, Ncheb, psi0, H.sites;
scale = H.scale,
center = H.center,
+ identity_mpo = physical_projector(H),
maxdim = maxdim,
cutoff = cutoff,
verbose = verbose)
@@ -251,6 +254,7 @@ Returns `(Tn_mps_list, scale, center)` where `Tn_mps_list[n+1]` = |φ_n⟩.
function KPM_Tn_mps(H_mpo::MPO, N::Int, psi0::MPS, sites;
scale::Union{Real, Nothing} = nothing,
center::Real = 0.0,
+ identity_mpo::Union{MPO,Nothing} = nothing,
maxdim::Int = 40,
dmrg_nsweeps::Int = 5,
dmrg_maxdim = [10, 20, 40],
@@ -267,7 +271,7 @@ function KPM_Tn_mps(H_mpo::MPO, N::Int, psi0::MPS, sites;
end
# ── Scaled Hamiltonian: (H − center·I) / scale ────────────────────────
- I_mpo = MPO(sites, "Id")
+ I_mpo = isnothing(identity_mpo) ? MPO(sites, "Id") : copy(identity_mpo)
Ham_n = (1 / scale) * +(H_mpo, (-center) * I_mpo; cutoff = cutoff)
# ── Chebyshev recursion T_0 = |ψ₀⟩, |T_1⟩ = H_scaled|ψ₀⟩, |T_k⟩ = 2H_scaled|ψ_{k-1}⟩ − |ψ_{k-2}⟩
@@ -305,6 +309,7 @@ function KPM_Tn_mps(H::TBHamiltonian, N::Int, psi0::MPS;
Tn_mps, _, _ = KPM_Tn_mps(H.mpo, N, psi0, H.sites;
scale = H.scale,
center = H.center,
+ identity_mpo = physical_projector(H),
maxdim = maxdim,
cutoff = cutoff,
verbose = verbose)
@@ -574,7 +579,7 @@ function get_ldos_online(H::TBHamiltonian, Ncheb::Int, X::Int, ω_phys_vals;
nambu_proj, spin_proj, layer_proj, sublat_proj =
_autoenable_proj(H, nambu_proj, spin_proj, layer_proj, sublat_proj)
- I_mpo = MPO(H.sites, "Id")
+ I_mpo = physical_projector(H)
Ham_n = (1 / H.scale) * +(H.mpo, (-H.center) * I_mpo; cutoff=cutoff)
ω_vals = (collect(ω_phys_vals) .- H.center) ./ H.scale
@@ -592,7 +597,7 @@ function get_ldos_online(H::TBHamiltonian, Ncheb::Int, X::Int, ω_phys_vals;
for σ_n in nambu_range, σ_s in spin_range, σ_l in layer_range, σ_sl in sl_range
psi0 = any_aux_proj ?
_ldos_make_psi0(H, X, σ_n, σ_s, σ_l, σ_sl) :
- (L_tot == H.L ? binary_to_MPS(X - 1, H.L, H.sites) :
+ (L_tot == H.L ? physical_site_state(H, X) :
mpsexciton(X, H.sites))
_run_kpm_mps!(Ham_n, psi0, Ncheb, W, valid, accum;
cutoff=cutoff, maxdim=maxdim,
@@ -715,6 +720,12 @@ With no sublattice DOF the shape is always `(Nω × ng)`, `ng = num_x`.
**Other auxiliary DOF projections** (same interface as `get_bands`):
`nambu_proj`/`proj_nambu`, `spin_proj`/`proj_s`, `layer_proj`/`proj_layer`.
+For a Fibonacci position space, `ordering=:conumber` requires full-resolution
+point sampling. `conumber_alignment=:atomic` (default) places the `AA` sites in
+one central block; `:raw` exposes the unshifted modular residues. A recursive
+atomic zoom should slice the interval returned by `fibonacci_rg_partition`
+rather than re-conumbering its sites with a reduced `L`.
+
Examples
--------
```julia
@@ -755,6 +766,11 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
maxdim::Int = 100,
cutoff::Real = 1e-8,
verbose::Bool = false,
+ ordering::Symbol = :physical,
+ conumber_orientation::Symbol = :standard,
+ conumber_centered::Bool = true,
+ conumber_origin::Integer = 0,
+ conumber_alignment::Symbol = :atomic,
nambu_proj::Bool = false,
proj_nambu = nothing,
spin_proj::Bool = false,
@@ -764,6 +780,29 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
sublat_proj::Bool = false, # kept for backward compat; auto-on when H.sublattice_s is set
proj_sl = nothing)
+ ordering in (:physical, :conumber) ||
+ throw(ArgumentError("ordering must be :physical or :conumber"))
+ x_groups_effective = x_groups
+ if ordering === :conumber
+ full_resolution = num_x == H.N && num_y === nothing && num_avg == 1 &&
+ x_start == 1 && x_end == H.N && x_groups === nothing && !grid &&
+ xwin === nothing && ywin === nothing && box_half == 0 && reduce === :point &&
+ H.spin_s === nothing && H.nambu_s === nothing && H.layer_s === nothing &&
+ H.sublattice_s === nothing
+ full_resolution || throw(ArgumentError(
+ "ordering=:conumber currently requires full-resolution 1D point sampling " *
+ "without averaging, grids, blocks, custom groups, or auxiliary degrees of freedom"
+ ))
+ permutation = site_permutation(
+ H; ordering=:conumber,
+ orientation=conumber_orientation,
+ centered=conumber_centered,
+ origin=conumber_origin,
+ alignment=conumber_alignment,
+ )
+ x_groups_effective = [[x] for x in permutation]
+ end
+
# ── Geometry-aware sampling plan (unit-cell groups + sublattice decision) ──
if box_half > 0 || grid || xwin !== nothing || ywin !== nothing || reduce === :block
isnothing(H.geometry) &&
@@ -784,7 +823,7 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
num_x = num_x, num_y = num_y, num_avg = num_avg,
x_start = x_start, x_end = x_end,
xwin = xwin, ywin = ywin,
- x_groups = x_groups, box_half = box_half,
+ x_groups = x_groups_effective, box_half = box_half,
sublattice = sublattice)
groups = plan.groups
is_block = plan.reduce === :block
@@ -834,7 +873,7 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
(isnothing(proj_sl) ? (1:n_sub) : (proj_sl:proj_sl)) :
(1:1)
- I_mpo = MPO(H.sites, "Id")
+ I_mpo = physical_projector(H)
Ham_n = (1 / H.scale) * +(H.mpo, (-H.center) * I_mpo; cutoff=cutoff)
ω_vals = (collect(ω_phys_vals) .- H.center) ./ H.scale
@@ -861,7 +900,7 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
σ_sl in sl_fill
psi0 = any_aux_proj ?
_ldos_make_psi0(H, x, σ_n, σ_s, σ_l, σ_sl) :
- (L_tot == H.L ? binary_to_MPS(x - 1, H.L, H.sites) :
+ (L_tot == H.L ? physical_site_state(H, x) :
mpsexciton(x, H.sites))
accum_loc = zeros(Float64, Nω)
_run_kpm_mps!(Ham_n, psi0, Ncheb, W, valid, accum_loc;
@@ -901,7 +940,7 @@ function get_ldos_spatial(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
pos_sites = filter(s -> s ∉ aux_to_drop, H.sites)
psi_dict = if isempty(aux_to_drop)
- Dict(x => (L_tot == H.L ? binary_to_MPS(x - 1, H.L, H.sites) :
+ Dict(x => (L_tot == H.L ? physical_site_state(H, x) :
mpsexciton(x, H.sites)) for x in all_xs)
else
@assert length(pos_sites) == H.L "get_ldos_spatial: $(length(pos_sites)) position sites after dropping aux but expected H.L=$(H.L)."
@@ -1102,10 +1141,11 @@ function get_dos_stochastic(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
dos_weighting in (:trace, :sample) ||
error("get_dos_stochastic: dos_weighting must be :trace or :sample.")
- I_mpo = MPO(H.sites, "Id")
+ I_mpo = physical_projector(H)
Ham_n = (1 / H.scale) * +(H.mpo, (-H.center) * I_mpo; cutoff=cutoff)
- D = prod(ITensors.dim(s) for s in H.sites)
+ projected_position_space = !_is_binary_position_space(H)
+ D = projected_position_space ? H.N : prod(ITensors.dim(s) for s in H.sites)
N_phys = H.N
is_exc = length(H.sites) == 2 * H.L
@@ -1158,9 +1198,11 @@ function get_dos_stochastic(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
end
# ── Full Hilbert space samples (weight = D / N_sample per sample) ─────────
- samples = rand(rng, 0:(D - 1), N_sample)
- for (i, k) in enumerate(samples)
- psi0 = _basis_state_mps(k, H.sites)
+ samples = projected_position_space ?
+ rand(rng, 1:H.N, N_sample) : rand(rng, 0:(D - 1), N_sample)
+ for (i, sample) in enumerate(samples)
+ psi0 = projected_position_space ?
+ physical_site_state(H, sample) : _basis_state_mps(sample, H.sites)
χ = _run_kpm_mps!(Ham_n, psi0, Ncheb, W, valid, accum_full;
weight=1.0/N_sample, cutoff=cutoff, maxdim=maxdim)
verbose && i % 15 == 0 && println("Full sample $i/$N_sample maxlinkdim=$χ")
@@ -1197,6 +1239,70 @@ function get_dos_stochastic(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
end
+"""
+ get_dos_trace(H, Ncheb, ω_phys_vals; normalize=false, kernel=:jackson,
+ lambda=4.0, eta=0.0, m_order=4,
+ maxdim=100, cutoff=1e-8, verbose=false) -> Vector{Float64}
+
+Deterministic total DOS from the exact tensor-network trace of each online
+Chebyshev MPO. Only three MPOs are retained. At every order the diagonal MPO is
+converted to an MPS and contracted with the product MPS `|1,1,...>`, giving
+`Tr[T_n(H_tilde)]` without summing LDOS curves or integrating a spectrum.
+
+For projected position spaces, `T_0` is `physical_projector(H)` and the trace is
+therefore over physical states only. `normalize=true` divides by `Tr(T_0)`;
+otherwise the spectral weight corresponds to the total traced state count.
+"""
+function get_dos_trace(H::TBHamiltonian, Ncheb::Int, ω_phys_vals;
+ normalize::Bool=false,
+ kernel::Symbol=:jackson,
+ lambda::Real=4.0,
+ eta::Real=0.0,
+ m_order::Int=4,
+ maxdim::Int=100,
+ cutoff::Real=1e-8,
+ verbose::Bool=false)
+ Ncheb >= 2 || throw(ArgumentError("Ncheb must be at least 2"))
+ _ensure_scale!(H)
+ P = physical_projector(H)
+ Ham_n = (1 / H.scale) * +(H.mpo, (-H.center) * P; cutoff=cutoff)
+
+ function trace_diagonal(Tn::MPO)
+ diagonal = extract_diagonal_to_mps(Tn)
+ ITensorMPS.truncate!(diagonal; cutoff=cutoff, maxdim=maxdim)
+ ones_state = MPS([ITensor(ones(Float64, dim(s)), s)
+ for s in siteinds(diagonal)])
+ return real(inner(ones_state, diagonal))
+ end
+
+ moments = zeros(Float64, Ncheb)
+ Tkm2 = P
+ Tkm1 = Ham_n
+ moments[1] = trace_diagonal(Tkm2)
+ moments[2] = trace_diagonal(Tkm1)
+ for k in 3:Ncheb
+ Tk = +(2 * apply(Ham_n, Tkm1; cutoff=cutoff), -Tkm2;
+ cutoff=cutoff, maxdim=maxdim)
+ ITensorMPS.truncate!(Tk; cutoff=cutoff, maxdim=maxdim)
+ moments[k] = trace_diagonal(Tk)
+ Tkm2, Tkm1 = Tkm1, Tk
+ verbose && (k % 10 == 0 || k == Ncheb) &&
+ println("get_dos_trace step $k/$Ncheb maxlinkdim=$(maxlinkdim(Tkm1))")
+ end
+
+ ω_vals = (collect(ω_phys_vals) .- H.center) ./ H.scale
+ W, denom = _dos_weight_matrix(Ncheb, ω_vals;
+ kernel, lambda, eta, m_order)
+ result = zeros(Float64, length(ω_vals))
+ for iω in eachindex(ω_vals)
+ abs(ω_vals[iω]) < 1 || continue
+ result[iω] = dot(moments, view(W, :, iω)) / denom[iω]
+ end
+ normalize && (result ./= moments[1])
+ return result
+end
+
+
"""
get_ldos_from_mun(mun_list, N, E; kernel=:jackson, lambda=4.0) -> Real
@@ -1633,6 +1739,11 @@ are generated over `x_start:x_end`, with `num_avg` subpositions per group.
`kernel=:hodc` uses the HODC reconstruction (`eta`, `m_order`); otherwise the
standard KPM kernels are available (`:jackson`, `:lorentz`, `:fejer`,
`:dirichlet`).
+
+`return_maxlinkdim=true` returns `(result, linkdims)` instead of just `result`,
+where `linkdims::Vector{Int}` is the reached MPS bond dimension per output column
+(the χ the Chebyshev recursion hit under the given `maxdim`/`cutoff`). Mirrors the
+GPU entry point; useful for cutoff/tolerance studies where χ is the observable.
"""
function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
X_list = nothing,
@@ -1649,7 +1760,8 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
maxdim::Int = 100,
cutoff::Real = 1e-8,
verbose::Bool = false,
- printinfo::Bool = false)
+ printinfo::Bool = false,
+ return_maxlinkdim::Bool = false)
_ensure_scale!(H)
length(H.sites) == 2 * H.L ||
error("get_exciton_ldos_spatial: H is not an exciton Hamiltonian (expected length(H.sites) == 2*H.L).")
@@ -1661,9 +1773,7 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
group_arg = X_groups !== nothing ? X_groups : x_groups
groups = if group_arg !== nothing
- group_arg isa AbstractVector{<:AbstractVector} ?
- [collect(Int, grp) for grp in group_arg] :
- [[Int(x)] for x in group_arg]
+ spatial_sampling_plan(H.L; x_groups=group_arg).groups
elseif X_list !== nothing
[[Int(x)] for x in X_list]
else
@@ -1674,12 +1784,9 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
window = x_end - x_start + 1
num_x <= window ||
error("get_exciton_ldos_spatial: num_x=$num_x exceeds sampling window length $window.")
- dx = div(window, num_x)
- dx_sub = max(1, div(dx, num_avg))
- [[x_start + (i - 1) * dx + k * dx_sub
- for k in 0:num_avg-1
- if x_start + (i - 1) * dx + k * dx_sub <= x_end]
- for i in 1:num_x]
+ # 1D point layout of the shared planner (core/Utils.jl): stride
+ # window ÷ num_x with num_avg sub-probes per coarse cell.
+ spatial_sampling_plan(H.L; num_x, num_avg, x_start, x_end).groups
end
isempty(groups) && error("get_exciton_ldos_spatial: no spatial groups were selected.")
@@ -1702,6 +1809,7 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
nX = length(groups)
Xs = first.(groups)
result = zeros(Float64, Nomega, nX)
+ linkdims = zeros(Int, nX) # reached MPS bond dim per output column (see return_maxlinkdim)
for (j, group) in enumerate(groups)
last_linkdim = 0
@@ -1718,12 +1826,13 @@ function get_exciton_ldos_spatial(H::TBHamiltonian, Ncheb::Int, omega_phys_vals;
valid[iomega] || continue
result[iomega, j] = accum_group[iomega] / denom[iomega]
end
+ linkdims[j] = last_linkdim
(verbose || printinfo) && (j % 5 == 0 || j == nX) &&
println(" exciton ldos $j/$nX (X=$(Xs[j]), n_avg=$(length(group))) maxlinkdim=$last_linkdim")
end
- return result
+ return return_maxlinkdim ? (result, linkdims) : result
end
function get_exciton_ldos(H::TBHamiltonian, X::Int, omega_phys::Real;
diff --git a/test/fibonacci.jl b/test/fibonacci.jl
new file mode 100644
index 0000000..2a78ed9
--- /dev/null
+++ b/test/fibonacci.jl
@@ -0,0 +1,221 @@
+using Test
+using LinearAlgebra
+using TensorBinding
+using ITensors
+using ITensorMPS
+
+function _projected_fibonacci_matrix(H)
+ basis = [TensorBinding.physical_site_state(H, x) for x in 1:H.N]
+ matrix = zeros(ComplexF64, H.N, H.N)
+ for j in eachindex(basis)
+ Hket = apply(H.mpo, basis[j]; cutoff=1e-13, maxdim=300)
+ for i in eachindex(basis)
+ matrix[i, j] = inner(basis[i], Hket)
+ end
+ end
+ return matrix
+end
+
+function _dense_kpm_moments(decomposition, Ncheb, center, scale; site=nothing)
+ scaled = clamp.((decomposition.values .- center) ./ scale, -1.0, 1.0)
+ angles = acos.(scaled)
+ weights = isnothing(site) ? ones(length(scaled)) :
+ abs2.(decomposition.vectors[site, :])
+ return [sum(weights .* cos.(n .* angles)) for n in 0:(Ncheb - 1)]
+end
+
+@testset "Fibonacci projected position space" begin
+ H = TensorBinding.fibonacci_hamiltonian(4; A=1.0, B=2.0)
+ H_dispatch = TensorBinding.get_Hamiltonian(
+ "fibonacci", (A=1.0, B=2.0); L=4,
+ )
+
+ @test H.N == 8
+ @test TensorBinding.ambient_dimension(H) == 16
+ @test TensorBinding.ambient_dimension(H) isa BigInt
+ @test H.position_space isa TensorBinding.FibonacciPositionSpace
+ @test H_dispatch.N == H.N
+ @test H.scale > 0
+ @test H.center == 0
+
+ P = TensorBinding.physical_projector(H)
+ @test real(tr(P)) ≈ H.N atol=1e-12
+ @test norm(apply(P, P; cutoff=1e-13) - P) / norm(P) < 1e-12
+
+ @testset "TN and dense construction" begin
+ for model in (:onsite, :hopping), boundary in (:open, :periodic)
+ parameters = (; A=1.2, B=0.7, t=0.9, onsite=0.2,
+ model, boundary)
+ Htn = TensorBinding.fibonacci_hamiltonian(4; parameters...)
+ Hdense = TensorBinding._dense_fibonacci_hamiltonian(4; parameters...)
+ matrix = _projected_fibonacci_matrix(Htn)
+ @test maximum(abs.(matrix .- Hdense)) < 1e-11
+ @test norm(matrix - matrix') < 1e-11
+
+ expected_wrap = boundary === :open ? 0.0 :
+ (model === :onsite ? parameters.t :
+ (iszero(TensorBinding.fibonacci_zeckendorf_digits(Htn.N - 1, 4)[end]) ?
+ parameters.A : parameters.B))
+ @test matrix[end, 1] ≈ expected_wrap atol=1e-11
+ end
+ end
+
+ @test_throws ArgumentError TensorBinding.fibonacci_hamiltonian(
+ 4; A=1 + 1im, B=2.0, model=:onsite,
+ )
+ Hcomplex = TensorBinding.fibonacci_hamiltonian(
+ 4; A=1 + 0.2im, B=2 - 0.1im, model=:hopping,
+ )
+ complex_matrix = _projected_fibonacci_matrix(Hcomplex)
+ @test norm(complex_matrix - complex_matrix') < 1e-11
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "fibonacci", (A=1.0,); L=4,
+ )
+
+ @testset "KPM projector and spectra" begin
+ Ncheb = 8
+ Tn, _, _ = TensorBinding.KPM_Tn(
+ H, Ncheb; maxdim=100, cutoff=1e-12,
+ )
+ @test real(tr(Tn[1])) ≈ H.N atol=1e-10
+ @test norm(Tn[1] - P) < 1e-12
+
+ dense = TensorBinding._dense_fibonacci_hamiltonian(
+ 4; A=1.0, B=2.0,
+ )
+ decomposition = eigen(Hermitian(dense))
+ dense_moments = _dense_kpm_moments(
+ decomposition, Ncheb, H.center, H.scale,
+ )
+ tn_moments = real.(tr.(Tn[1:Ncheb]))
+ @test maximum(abs.(tn_moments .- dense_moments)) < 1e-8
+
+ energies = collect(range(-3.5, 3.5; length=7))
+ dos_tn = TensorBinding.get_dos_trace(
+ H, Ncheb, energies; maxdim=100, cutoff=1e-12,
+ )
+ dos_dense = [
+ TensorBinding.get_ldos_from_mun(
+ dense_moments, Ncheb, (energy - H.center) / H.scale,
+ ) for energy in energies
+ ]
+ @test maximum(abs.(dos_tn .- dos_dense)) < 1e-8
+
+ ldos_tn = TensorBinding.get_ldos_spatial(
+ H, Ncheb, energies;
+ mode=:mps, ordering=:conumber,
+ maxdim=100, cutoff=1e-12,
+ )
+ permutation = TensorBinding.site_permutation(H; ordering=:conumber)
+ ldos_dense = zeros(length(energies), H.N)
+ for site in 1:H.N
+ moments = _dense_kpm_moments(
+ decomposition, Ncheb, H.center, H.scale; site,
+ )
+ for (i, energy) in pairs(energies)
+ ldos_dense[i, site] = TensorBinding.get_ldos_from_mun(
+ moments, Ncheb, (energy - H.center) / H.scale,
+ )
+ end
+ end
+ @test maximum(abs.(ldos_tn .- ldos_dense[:, permutation])) < 1e-8
+
+ dos_stochastic = TensorBinding.get_dos_stochastic(
+ H, 4, [0.0]; N_sample=100, seed=7,
+ maxdim=60, cutoff=1e-12,
+ )
+ dos_exact = TensorBinding.get_dos_trace(
+ H, 4, [0.0]; maxdim=60, cutoff=1e-12,
+ )
+ @test all(isfinite, dos_stochastic)
+ @test abs(dos_stochastic[1] - dos_exact[1]) < 0.35 * max(abs(dos_exact[1]), 1.0)
+ end
+
+ @testset "Conumbering and guards" begin
+ axis = TensorBinding.site_axis(H; ordering=:conumber)
+ permutation = TensorBinding.site_permutation(H; ordering=:conumber)
+ reflected = TensorBinding.site_permutation(
+ H; ordering=:conumber, orientation=:reversed,
+ )
+ @test axis == collect(-4:3)
+ @test sort(permutation) == collect(1:H.N)
+ @test sort(reflected) == collect(1:H.N)
+ @test permutation != reflected
+
+ classes = TensorBinding.fibonacci_site_environment.(4, permutation)
+ @test all(!=(:atomic), classes[1:3])
+ @test classes[4:5] == fill(:atomic, 2)
+ @test all(!=(:atomic), classes[6:8])
+ for Ltest in 4:9
+ Ntest = TensorBinding.fibonacci_site_count(Ltest)
+ sites_by_c = sortperm([
+ TensorBinding.fibonacci_conumber(
+ Ltest, site; centered=false,
+ ) for site in 1:Ntest
+ ])
+ environments = TensorBinding.fibonacci_site_environment.(
+ Ltest, sites_by_c,
+ )
+ molecular_count = Int(TensorBinding.fibonacci_number(Ltest))
+ atomic_count = Int(TensorBinding.fibonacci_number(Ltest - 1))
+ @test all(!=(:atomic), environments[1:molecular_count])
+ @test environments[(molecular_count + 1):(molecular_count + atomic_count)] ==
+ fill(:atomic, atomic_count)
+ @test all(!=(:atomic), environments[(molecular_count + atomic_count + 1):end])
+ @test all(site -> TensorBinding.fibonacci_site_from_conumber(
+ Ltest,
+ TensorBinding.fibonacci_conumber(Ltest, site),
+ ) == site, 1:Ntest)
+ end
+
+ Llarge = 43
+ large_partition = TensorBinding.fibonacci_rg_partition(Llarge)
+ @test TensorBinding.fibonacci_site_count(Llarge) == 1_134_903_170
+ @test large_partition.molecular_count == 433_494_437
+ @test large_partition.atomic_count == 267_914_296
+
+ # The largest Fibonacci register whose physical site count fits Int64
+ # still needs overflow-safe modular shifts and a BigInt ambient size.
+ Lmax = 90
+ Nmax = TensorBinding.fibonacci_site_count(Lmax)
+ site_at_high_raw_rank = TensorBinding.fibonacci_site_from_conumber(
+ Lmax, Nmax - 1; alignment=:raw, centered=false,
+ )
+ atomic_conumber = TensorBinding.fibonacci_conumber(
+ Lmax, site_at_high_raw_rank; alignment=:atomic, centered=false,
+ )
+ @test 0 <= atomic_conumber < Nmax
+ @test TensorBinding.fibonacci_site_from_conumber(
+ Lmax, atomic_conumber; alignment=:atomic, centered=false,
+ ) == site_at_high_raw_rank
+
+ Hwide = deepcopy(H)
+ Hwide.L = Lmax
+ @test TensorBinding.ambient_dimension(Hwide) == big(2)^Lmax
+
+ deep = TensorBinding.fibonacci_rg_partition(Llarge; depth=13)
+ @test deep.effective_L == 4
+ @test deep.window_count == 8
+ @test deep.molecular_count == 3
+ @test deep.atomic_count == 2
+ deep_sites = [TensorBinding.fibonacci_site_from_conumber(
+ Llarge, c; centered=false,
+ ) for c in deep.window_ranks]
+ @test all(==(:atomic), TensorBinding.fibonacci_site_environment.(
+ Llarge, deep_sites,
+ ))
+ @test all(>=(13), TensorBinding.fibonacci_atomic_depth.(
+ Llarge, deep_sites,
+ ))
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial(
+ H, 4, [0.0]; ordering=:conumber, num_x=4,
+ )
+ @test_throws ArgumentError TensorBinding.add_onsite!(H, 0.1)
+ @test_throws ArgumentError TensorBinding.get_bands(H, 4, 1, [0.0])
+ end
+
+ Hbinary = TensorBinding.get_Hamiltonian("chain_1d", 1.0; L=3)
+ @test Hbinary.position_space isa TensorBinding.BinaryPositionSpace
+ @test real(tr(TensorBinding.physical_projector(Hbinary))) ≈ Hbinary.N atol=1e-12
+ @test TensorBinding.site_axis(Hbinary) == collect(0:7)
+end
diff --git a/test/fibonacci_sampling.jl b/test/fibonacci_sampling.jl
new file mode 100644
index 0000000..27d0e3a
--- /dev/null
+++ b/test/fibonacci_sampling.jl
@@ -0,0 +1,206 @@
+using Test
+
+@testset "Fibonacci LDOS sampling planner" begin
+ @testset "exact balanced intervals and samples" begin
+ plan = TensorBinding.fibonacci_ldos_sampling_plan(
+ 4; num_x=3, num_avg=2,
+ )
+
+ @test plan.L == 4
+ @test plan.N == 8
+ @test plan.depth == 0
+ @test plan.effective_L == 4
+ @test plan.num_x == 3
+ @test plan.num_x_requested == 3
+ @test plan.num_avg == 2
+ @test plan.intervals == [(-4):(-2), (-1):1, 2:3]
+ @test plan.interval_rank_first == [0, 3, 6]
+ @test plan.interval_rank_last == [2, 5, 7]
+ @test plan.sample_ranks == [[0, 2], [3, 5], [6, 7]]
+ @test plan.sample_conumbers == [[-4, -2], [-1, 1], [2, 3]]
+ @test plan.conumber_axis == [-3, 0, 2]
+ @test plan.groups == [[8, 6], [1, 7], [2, 5]]
+ @test plan.centers == [3, 4, 2]
+ @test plan.group_sizes == [2, 2, 2]
+ @test plan.group_offsets == [1, 3, 5, 7]
+ @test plan.group_offsets_zero == [0, 2, 4, 6]
+ @test plan.group_offsets_base == 1
+ @test plan.metadata.group_offsets_base == 1
+ @test plan.sample_sites_flat == [8, 6, 1, 7, 2, 5]
+ @test plan.sample_conumbers_flat == [-4, -2, -1, 1, 2, 3]
+ @test plan.column_indices == [1, 2, 3]
+ @test plan.total_samples == 6
+
+ for column in plan.column_indices
+ stored = plan.group_offsets[column]:(plan.group_offsets[column + 1] - 1)
+ @test plan.sample_sites_flat[stored] == plan.groups[column]
+ @test plan.sample_conumbers_flat[stored] == plan.sample_conumbers[column]
+
+ # Python/HDF5 consumers use the zero-based half-open slice
+ # flat[offsets[i]:offsets[i+1]]. Translate it by one for Julia.
+ stored_zero = (plan.group_offsets_zero[column] + 1):plan.group_offsets_zero[column + 1]
+ @test plan.sample_sites_flat[stored_zero] == plan.groups[column]
+ @test plan.sample_conumbers_flat[stored_zero] ==
+ plan.sample_conumbers[column]
+ end
+ end
+
+ @testset "inherited zoom coordinates" begin
+ # L=10 -> L=4 after two atomic deflations. The original uncentered
+ # ranks are 68:75, not a freshly assigned 0:7 reduced-chain axis.
+ zoom = TensorBinding.fibonacci_ldos_sampling_plan(
+ 10; depth=2, num_x=3, num_avg=3, centered=false,
+ )
+ @test zoom.effective_L == 4
+ @test zoom.window_count == 8
+ @test zoom.window_first == 68
+ @test zoom.window_last == 75
+ @test zoom.intervals == [68:70, 71:73, 74:75]
+ @test zoom.sample_conumbers == [[68, 69, 70], [71, 72, 73], [74, 75]]
+ @test zoom.conumber_axis == [69, 72, 74]
+ @test all(TensorBinding.fibonacci_conumber(
+ 10, zoom.groups[column][sample]; centered=false,
+ ) == zoom.sample_conumbers[column][sample]
+ for column in eachindex(zoom.groups)
+ for sample in eachindex(zoom.groups[column]))
+
+ # More requested columns/samples than sites produces singleton groups,
+ # never duplicate samples or empty intervals.
+ tiny = TensorBinding.fibonacci_ldos_sampling_plan(
+ 4; depth=0, num_x=100, num_avg=100,
+ )
+ @test tiny.num_x == 8
+ @test tiny.group_sizes == ones(Int, 8)
+ @test tiny.intervals == [conumber:conumber for conumber in (-4):3]
+ @test tiny.groups == [[TensorBinding.fibonacci_site_from_conumber(4, c)]
+ for c in (-4):3]
+ end
+
+ @testset "conumber conventions" begin
+ standard = TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=1, num_x=7, num_avg=4,
+ orientation=:standard, alignment=:atomic,
+ centered=true, origin=0,
+ )
+ reversed = TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=1, num_x=7, num_avg=4,
+ orientation=:reversed, alignment=:atomic,
+ centered=true, origin=0,
+ )
+ @test standard.intervals == reversed.intervals
+ @test standard.sample_conumbers == reversed.sample_conumbers
+ @test standard.groups != reversed.groups
+ @test all(==(:atomic), TensorBinding.fibonacci_site_environment.(
+ standard.L, standard.sample_sites_flat,
+ ))
+ @test all(==(:atomic), TensorBinding.fibonacci_site_environment.(
+ reversed.L, reversed.sample_sites_flat,
+ ))
+
+ for plan in (standard, reversed)
+ @test all(TensorBinding.fibonacci_conumber(
+ plan.L, plan.groups[column][sample];
+ orientation=plan.orientation,
+ alignment=plan.alignment,
+ centered=plan.centered,
+ origin=plan.origin,
+ ) == plan.sample_conumbers[column][sample]
+ for column in eachindex(plan.groups)
+ for sample in eachindex(plan.groups[column]))
+ end
+
+ # A noncanonical phase is meaningful for a complete depth-zero view,
+ # but cannot be described as the nested atomic RG window.
+ canonical_full = TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=0, num_x=7, num_avg=4,
+ )
+ raw = TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=0, num_x=7, num_avg=4, alignment=:raw,
+ )
+ shifted = TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=0, num_x=7, num_avg=4, origin=3,
+ )
+ @test canonical_full.intervals == raw.intervals == shifted.intervals
+ @test canonical_full.groups != raw.groups
+ @test canonical_full.groups != shifted.groups
+ for plan in (raw, shifted)
+ @test all(TensorBinding.fibonacci_conumber(
+ plan.L, plan.groups[column][sample];
+ orientation=plan.orientation,
+ alignment=plan.alignment,
+ centered=plan.centered,
+ origin=plan.origin,
+ ) == plan.sample_conumbers[column][sample]
+ for column in eachindex(plan.groups)
+ for sample in eachindex(plan.groups[column]))
+ end
+
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=1, alignment=:raw,
+ )
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(
+ 8; depth=1, origin=1,
+ )
+ end
+
+ @testset "L=43 remains output-sized" begin
+ # Warm the exact specialization before measuring allocations.
+ plan = TensorBinding.fibonacci_ldos_sampling_plan(
+ 43; depth=3, num_x=100, num_avg=5,
+ )
+ bytes = @allocated TensorBinding.fibonacci_ldos_sampling_plan(
+ 43; depth=3, num_x=100, num_avg=5,
+ )
+
+ @test plan.N == 1_134_903_170
+ @test plan.effective_L == 34
+ @test plan.num_x == 100
+ @test length(plan.groups) == 100
+ @test plan.total_samples == 500
+ @test length(plan.sample_sites_flat) == 500
+ @test all(length(group) == 5 for group in plan.groups)
+ @test all(1 <= site <= plan.N for site in plan.sample_sites_flat)
+ @test first(first(plan.intervals)) == plan.window_first
+ @test last(last(plan.intervals)) == plan.window_last
+ @test sum(length, plan.intervals) == plan.window_count
+ @test all(last(plan.intervals[i]) + 1 == first(plan.intervals[i + 1])
+ for i in 1:(plan.num_x - 1))
+ @test bytes < 20_000_000
+
+ deepest_requested_view = TensorBinding.fibonacci_ldos_sampling_plan(
+ 43; depth=12, num_x=100, num_avg=5,
+ )
+ @test deepest_requested_view.effective_L == 7
+ @test deepest_requested_view.window_count == 34
+ @test deepest_requested_view.num_x == 34
+ @test deepest_requested_view.group_sizes == ones(Int, 34)
+ end
+
+ @testset "Hamiltonian overload and validation" begin
+ Hfib = TensorBinding.fibonacci_hamiltonian(
+ 2; A=1.0, B=2.0, boundary=:open,
+ )
+ from_H = TensorBinding.fibonacci_ldos_sampling_plan(
+ Hfib; num_x=2, num_avg=1,
+ )
+ from_L = TensorBinding.fibonacci_ldos_sampling_plan(
+ 2; num_x=2, num_avg=1,
+ )
+ @test from_H.groups == from_L.groups
+ @test from_H.metadata == from_L.metadata
+
+ Hbinary = TensorBinding.get_Hamiltonian("chain_1d", 1.0; L=2)
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(Hbinary)
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(4; num_x=0)
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(4; num_avg=0)
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(
+ 4; orientation=:sideways,
+ )
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(
+ 4; alignment=:molecular,
+ )
+ @test_throws ArgumentError TensorBinding.fibonacci_ldos_sampling_plan(
+ 4; depth=1,
+ )
+ end
+end
diff --git a/test/gpu_mps_ldos.jl b/test/gpu_mps_ldos.jl
new file mode 100644
index 0000000..40a40a5
--- /dev/null
+++ b/test/gpu_mps_ldos.jl
@@ -0,0 +1,152 @@
+using Test
+using LinearAlgebra
+using TensorBinding
+
+@testset "GPU MPS spatial LDOS interface" begin
+ H = TensorBinding.get_Hamiltonian("chain_1d", 1.0; L=3, scale=2.5)
+ energies = [-1.0, 0.0, 1.0]
+
+ @testset "moment-column reconstruction" begin
+ moments = [1.0 2.0; 0.5 -1.0; -0.25 0.75]
+ weights = [1.0 2.0 3.0; 0.5 -1.0 4.0; 2.0 0.25 -2.0]
+ denom = [2.0, 4.0, 0.0]
+ valid = [true, true, false]
+ reconstructed = TensorBinding._reconstruct_ldos_moment_columns(
+ moments, weights, denom, valid,
+ )
+ expected = transpose(weights) * moments
+ expected[1, :] ./= denom[1]
+ expected[2, :] ./= denom[2]
+ expected[3, :] .= 0.0
+ @test reconstructed == expected
+ @test size(reconstructed) == (size(weights, 2), size(moments, 2))
+ @test_throws DimensionMismatch TensorBinding._reconstruct_ldos_moment_columns(
+ moments[1:2, :], weights, denom, valid,
+ )
+ @test_throws DimensionMismatch TensorBinding._reconstruct_ldos_moment_columns(
+ moments, weights, denom[1:2], valid,
+ )
+ @test_throws DimensionMismatch TensorBinding._reconstruct_ldos_moment_columns(
+ moments, weights, denom, valid[1:2],
+ )
+ end
+
+ # These checks deliberately run before CUDA discovery: unsupported requests
+ # should fail at the public API boundary, not deep inside the GPU backend.
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 1, energies; x_groups=[1],
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=[1], reduce=:block,
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=[1], grid=true,
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=[1], ordering=:conumber,
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=[1], spin_proj=true,
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=Vector{Vector{Int}}(),
+ )
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 4, energies; x_groups=[[0]],
+ )
+ Haux = TensorBinding.get_Hamiltonian("ssh_sublattice", (t=1.0, d=0.2); L=3)
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial_mps_gpu(
+ Haux, 4, energies; x_groups=[1],
+ )
+
+ cuda_functional = false
+ cuda_error = nothing
+ if Base.find_package("CUDA") !== nothing
+ try
+ @eval using CUDA
+ cuda_functional = CUDA.functional()
+ catch err
+ cuda_error = err
+ end
+ end
+
+ if cuda_functional
+ @testset "projected Fibonacci CPU/GPU agreement" begin
+ Hf = TensorBinding.fibonacci_hamiltonian(
+ 4; A=1.0, B=2.0, model=:onsite, t=0.6,
+ boundary=:open,
+ cutoff=1e-12, maxdim=100,
+ )
+ @test Hf.center != 0.0 # exercises the projected shift H - center*P
+ groups = [[1, 2], [4], [7, 8]]
+ Ncheb = 8
+ ω = collect(range(-0.1, 3.1; length=9))
+
+ cpu = TensorBinding.get_ldos_spatial(
+ Hf, Ncheb, ω;
+ mode=:mps, x_groups=groups,
+ maxdim=100, cutoff=1e-10,
+ )
+ gpu, moments, linkdims = TensorBinding.get_ldos_spatial_mps_gpu(
+ Hf, Ncheb, ω;
+ x_groups=groups,
+ type=ComplexF32, maxdim=100, cutoff=1e-6,
+ return_maxlinkdim=true,
+ return_moments=true,
+ )
+
+ @test size(gpu) == (length(ω), length(groups))
+ @test size(moments) == (Ncheb, length(groups))
+ @test length(linkdims) == length(groups)
+ @test all(>=(1), linkdims)
+ @test gpu ≈ cpu rtol=5e-4 atol=5e-5
+
+ dense = TensorBinding._dense_fibonacci_hamiltonian(
+ 4; A=1.0, B=2.0, model=:onsite, t=0.6, boundary=:open,
+ )
+ decomposition = eigen(Hermitian(dense))
+ scaled_eigenvalues = clamp.(
+ (decomposition.values .- Hf.center) ./ Hf.scale, -1.0, 1.0,
+ )
+ eigenangles = acos.(scaled_eigenvalues)
+ dense_site_moments(site) = [
+ sum(
+ abs2.(decomposition.vectors[site, :]) .*
+ cos.(n .* eigenangles)
+ ) for n in 0:(Ncheb - 1)
+ ]
+ dense_group_moments = hcat([
+ sum(
+ dense_site_moments(x) for x in group
+ ) ./ length(group)
+ for group in groups
+ ]...)
+ @test moments ≈ dense_group_moments rtol=5e-4 atol=5e-5
+
+ ω_scaled = (ω .- Hf.center) ./ Hf.scale
+ W, denom = TensorBinding._dos_weight_matrix(Ncheb, ω_scaled)
+ reconstructed = TensorBinding._reconstruct_ldos_moment_columns(
+ moments, W, denom, abs.(ω_scaled) .< 1.0,
+ )
+ @test reconstructed ≈ gpu rtol=5e-13 atol=5e-13
+ end
+
+ @testset "ordinary binary position space" begin
+ groups = [[1], [3, 4]]
+ cpu = TensorBinding.get_ldos_spatial(
+ H, 6, energies;
+ mode=:mps, x_groups=groups,
+ maxdim=40, cutoff=1e-10,
+ )
+ gpu = TensorBinding.get_ldos_spatial_mps_gpu(
+ H, 6, energies;
+ x_groups=groups,
+ type=ComplexF32, maxdim=40, cutoff=1e-6,
+ )
+ @test gpu ≈ cpu rtol=5e-4 atol=5e-5
+ end
+ else
+ @info "Skipping CUDA-functional GPU MPS LDOS comparisons" exception=cuda_error
+ @test true
+ end
+end
diff --git a/test/kbonacci.jl b/test/kbonacci.jl
new file mode 100644
index 0000000..e623e69
--- /dev/null
+++ b/test/kbonacci.jl
@@ -0,0 +1,250 @@
+using Test
+using LinearAlgebra
+using TensorBinding
+using ITensors
+using ITensorMPS
+
+# Substitution fixed point a_i -> a_1 a_(i+1) (i < k), a_k -> a_1 on the letters
+# A, B, C, …: an oracle independent of the numeration system used by the package.
+function _kbonacci_word(k::Integer, n::Integer)
+ letters = ['A' + i for i in 0:(k - 1)]
+ word = "A"
+ while length(word) < n
+ word = join(c == letters[k] ? "A" :
+ "A" * letters[findfirst(==(c), letters) + 1] for c in word)
+ end
+ return word[1:n]
+end
+
+_kb_trailing_ones(d) = (r = 0; for x in reverse(d); x == 1 ? (r += 1) : break; end; r)
+
+function _projected_kb_matrix(H)
+ basis = [TensorBinding.physical_site_state(H, x) for x in 1:H.N]
+ matrix = zeros(ComplexF64, H.N, H.N)
+ for j in eachindex(basis)
+ Hket = apply(H.mpo, basis[j]; cutoff=1e-13, maxdim=300)
+ for i in eachindex(basis)
+ matrix[i, j] = inner(basis[i], Hket)
+ end
+ end
+ return matrix
+end
+
+function _kb_dense_kpm_moments(decomposition, Ncheb, center, scale; site=nothing)
+ scaled = clamp.((decomposition.values .- center) ./ scale, -1.0, 1.0)
+ angles = acos.(scaled)
+ weights = isnothing(site) ? ones(length(scaled)) :
+ abs2.(decomposition.vectors[site, :])
+ return [sum(weights .* cos.(n .* angles)) for n in 0:(Ncheb - 1)]
+end
+
+@testset "k-bonacci projected position space" begin
+ @testset "numeration system and word" begin
+ @test [Int(TensorBinding.kbonacci_number(2, n)) for n in 0:8] == [1, 2, 3, 5, 8, 13, 21, 34, 55]
+ @test [Int(TensorBinding.kbonacci_number(3, n)) for n in 0:8] == [1, 2, 4, 7, 13, 24, 44, 81, 149]
+ @test [Int(TensorBinding.kbonacci_number(4, n)) for n in 0:8] == [1, 2, 4, 8, 15, 29, 56, 108, 208]
+ @test TensorBinding.kbonacci_number(3, 60) isa BigInt
+ for k in 2:5, L in 2:6
+ N = TensorBinding.kbonacci_site_count(k, L)
+ @test N == (L < k ? 2^L : sum(TensorBinding.kbonacci_site_count(k, L - i) for i in 1:k))
+ weights = [Int(TensorBinding.kbonacci_number(k, L - p)) for p in 1:L]
+ strings = [TensorBinding.kbonacci_digits(k, n, L) for n in 0:(N - 1)]
+ @test all(d -> all(0 .<= d .<= 1), strings)
+ @test all(d -> !occursin("1"^k, join(d)), strings)
+ @test [sum(d .* weights) for d in strings] == collect(0:(N - 1))
+ @test issorted(strings)
+ brute = sort([[(s >> (L - p)) & 1 for p in 1:L] for s in 0:(2^L - 1)
+ if !occursin("1"^k, string(s; base=2, pad=L))])
+ @test brute == strings
+ @test strings[end] == [(p % k == 0) ? 0 : 1 for p in 1:L]
+
+ letters = [TensorBinding.kbonacci_bond_symbol(k, L, b) for b in 1:N]
+ @test join(string.(letters)) == _kbonacci_word(k, N)
+ @test letters == [TensorBinding.kbonacci_letters(k)[_kb_trailing_ones(d) + 1] for d in strings]
+ end
+ @test TensorBinding.kbonacci_letters(3) == [:A, :B, :C]
+ @test_throws ArgumentError TensorBinding.kbonacci_digits(3, 13, 4)
+ @test_throws ArgumentError TensorBinding.kbonacci_number(1, 3)
+ @test_throws ArgumentError TensorBinding.kbonacci_number(3, -1)
+ @test_throws ArgumentError TensorBinding.kbonacci_letters(1)
+ @test_throws BoundsError TensorBinding.kbonacci_bond_symbol(3, 4, 14)
+ end
+
+ @testset "k = 2 reproduces the Fibonacci chain" begin
+ for L in 2:7
+ N = TensorBinding.fibonacci_site_count(L)
+ @test TensorBinding.kbonacci_site_count(2, L) == N
+ @test all(n -> TensorBinding.kbonacci_digits(2, n, L) ==
+ TensorBinding.fibonacci_zeckendorf_digits(n, L), 0:(N - 1))
+ end
+ for model in (:onsite, :hopping), boundary in (:open, :periodic)
+ dense_kb = TensorBinding._dense_kbonacci_hamiltonian(
+ 2, 4; values=(1.2, 0.7), t=0.9, onsite=0.2, model, boundary,
+ )
+ dense_fib = TensorBinding._dense_fibonacci_hamiltonian(
+ 4; A=1.2, B=0.7, t=0.9, onsite=0.2, model, boundary,
+ )
+ @test dense_kb == dense_fib
+ Hkb = TensorBinding.kbonacci_hamiltonian(
+ 2, 4; values=(1.2, 0.7), t=0.9, onsite=0.2, model, boundary,
+ )
+ Hfib = TensorBinding.fibonacci_hamiltonian(
+ 4; A=1.2, B=0.7, t=0.9, onsite=0.2, model, boundary,
+ )
+ @test Hkb.N == 8
+ @test all(s -> dim(s) == 2, Hkb.sites)
+ @test Hkb.scale ≈ Hfib.scale && Hkb.center ≈ Hfib.center
+ @test maximum(abs.(_projected_kb_matrix(Hkb) .- dense_fib)) < 1e-11
+ end
+ end
+
+ @testset "TN construction matches the dense oracle" begin
+ for k in (3, 4), model in (:onsite, :hopping), boundary in (:open, :periodic)
+ values = (1.2, 0.7, 0.4, 1.5)[1:k]
+ parameters = (; values, t=0.9, onsite=0.2, model, boundary)
+ Htn = TensorBinding.kbonacci_hamiltonian(k, 4; parameters...)
+ Hdense = TensorBinding._dense_kbonacci_hamiltonian(k, 4; parameters...)
+ @test Htn.N == TensorBinding.kbonacci_site_count(k, 4)
+ @test Htn.N == (k == 3 ? 13 : 15)
+ @test Htn.position_space isa TensorBinding.KBonacciPositionSpace
+ @test Htn.position_space.k == k
+ @test TensorBinding.ambient_dimension(Htn) == big(2)^4
+ @test TensorBinding.ambient_dimension(Htn) isa BigInt
+ @test length(Htn.sites) == 4 && all(s -> dim(s) == 2, Htn.sites)
+ @test Htn.scale > 0
+ matrix = _projected_kb_matrix(Htn)
+ @test maximum(abs.(matrix .- Hdense)) < 1e-11
+ @test norm(matrix - matrix') < 1e-11
+ last_letter = _kb_trailing_ones(TensorBinding.kbonacci_digits(k, Htn.N - 1, 4)) + 1
+ expected_wrap = boundary === :open ? 0.0 :
+ (model === :onsite ? parameters.t : values[last_letter])
+ @test matrix[end, 1] ≈ expected_wrap atol=1e-11
+ end
+
+ # k > L: every binary string is admissible and the chain has 2^L sites.
+ Hall = TensorBinding.kbonacci_hamiltonian(5, 3; values=(1, 2, 3, 4, 5))
+ @test Hall.N == 8
+ @test maximum(abs.(_projected_kb_matrix(Hall) .-
+ TensorBinding._dense_kbonacci_hamiltonian(5, 3; values=(1, 2, 3, 4, 5)))) < 1e-11
+
+ Hcomplex = TensorBinding.kbonacci_hamiltonian(
+ 3, 4; values=(1 + 0.2im, 2 - 0.1im, 0.5im), model=:hopping,
+ )
+ complex_matrix = _projected_kb_matrix(Hcomplex)
+ @test norm(complex_matrix - complex_matrix') < 1e-11
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(
+ 3, 4; values=(1 + 1im, 2.0, 3.0), model=:onsite,
+ )
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(1, 4; values=(1.0,))
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(3, 1; values=(1.0, 2.0, 3.0))
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(3, 4; values=(1.0, 2.0))
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(
+ 3, 4; values=(1.0, 2.0, 3.0), boundary=:twisted,
+ )
+ @test_throws ArgumentError TensorBinding.kbonacci_hamiltonian(
+ 3, 4; values=(1.0, 2.0, 3.0), model=:mixed,
+ )
+ end
+
+ # Tribonacci hopping chain with the paper's t_A/t_B = t_B/t_C = 0.8, t_C = 1:
+ # 24 physical sites inside a 32-state register.
+ H = TensorBinding.kbonacci_hamiltonian(3, 5; values=(0.64, 0.8, 1.0))
+ P = TensorBinding.physical_projector(H)
+ @test H.N == 24
+ @test real(tr(P)) ≈ H.N atol=1e-12
+ @test norm(apply(P, P; cutoff=1e-13) - P) / norm(P) < 1e-12
+ @test H.center == 0
+ @test occursin("TBHamiltonian", sprint(show, H))
+
+ @testset "projector-aware KPM" begin
+ Ncheb = 8
+ Tn, _, _ = TensorBinding.KPM_Tn(H, Ncheb; maxdim=100, cutoff=1e-12)
+ @test real(tr(Tn[1])) ≈ H.N atol=1e-10
+ @test norm(Tn[1] - P) < 1e-12
+
+ dense = TensorBinding._dense_kbonacci_hamiltonian(3, 5; values=(0.64, 0.8, 1.0))
+ decomposition = eigen(Hermitian(dense))
+ # Hopping model, zero onsite, even N: the spectrum is exactly chiral.
+ @test maximum(abs.(decomposition.values .+ reverse(decomposition.values))) < 1e-12
+ dense_moments = _kb_dense_kpm_moments(decomposition, Ncheb, H.center, H.scale)
+ @test maximum(abs.(real.(tr.(Tn[1:Ncheb])) .- dense_moments)) < 1e-8
+
+ energies = collect(range(-1.9, 1.9; length=7))
+ dos_tn = TensorBinding.get_dos_trace(H, Ncheb, energies; maxdim=100, cutoff=1e-12)
+ dos_dense = [
+ TensorBinding.get_ldos_from_mun(
+ dense_moments, Ncheb, (energy - H.center) / H.scale,
+ ) for energy in energies
+ ]
+ @test maximum(abs.(dos_tn .- dos_dense)) < 1e-8
+
+ ldos_dense = zeros(length(energies), H.N)
+ for site in 1:H.N
+ moments = _kb_dense_kpm_moments(decomposition, Ncheb, H.center, H.scale; site)
+ for (i, energy) in pairs(energies)
+ ldos_dense[i, site] = TensorBinding.get_ldos_from_mun(
+ moments, Ncheb, (energy - H.center) / H.scale,
+ )
+ end
+ end
+ for mode in (:mps, :mpo)
+ ldos_tn = TensorBinding.get_ldos_spatial(
+ H, Ncheb, energies; mode, maxdim=100, cutoff=1e-12,
+ )
+ @test size(ldos_tn) == (length(energies), H.N)
+ @test maximum(abs.(ldos_tn .- ldos_dense)) < 1e-8
+ end
+ probe = 7
+ ldos_online = TensorBinding.get_ldos_online(
+ H, Ncheb, probe, energies; maxdim=100, cutoff=1e-12,
+ )
+ @test maximum(abs.(ldos_online .- ldos_dense[:, probe])) < 1e-8
+
+ dos_stochastic = TensorBinding.get_dos_stochastic(
+ H, 4, [0.0]; N_sample=60, seed=7, maxdim=60, cutoff=1e-12,
+ )
+ dos_exact = TensorBinding.get_dos_trace(H, 4, [0.0]; maxdim=60, cutoff=1e-12)
+ @test all(isfinite, dos_stochastic)
+ @test abs(dos_stochastic[1] - dos_exact[1]) < 0.35 * max(abs(dos_exact[1]), 1.0)
+ end
+
+ @testset "interface and guards" begin
+ Hletters = TensorBinding.get_Hamiltonian("kbonacci", (A=0.64, B=0.8, C=1.0); L=5, k=3)
+ Hvalues = TensorBinding.get_Hamiltonian("kbonacci", (values=(0.64, 0.8, 1.0),); L=5, k=3)
+ @test Hletters.N == H.N && Hvalues.N == H.N
+ @test Hletters.position_space isa TensorBinding.KBonacciPositionSpace
+ @test Hletters.position_space.k == 3
+ @test maximum(abs.(_projected_kb_matrix(Hletters) .- _projected_kb_matrix(Hvalues))) < 1e-11
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8, C=1.0); L=5,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8); L=5, k=3,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8, C=1.0, values=(1, 2, 3)); L=5, k=3,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8, C=1.0, foo=1); L=5, k=3,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "kbonacci", (A=0.64, B=0.8, C=1.0); L=5, k=3, ref_sites=siteinds("Qubit", 5),
+ )
+
+ @test TensorBinding.site_axis(H) == collect(0:(H.N - 1))
+ @test TensorBinding.site_permutation(H) == collect(1:H.N)
+ @test_throws ArgumentError TensorBinding.site_axis(H; ordering=:conumber)
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial(
+ H, 4, [0.0]; ordering=:conumber,
+ )
+ @test_throws ArgumentError TensorBinding.add_onsite!(H, 0.1)
+ @test_throws ArgumentError TensorBinding.add_hopping!(H, 0.1)
+ @test_throws ArgumentError TensorBinding.get_bands(H, 4, 1, [0.0])
+ @test_throws BoundsError TensorBinding.physical_site_state(H, H.N + 1)
+ for x in 1:H.N
+ psi = TensorBinding.physical_site_state(H, x)
+ @test norm(psi) ≈ 1
+ @test abs(inner(psi, apply(P, psi))) ≈ 1 atol=1e-12
+ end
+ end
+end
diff --git a/test/metallic_mean.jl b/test/metallic_mean.jl
new file mode 100644
index 0000000..1df078f
--- /dev/null
+++ b/test/metallic_mean.jl
@@ -0,0 +1,213 @@
+using Test
+using LinearAlgebra
+using TensorBinding
+using ITensors
+using ITensorMPS
+
+# Substitution fixed point A -> A^m B, B -> A: an oracle independent of the
+# numeration system used by the package.
+function _metallic_mean_word(m::Integer, n::Integer)
+ word = "A"
+ while length(word) < n
+ word = join(c == 'A' ? "A"^m * "B" : "A" for c in word)
+ end
+ return word[1:n]
+end
+
+function _projected_mm_matrix(H)
+ basis = [TensorBinding.physical_site_state(H, x) for x in 1:H.N]
+ matrix = zeros(ComplexF64, H.N, H.N)
+ for j in eachindex(basis)
+ Hket = apply(H.mpo, basis[j]; cutoff=1e-13, maxdim=300)
+ for i in eachindex(basis)
+ matrix[i, j] = inner(basis[i], Hket)
+ end
+ end
+ return matrix
+end
+
+function _mm_dense_kpm_moments(decomposition, Ncheb, center, scale; site=nothing)
+ scaled = clamp.((decomposition.values .- center) ./ scale, -1.0, 1.0)
+ angles = acos.(scaled)
+ weights = isnothing(site) ? ones(length(scaled)) :
+ abs2.(decomposition.vectors[site, :])
+ return [sum(weights .* cos.(n .* angles)) for n in 0:(Ncheb - 1)]
+end
+
+@testset "Metallic-mean projected position space" begin
+ @testset "numeration system and word" begin
+ for m in 1:4, L in 2:5
+ N = TensorBinding.metallic_mean_site_count(m, L)
+ @test N == m * TensorBinding.metallic_mean_site_count(m, L - 1) +
+ TensorBinding.metallic_mean_site_count(m, L - 2)
+ weights = [Int(TensorBinding.metallic_mean_number(m, L - p)) for p in 1:L]
+ strings = [TensorBinding.metallic_mean_digits(m, n, L) for n in 0:(N - 1)]
+ @test all(d -> all(0 .<= d .<= m), strings)
+ @test all(d -> all(i -> !(d[i] == m && d[i + 1] != 0), 1:(L - 1)), strings)
+ @test [sum(d .* weights) for d in strings] == collect(0:(N - 1))
+ @test issorted(strings)
+ brute = [digits_vec for s in 0:((m + 1)^L - 1)
+ for digits_vec in ([(s ÷ (m + 1)^(L - p)) % (m + 1) for p in 1:L],)
+ if all(i -> !(digits_vec[i] == m && digits_vec[i + 1] != 0), 1:(L - 1))]
+ @test sort(brute) == strings
+
+ letters = [d[end] == m ? 'B' : 'A' for d in strings]
+ @test join(letters) == _metallic_mean_word(m, N)
+ @test [TensorBinding.metallic_mean_bond_symbol(m, L, b) for b in 1:N] ==
+ [c == 'A' ? :A : :B for c in letters]
+ end
+ @test_throws ArgumentError TensorBinding.metallic_mean_digits(2, 17, 3)
+ @test_throws ArgumentError TensorBinding.metallic_mean_number(0, 3)
+ @test_throws BoundsError TensorBinding.metallic_mean_bond_symbol(2, 3, 18)
+ end
+
+ @testset "m = 1 reproduces the Fibonacci chain" begin
+ for L in 2:6
+ N = TensorBinding.fibonacci_site_count(L)
+ @test TensorBinding.metallic_mean_site_count(1, L) == N
+ @test all(n -> TensorBinding.metallic_mean_digits(1, n, L) ==
+ TensorBinding.fibonacci_zeckendorf_digits(n, L), 0:(N - 1))
+ end
+ params = (; A=1.2, B=0.7, t=0.9, onsite=0.2)
+ for model in (:onsite, :hopping), boundary in (:open, :periodic)
+ dense_mm = TensorBinding._dense_metallic_mean_hamiltonian(
+ 1, 4; params..., model, boundary,
+ )
+ dense_fib = TensorBinding._dense_fibonacci_hamiltonian(
+ 4; params..., model, boundary,
+ )
+ @test dense_mm == dense_fib
+ Hmm = TensorBinding.metallic_mean_hamiltonian(1, 4; params..., model, boundary)
+ @test Hmm.N == 8
+ @test all(s -> dim(s) == 2, Hmm.sites)
+ @test maximum(abs.(_projected_mm_matrix(Hmm) .- dense_fib)) < 1e-11
+ end
+ end
+
+ @testset "TN construction matches the dense oracle" begin
+ for m in (2, 3), model in (:onsite, :hopping), boundary in (:open, :periodic)
+ parameters = (; A=1.2, B=0.7, t=0.9, onsite=0.2, model, boundary)
+ Htn = TensorBinding.metallic_mean_hamiltonian(m, 3; parameters...)
+ Hdense = TensorBinding._dense_metallic_mean_hamiltonian(m, 3; parameters...)
+ @test Htn.N == TensorBinding.metallic_mean_site_count(m, 3)
+ @test Htn.position_space isa TensorBinding.MetallicMeanPositionSpace
+ @test Htn.position_space.m == m
+ @test TensorBinding.ambient_dimension(Htn) == big(m + 1)^3
+ @test TensorBinding.ambient_dimension(Htn) isa BigInt
+ @test length(Htn.sites) == 3 && all(s -> dim(s) == m + 1, Htn.sites)
+ @test Htn.scale > 0
+ matrix = _projected_mm_matrix(Htn)
+ @test maximum(abs.(matrix .- Hdense)) < 1e-11
+ @test norm(matrix - matrix') < 1e-11
+ expected_wrap = boundary === :open ? 0.0 :
+ (model === :onsite ? parameters.t :
+ (TensorBinding.metallic_mean_digits(m, Htn.N - 1, 3)[end] == m ?
+ parameters.B : parameters.A))
+ @test matrix[end, 1] ≈ expected_wrap atol=1e-11
+ end
+
+ Hcomplex = TensorBinding.metallic_mean_hamiltonian(
+ 2, 3; A=1 + 0.2im, B=2 - 0.1im, model=:hopping,
+ )
+ complex_matrix = _projected_mm_matrix(Hcomplex)
+ @test norm(complex_matrix - complex_matrix') < 1e-11
+ @test_throws ArgumentError TensorBinding.metallic_mean_hamiltonian(
+ 2, 3; A=1 + 1im, B=2.0, model=:onsite,
+ )
+ @test_throws ArgumentError TensorBinding.metallic_mean_hamiltonian(0, 3; A=1.0, B=2.0)
+ @test_throws ArgumentError TensorBinding.metallic_mean_hamiltonian(2, 1; A=1.0, B=2.0)
+ @test_throws ArgumentError TensorBinding.metallic_mean_hamiltonian(
+ 2, 3; A=1.0, B=2.0, boundary=:twisted,
+ )
+ end
+
+ # Silver-mean hopping chain: 17 physical sites inside a 27-state register.
+ H = TensorBinding.metallic_mean_hamiltonian(2, 3; A=1.0, B=2.0)
+ P = TensorBinding.physical_projector(H)
+ @test H.N == 17
+ @test real(tr(P)) ≈ H.N atol=1e-12
+ @test norm(apply(P, P; cutoff=1e-13) - P) / norm(P) < 1e-12
+ @test H.center == 0
+ @test occursin("TBHamiltonian", sprint(show, H))
+
+ @testset "projector-aware KPM" begin
+ Ncheb = 8
+ Tn, _, _ = TensorBinding.KPM_Tn(H, Ncheb; maxdim=100, cutoff=1e-12)
+ @test real(tr(Tn[1])) ≈ H.N atol=1e-10
+ @test norm(Tn[1] - P) < 1e-12
+
+ dense = TensorBinding._dense_metallic_mean_hamiltonian(2, 3; A=1.0, B=2.0)
+ decomposition = eigen(Hermitian(dense))
+ dense_moments = _mm_dense_kpm_moments(decomposition, Ncheb, H.center, H.scale)
+ @test maximum(abs.(real.(tr.(Tn[1:Ncheb])) .- dense_moments)) < 1e-8
+
+ energies = collect(range(-3.5, 3.5; length=7))
+ dos_tn = TensorBinding.get_dos_trace(H, Ncheb, energies; maxdim=100, cutoff=1e-12)
+ dos_dense = [
+ TensorBinding.get_ldos_from_mun(
+ dense_moments, Ncheb, (energy - H.center) / H.scale,
+ ) for energy in energies
+ ]
+ @test maximum(abs.(dos_tn .- dos_dense)) < 1e-8
+
+ ldos_dense = zeros(length(energies), H.N)
+ for site in 1:H.N
+ moments = _mm_dense_kpm_moments(decomposition, Ncheb, H.center, H.scale; site)
+ for (i, energy) in pairs(energies)
+ ldos_dense[i, site] = TensorBinding.get_ldos_from_mun(
+ moments, Ncheb, (energy - H.center) / H.scale,
+ )
+ end
+ end
+ for mode in (:mps, :mpo)
+ ldos_tn = TensorBinding.get_ldos_spatial(
+ H, Ncheb, energies; mode, maxdim=100, cutoff=1e-12,
+ )
+ @test size(ldos_tn) == (length(energies), H.N)
+ @test maximum(abs.(ldos_tn .- ldos_dense)) < 1e-8
+ end
+ probe = 5
+ ldos_online = TensorBinding.get_ldos_online(
+ H, Ncheb, probe, energies; maxdim=100, cutoff=1e-12,
+ )
+ @test maximum(abs.(ldos_online .- ldos_dense[:, probe])) < 1e-8
+
+ dos_stochastic = TensorBinding.get_dos_stochastic(
+ H, 4, [0.0]; N_sample=60, seed=7, maxdim=60, cutoff=1e-12,
+ )
+ dos_exact = TensorBinding.get_dos_trace(H, 4, [0.0]; maxdim=60, cutoff=1e-12)
+ @test all(isfinite, dos_stochastic)
+ @test abs(dos_stochastic[1] - dos_exact[1]) < 0.35 * max(abs(dos_exact[1]), 1.0)
+ end
+
+ @testset "interface and guards" begin
+ Hd = TensorBinding.get_Hamiltonian("metallic_mean", (A=1.0, B=2.0); L=3, m=2)
+ @test Hd.N == H.N
+ @test Hd.position_space isa TensorBinding.MetallicMeanPositionSpace
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "metallic_mean", (A=1.0, B=2.0); L=3,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "metallic_mean", (A=1.0,); L=3, m=2,
+ )
+ @test_throws ArgumentError TensorBinding.get_Hamiltonian(
+ "metallic_mean", (A=1.0, B=2.0, foo=1); L=3, m=2,
+ )
+
+ @test TensorBinding.site_axis(H) == collect(0:(H.N - 1))
+ @test TensorBinding.site_permutation(H) == collect(1:H.N)
+ @test_throws ArgumentError TensorBinding.site_axis(H; ordering=:conumber)
+ @test_throws ArgumentError TensorBinding.get_ldos_spatial(
+ H, 4, [0.0]; ordering=:conumber,
+ )
+ @test_throws ArgumentError TensorBinding.add_onsite!(H, 0.1)
+ @test_throws ArgumentError TensorBinding.add_hopping!(H, 0.1)
+ @test_throws ArgumentError TensorBinding.get_bands(H, 4, 1, [0.0])
+ @test_throws BoundsError TensorBinding.physical_site_state(H, H.N + 1)
+ for x in 1:H.N
+ psi = TensorBinding.physical_site_state(H, x)
+ @test norm(psi) ≈ 1
+ @test abs(inner(psi, apply(P, psi))) ≈ 1 atol=1e-12
+ end
+ end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 591146f..c714e6a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -99,3 +99,9 @@ using TensorBinding: get_Hamiltonian, KPM_Tn, get_density_from_Tn,
end
end
+
+include("fibonacci.jl")
+include("fibonacci_sampling.jl")
+include("metallic_mean.jl")
+include("kbonacci.jl")
+include("gpu_mps_ldos.jl")