Anouar

.gitignore

@@ -12,6 +12,10 @@ examples/local_testing3.ipynb
 src/RSI_tk.jl
 RhombGraph.ipynb
 logo_tensorBinding.svg
+examples/nontracked/
+GPU_tk.jl
+Conductivity_tk.jl
+TensorBinding_JOSS_AI.txt
 docs/src/joss.txt
 docs/src/oam_rabi.tex
-docs/src/bubble_lehmann.tex
+docs/src/bubble_lehmann.tex

Project.toml

@@ -10,13 +10,9 @@ ITensorMPS = "0d1a4710-d33b-49a5-8f18-73bdf49b47e2"
 ITensors = "9136182c-28ba-11e9-034c-db9fb085ebd5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NDTensors = "23ae76d9-e61a-49c4-8f12-3f1a16adf9cf"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
-PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Quantics = "87f76fb3-a40a-40c9-a63c-29fcfe7b7547"
 QuanticsGrids = "634c7f73-3e90-4749-a1bd-001b8efc642d"
 QuanticsTCI = "b11687fd-3a1c-4c41-97d0-998ab401d50e"
-TCIITensorConversion = "9f0aa9f4-9415-4e6a-8795-331ebf40aa04"
 TensorCrossInterpolation = "b261b2ec-6378-4871-b32e-9173bb050604"
 
 [compat]

src/2Dlattice_tk.jl

@@ -1573,172 +1573,4 @@ _geom_n_sub(::Val{:kagome})    = 3
 _geom_n_sub(::Val{:lieb})      = 3
 _geom_n_sub(::Val{:dice})      = 4
 
