diff --git a/.gitignore b/.gitignore index f4dbd51..d7e741e 100644 --- a/.gitignore +++ b/.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 \ No newline at end of file +docs/src/bubble_lehmann.tex diff --git a/Project.toml b/Project.toml index 515e642..457a649 100644 --- a/Project.toml +++ b/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] diff --git a/logo_TensorBinding_v0.png b/logo_TensorBinding_v0.png new file mode 100644 index 0000000..7bc4070 Binary files /dev/null and b/logo_TensorBinding_v0.png differ diff --git a/src/2Dlattice_tk.jl b/src/2Dlattice_tk.jl index 864636a..afe8f78 100644 --- a/src/2Dlattice_tk.jl +++ b/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 diff --git a/src/KPM_tk.jl b/src/KPM_tk.jl index 0556ca3..cbb29e1 100644 --- a/src/KPM_tk.jl +++ b/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 diff --git a/src/QFT_tk.jl b/src/QFT_tk.jl index 9d3dff8..7c9df72 100644 --- a/src/QFT_tk.jl +++ b/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) # diff --git a/src/RPA_tk.jl b/src/RPA_tk.jl index 7c57fd3..83b5267 100644 --- a/src/RPA_tk.jl +++ b/src/RPA_tk.jl @@ -68,14 +68,14 @@ Replace the physical (bra + ket) indices of each site in `MPOin` with the corresponding index from `newsites`, preserving prime levels. """ function replace_sites(MPOin::MPO, newsites) - L = length(MPOin) + L = length(MPOin) indsMPO = siteinds(MPOin) T = MPO(L) for n in 1:L bra_old, ket_old = _bra_ket(indsMPO[n]) - T[n] = MPOin[n] * - delta(bra_old, prime(newsites[n])) * - delta(ket_old, newsites[n]) + T[n] = replaceinds(MPOin[n], + [bra_old, ket_old], + [prime(newsites[n]), newsites[n]]) end return T end @@ -576,10 +576,14 @@ function get_bubble_mpo(H1::TBHamiltonian, H2::TBHamiltonian, ω::Real; P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - verbose && println("Polarization bubble: computing P2...") - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + if H1 === H2 + P2 = P1 + else + verbose && println("Polarization bubble: computing P2...") + P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) + end # ---- Numerator: I₁⊗P₂ − P₁⊗I₂ ---- id1 = MPO(sites1, "Id") @@ -1497,20 +1501,28 @@ function get_bubble_mpo_cheb2d(H1::TBHamiltonian, H2::TBHamiltonian, Tn1, _, _ = KPM_Tn(H1.mpo, Ncheb, H1.sites; scale=scale1, center=center1, maxdim=maxdim, cutoff=cutoff, verbose=false) - verbose && println("cheb2d: building T_n(H2) moments...") - Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; - scale=scale2, center=center2, - maxdim=maxdim, cutoff=cutoff, verbose=false) + if H1 === H2 + Tn2 = Tn1 + else + verbose && println("cheb2d: building T_n(H2) moments...") + Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; + scale=scale2, center=center2, + maxdim=maxdim, cutoff=cutoff, verbose=false) + end N = length(Tn1) # = Ncheb + 1 (T_0 … T_Ncheb) verbose && println("cheb2d: computing P1...") P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - verbose && println("cheb2d: computing P2...") - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + if H1 === H2 + P2 = P1 + else + verbose && println("cheb2d: computing P2...") + P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) + end verbose && println("cheb2d: precomputing T_m(H1)·P1 and T_n(H2)·P2...") TP1 = [ITensorMPS.truncate!( @@ -1639,18 +1651,22 @@ function get_bubble_mpo_cheb2d_tucker(H1::TBHamiltonian, H2::TBHamiltonian, Tn1, _, _ = KPM_Tn(H1.mpo, Ncheb, H1.sites; scale=scale1, center=center1, maxdim=maxdim, cutoff=cutoff, verbose=false) - Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; - scale=scale2, center=center2, - maxdim=maxdim, cutoff=cutoff, verbose=false) + if H1 === H2 + Tn2 = Tn1 + else + Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; + scale=scale2, center=center2, + maxdim=maxdim, cutoff=cutoff, verbose=false) + end N = length(Tn1) verbose && println("cheb2d_mpo_tucker: computing density matrices...") P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + P2 = H1 === H2 ? P1 : _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) out_sites = siteinds("Qubit", L) @@ -1716,8 +1732,11 @@ function get_bubble_mpo_cheb2d_tucker(H1::TBHamiltonian, H2::TBHamiltonian, had_A = _hadamard_mpo(A_tuck[s1], B_tuck[s2], out_sites; maxdim=maxdim, cutoff=cutoff) had_B = _hadamard_mpo(C_tuck[s1], E_tuck[s2], out_sites; maxdim=maxdim, cutoff=cutoff) D_tuck[s1, s2] = ITensorMPS.truncate!(+(had_A, -1 * had_B; maxdim=maxdim); cutoff=cutoff) - - verbose && println(" ($s1,$s2)/($r_m,$r_n) done") + if verbose + idx = (s1 - 1) * r_n + s2 + (idx % 10 == 0 || idx == r_m * r_n) && + println(" ($s1,$s2)/($r_m,$r_n) done [$idx/$(r_m*r_n)]") + end end # ── Per-ω accumulation: scalar × MPO additions only ────────────────────── @@ -1800,20 +1819,28 @@ function get_bubble_diag_cheb2d(H1::TBHamiltonian, H2::TBHamiltonian, Tn1, _, _ = KPM_Tn(H1.mpo, Ncheb, H1.sites; scale=scale1, center=center1, maxdim=maxdim, cutoff=cutoff, verbose=false) - verbose && println("cheb2d_diag: building T_n(H2) moments...") - Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; - scale=scale2, center=center2, - maxdim=maxdim, cutoff=cutoff, verbose=false) + if H1 === H2 + Tn2 = Tn1 + else + verbose && println("cheb2d_diag: building T_n(H2) moments...") + Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; + scale=scale2, center=center2, + maxdim=maxdim, cutoff=cutoff, verbose=false) + end N = length(Tn1) verbose && println("cheb2d_diag: computing P1...") P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - verbose && println("cheb2d_diag: computing P2...") - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + if H1 === H2 + P2 = P1 + else + verbose && println("cheb2d_diag: computing P2...") + P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) + end verbose && println("cheb2d_diag: precomputing T_m(H1)·P1 and T_n(H2)·P2...") TP1 = [ITensorMPS.truncate!( @@ -1970,18 +1997,22 @@ function get_bubble_diag_cheb2d_svd(H1::TBHamiltonian, H2::TBHamiltonian, Tn1, _, _ = KPM_Tn(H1.mpo, Ncheb, H1.sites; scale=scale1, center=center1, maxdim=maxdim, cutoff=cutoff, verbose=false) - Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; - scale=scale2, center=center2, - maxdim=maxdim, cutoff=cutoff, verbose=false) + if H1 === H2 + Tn2 = Tn1 + else + Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; + scale=scale2, center=center2, + maxdim=maxdim, cutoff=cutoff, verbose=false) + end N = length(Tn1) verbose && println("cheb2d_diag_svd: computing density matrices...") P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + P2 = H1 === H2 ? P1 : _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) out_sites = siteinds("Qubit", L) @@ -2121,18 +2152,22 @@ function get_bubble_diag_cheb2d_tucker(H1::TBHamiltonian, H2::TBHamiltonian, Tn1, _, _ = KPM_Tn(H1.mpo, Ncheb, H1.sites; scale=scale1, center=center1, maxdim=maxdim, cutoff=cutoff, verbose=false) - Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; - scale=scale2, center=center2, - maxdim=maxdim, cutoff=cutoff, verbose=false) + if H1 === H2 + Tn2 = Tn1 + else + Tn2, _, _ = KPM_Tn(H2.mpo, Ncheb, H2.sites; + scale=scale2, center=center2, + maxdim=maxdim, cutoff=cutoff, verbose=false) + end N = length(Tn1) verbose && println("cheb2d_tucker: computing density matrices...") P1 = _get_density_matrix(H1, ϵF, P_method, Ncheb, maxdim, cutoff, purify_method, purify_maxdim, purify_maxiters, purify_tol, verbose) - P2 = _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, - purify_method, purify_maxdim, purify_maxiters, - purify_tol, verbose) + P2 = H1 === H2 ? P1 : _get_density_matrix(H2, ϵF, P_method, Ncheb, maxdim, cutoff, + purify_method, purify_maxdim, purify_maxiters, + purify_tol, verbose) out_sites = siteinds("Qubit", L) @@ -2204,8 +2239,11 @@ function get_bubble_diag_cheb2d_tucker(H1::TBHamiltonian, H2::TBHamiltonian, D_phys = replace_sites(D, H1.sites) D_k = conjugate_by_qft(D_phys; tol=qft_tol, maxdim=qft_maxdim) diag_D[s1, s2] = ITensorMPS.truncate!(extract_diagonal_to_mps(D_k); cutoff=cutoff) - - verbose && println(" ($s1,$s2)/($r_m,$r_n) done") + if verbose + idx = (s1 - 1) * r_n + s2 + (idx % 10 == 0 || idx == r_m * r_n) && + println(" ($s1,$s2)/($r_m,$r_n) done [$idx/$(r_m*r_n)]") + end end # ── Accumulate per ω: scalar × MPS additions only ──────────────────────── diff --git a/src/TensorBinding.jl b/src/TensorBinding.jl index 660f4fc..c8c9fb8 100644 --- a/src/TensorBinding.jl +++ b/src/TensorBinding.jl @@ -8,12 +8,8 @@ using ITensorMPS using Quantics using QuanticsTCI using QuanticsGrids -using TCIITensorConversion using TensorCrossInterpolation import TensorCrossInterpolation as TCI -using PyCall -using PyPlot -using Plots using FFTW using Base.Threads diff --git a/src/Topology_tk.jl b/src/Topology_tk.jl index 8f973b5..e7714e6 100644 --- a/src/Topology_tk.jl +++ b/src/Topology_tk.jl @@ -565,11 +565,13 @@ function get_C(H::TBHamiltonian, xfunc=nothing, yfunc=nothing; fermi::Real = 0.0, l = nothing, Λ::Real = 10, + Lambda = nothing, # ASCII alias for Λ Nchebychev::Int = 300, maxdim::Int = 500, cutoff::Float64 = 1e-8, Nel = nothing, quenched::Bool = true) + Λ_val = Lambda !== nothing ? Float64(Lambda) : Λ if xfunc === nothing || yfunc === nothing geom = H.geometry_uc !== nothing ? H.geometry_uc : H.geometry !== nothing ? H.geometry : @@ -580,7 +582,7 @@ function get_C(H::TBHamiltonian, xfunc=nothing, yfunc=nothing; P = _get_projector(H; method=method, fermi=fermi, Nchebychev=Nchebychev, maxdim=maxdim, cutoff=cutoff, Nel=Nel) return get_C_op_MPO_from_P(P, H.L, H.sites, xfunc, yfunc; - l=l, Λ=Λ, maxdim=maxdim, cutoff=cutoff, + l=l, Λ=Λ_val, maxdim=maxdim, cutoff=cutoff, quenched=quenched) end