-"""
-    plot_ldos_2d(ldos_mat, ωlist, ω_target;
-                 geometry, Lx, Ly,
-                 markersize, colormap, colorbar, clims, title, kwargs...)
-        -> Plot
-
-Scatter plot of the spatial LDOS at the energy in `ωlist` nearest to `ω_target`.
-Each atom in the geometry is drawn as a coloured dot.
-
-**Arguments**
-
-- `ldos_mat` : `(Nω × n_atoms)` matrix as returned by `get_ldos_spatial` with
-  `num_x = H.N`.  Column `k` maps to atom `k` in the positions matrix of the
-  geometry (interleaved order `[A₀, B₀, …]`).  Pass the full-lattice matrix from
-  `get_ldos_spatial` directly — no sublattice filtering required:
-  - `proj_sl=k` result: only sublattice `k` atoms carry weight, others are zero.
-  - `proj_sl=nothing` result: every atom carries its own sublattice LDOS.
-- `geometry` : `:honeycomb`, `:kagome`, `:lieb`, or `:dice`.
-- `Lx`, `Ly` : log₂ grid dimensions (same values passed to the constructor).
-
-An assertion checks that `size(ldos_mat, 2) == n_atoms`; the error message
-reminds the user to call `get_ldos_spatial` with `num_x = H.N`.
-
-Examples
---------
-```julia
-# Single sublattice (only A atoms lit up, B atoms zero)
-p1 = plot_ldos_2d(ldos_A, ωlist, 0.5; geometry=:honeycomb, Lx=2, Ly=2)
-
-# All sublattices — every atom coloured by its own LDOS
-ldos_all = get_ldos_spatial(H_kg, 100, ωlist; num_x=H_kg.N)
-p  = plot_ldos_2d(ldos_all, ωlist, -2.0; geometry=:kagome, Lx=2, Ly=2,
-                  markersize=14, colormap=:plasma)
-```
-"""
-function plot_ldos_2d(ldos_mat::AbstractMatrix, ωlist, ω_target;
-                      geometry::Symbol = :honeycomb,
-                      Lx::Int,
-                      Ly::Int,
-                      markersize::Real = 12.0,
-                      colormap::Symbol = :inferno,
-                      colorbar::Bool   = true,
-                      clims            = nothing,
-                      title::String    = "",
-                      kwargs...)
-    positions_all = _geom_positions(Val(geometry), Lx, Ly)
-    n_atoms       = size(positions_all, 1)
-    n_cols        = size(ldos_mat, 2)
-    n_atoms == n_cols || error(
-        "plot_ldos_2d: ldos_mat has $n_cols columns but geometry has $n_atoms atoms. " *
-        "Call get_ldos_spatial with num_x=H.N so each column maps to one atom.")
-
-    ω_arr    = collect(ωlist)
-    ω_idx    = argmin(abs.(ω_arr .- ω_target))
-    ω_actual = ω_arr[ω_idx]
-    vals     = ldos_mat[ω_idx, :]   # one value per atom, already in positions order
-
-    cl   = isnothing(clims) ? (0.0, maximum(vals) + eps(Float64)) : clims
-    tstr = isempty(title)   ? "LDOS  ω ≈ $(round(ω_actual; digits=3))" : title
-
-    Plots.scatter(positions_all[:, 1], positions_all[:, 2];
-                  marker_z          = vals,
-                  color             = colormap,
-                  colorbar          = colorbar,
-                  clims             = cl,
-                  markersize        = markersize,
-                  markerstrokewidth = 0,
-                  xlabel            = "x",
-                  ylabel            = "y",
-                  aspect_ratio      = :equal,
-                  title             = tstr,
-                  label             = "",
-                  kwargs...)
-end
-
-
-"""
-    plot_ldos_multilayer(ldos_layers, ωlist, ω_target;
-                         stacking=:Bernal, geometry=:honeycomb, Lx, Ly,
-                         markersize, colormap, colorbar, clims, title, kwargs...)
-        -> Plot
-
-Scatter plot of the multilayer LDOS **as seen from above** at the energy
-nearest to `ω_target`.
-
-`ldos_layers` is a `Vector` of `(Nω × n_atoms)` matrices, one per layer, as
-returned by
-
-```julia
-ldos_layers = [get_ldos_spatial(H, Nc, ωlist;
-                   proj_layer=k, num_x=H.N) for k in 1:n_layers]
-```
-
-**Geometry**
-
-For `:AA` stacking all layers share the same 2D positions; the result is a
-single honeycomb with LDOS summed over all layers.
-
-For `:Bernal` stacking the odd and even layers form two physically distinct
-groups: odd layers (1, 3, …) sit at the standard honeycomb positions and even
-layers (2, 4, …) are displaced by the intra-cell A→B bond vector τ.  The
-function plots both groups together — a total of `2 × n_atoms` scatter points
-— giving a visual picture of the two interlocked honeycomb sublattices that
-make up the Bernal stack viewed from above.  LDOS is summed independently
-within each group (odd / even) without sublattice permutation, since the τ
-shift already places even-layer atoms at their correct visual registry.
-"""
-function plot_ldos_multilayer(ldos_layers::AbstractVector{<:AbstractMatrix},
-                               ωlist, ω_target;
-                               stacking::Symbol = :Bernal,
-                               geometry::Symbol = :honeycomb,
-                               Lx::Int,
-                               Ly::Int,
-                               markersize::Real = 10.0,
-                               colormap::Symbol = :inferno,
-                               colorbar::Bool   = true,
-                               clims            = nothing,
-                               title::String    = "",
-                               kwargs...)
-    stacking ∈ (:AA, :Bernal) ||
-        error("plot_ldos_multilayer: unknown stacking :$stacking. Supported: :AA, :Bernal.")
-
-    rs   = _geom_positions(Val(geometry), Lx, Ly)   # (n_atoms, 2)
-    Nω   = size(ldos_layers[1], 1)
-    ω_arr  = collect(ωlist)
-    ω_idx  = argmin(abs.(ω_arr .- ω_target))
-    ω_actual = ω_arr[ω_idx]
-
-    if stacking === :AA
-        # All layers share the same positions: simple sum.
-        vals = sum(ldos_k[ω_idx, :] for ldos_k in ldos_layers)
-        cl   = isnothing(clims) ? (0.0, maximum(vals) + eps(Float64)) : clims
-        tstr = isempty(title) ? "LDOS (AA top view)  ω ≈ $(round(ω_actual; digits=3))" : title
-        return Plots.scatter(rs[:, 1], rs[:, 2];
-                             marker_z=vals, color=colormap, clims=cl,
-                             colorbar=colorbar, markersize=markersize,
-                             markerstrokewidth=0, aspect_ratio=:equal,
-                             xlabel="x", ylabel="y", title=tstr, label="",
-                             kwargs...)
-    end
-
-    # ── Bernal: two groups of layers at two shifted honeycomb lattices ────────
-    # τ = intra-cell A→B bond vector (first unit cell).
-    τ = rs[2, :] - rs[1, :]          # (2,) displacement
-
-    vals_odd  = zeros(size(rs, 1))   # group 1: layers 1, 3, 5, …
-    vals_even = zeros(size(rs, 1))   # group 2: layers 2, 4, …
-    for (k, ldos_k) in enumerate(ldos_layers)
-        if isodd(k)
-            vals_odd  .+= ldos_k[ω_idx, :]
-        else
-            vals_even .+= ldos_k[ω_idx, :]
-        end
-    end
-
-    rs_even  = rs .+ τ'              # even-layer honeycomb, shifted by τ
-    all_pos  = vcat(rs, rs_even)     # (2*n_atoms, 2)
-    all_vals = vcat(vals_odd, vals_even)
-
-    cl   = isnothing(clims) ? (0.0, maximum(all_vals) + eps(Float64)) : clims
-    tstr = isempty(title) ? "LDOS (Bernal top view)  ω ≈ $(round(ω_actual; digits=3))" : title
-
-    Plots.scatter(all_pos[:, 1], all_pos[:, 2];
-                  marker_z=all_vals, color=colormap, clims=cl,
-                  colorbar=colorbar, markersize=markersize,
-                  markerstrokewidth=0, aspect_ratio=:equal,
-                  xlabel="x", ylabel="y", title=tstr, label="",
-                  kwargs...)
-end
+# plot_ldos_2d and plot_ldos_multilayer moved to examples/nontracked/TensorBinding_plots.jl

src/KPM_tk.jl

@@ -112,7 +112,7 @@ function KPM_Tn(H_mpo::MPO, N::Int, sites;
         if verbose
             if k%5 == 0 || k == N+1 # print info every 5 iterations and at the end
                 println("Computed T_$((k-1)) with maxlinkdim = ", ITensorMPS.maxlinkdim(T_k))
-            end 
+            end
         end
     end
 

src/QFT_tk.jl

@@ -108,6 +108,84 @@ function conjugate_by_qft(W; tol=1e-9, maxdim::Int=100)
 end
 
 
+"""
+    conjugate_by_qft(H::TBHamiltonian, W::MPO; tol, maxdim) -> MPO
+
+TBHamiltonian-aware version of `conjugate_by_qft`.  Applies `U·W·U†` where
+`U` is the QFT acting **only on the position (Qubit) sites** of `H`, with
+identity operators at all auxiliary sites (Layer, spin, sublattice, Nambu).
+
+Use this overload whenever `W` lives on the full `H.sites` space (including
+aux indices), as is the case in the bubble pipeline after `replace_sites`.
+"""
+function conjugate_by_qft(H::TBHamiltonian, W::MPO; tol=1e-9, maxdim::Int=100)
+    pos_s = _pos_sites(H)
+    R     = length(pos_s)
+    FTirev_pos = fix_sites(MPO(TCI.reverse(QuanticsTCI.quanticsfouriermpo(R; sign=-1.0, normalize=true))), pos_s)
+    FTrev_pos  = fix_sites(MPO(TCI.reverse(QuanticsTCI.quanticsfouriermpo(R; sign=+1.0, normalize=true))), pos_s)
+    FTirev = _embed_in_full_sites(H, FTirev_pos)
+    FTrev  = _embed_in_full_sites(H, FTrev_pos)
+    Op1    = apply(W,                        FTirev; cutoff=tol, maxdim=maxdim)
+    Op2    = apply(swapprime(FTrev, 0 => 1), Op1;   cutoff=tol, maxdim=maxdim)
+    return TCI.truncate(Op2; cutoff=tol, maxdim=maxdim)
+end
+
+
+"""
+    _embed_in_full_sites(H, mpo_pos) -> MPO
+
+Embed `mpo_pos` (which lives on `_pos_sites(H)`) into the full `H.sites`
+space by prepending/appending dim-1 identity tensors at each auxiliary site.
+"""
+function _embed_in_full_sites(H::TBHamiltonian, mpo_pos::MPO)
+    pos_set  = Set(_pos_sites(H))
+    first_pos = findfirst(s -> s ∈ pos_set, H.sites)
+    last_pos  = findlast( s -> s ∈ pos_set, H.sites)
+    pre_aux  = H.sites[1:first_pos-1]
+    post_aux = H.sites[last_pos+1:end]
+    result   = mpo_pos
+    for s in reverse(pre_aux)
+        result = mpo_kron(MPO([dense(delta(s, prime(s)))]), result)
+    end
+    for s in post_aux
+        result = mpo_kron(result, MPO([dense(delta(s, prime(s)))]))
+    end
+    return result
+end
+
+
+"""
+    _embed_displacement_in_full_sites(H, mpo_pos) -> MPO
+
+Like `_embed_in_full_sites` but pads auxiliary sites with all-ones matrices
+instead of identity.  Required for current-operator construction: the
+displacement (xᵣ − xᵣ′) depends only on position, so it must be broadcast
+uniformly across all auxiliary (sublattice, layer, spin, Nambu) index pairs,
+including off-diagonal ones where physical hoppings exist.
+
+Using identity at an aux site with off-diagonal hoppings (e.g. sublattice A↔B
+in honeycomb, or inter-layer tunneling in bilayers) would set those current
+matrix elements to zero and give σ = 0.
+"""
+function _embed_displacement_in_full_sites(H::TBHamiltonian, mpo_pos::MPO)
+    pos_set   = Set(_pos_sites(H))
+    first_pos = findfirst(s -> s ∈ pos_set, H.sites)
+    last_pos  = findlast( s -> s ∈ pos_set, H.sites)
+    pre_aux   = H.sites[1:first_pos-1]
+    post_aux  = H.sites[last_pos+1:end]
+    result    = mpo_pos
+    for s in reverse(pre_aux)
+        ones_t = dense(ITensor(ones(Float64, dim(s), dim(s)), prime(s), s))
+        result = mpo_kron(MPO([ones_t]), result)
+    end
+    for s in post_aux
+        ones_t = dense(ITensor(ones(Float64, dim(s), dim(s)), prime(s), s))
+        result = mpo_kron(result, MPO([ones_t]))
+    end
+    return result
+end
+
+
 # ============================================================
 # 2. Legacy sublattice projection  (mask sandwich)
 #