GPU

KPM_Tn_gpu(H_mpo, N, sites; scale, center, maxdim, cutoff, keep_indices, verbose)
    -> (Tn_list, scale, center)

GPU version of KPM_Tn. Moves the identity and scaled Hamiltonian MPOs to GPU (ComplexF32) before the recurrence so all Tn tensors stay on GPU. Requires using CUDA.

If keep_indices is provided (a Set{Int}, 1-based into the returned vector where index 1 = T0, 2 = T1, …), only those Tns are retained in memory. All other slots are set to nothing. The recurrence itself always runs to completion — keep_indices only controls which results are stored.

density_profile_from_dm_gpu(density_mpo, sites=nothing; mode=:direct) -> MPS

GPU-resident analogue of density_profile_from_dm. If density_mpo is a CPU MPO it is uploaded once; if it is already on GPU it is used in place. The returned profile is a GPU MPS. mode=:complement returns 1 - diag(D) on GPU.

get_C_gpu(H::TBHamiltonian, xfunc=nothing, yfunc=nothing; kwargs...) -> Function

GPU-accelerated real-space Chern marker. Mirrors get_C exactly but runs all MPO×MPO products (projector assembly and C1–C4 construction) on GPU.

Returns the same closure C_at(uc::Int) -> ComplexF64 as get_C.

Key differences from get_C

• All apply/truncate! operations run on GPU tensors.
• dtype (default ComplexF32) selects the GPU element type; the marker is intrinsically complex, so only ComplexF32 / ComplexF64 are accepted. Use dtype=ComplexF64 to avoid NaN from ComplexF32 eigendecompositions on large systems at tight cutoffs (a warning is emitted for ComplexF32 + cutoff < 1e-6).
• The projector is built on CPU first (via _get_projector), then moved to GPU. For method=:mcweeny this means the purification loop runs on GPU.
• sequential mode is not supported (non-sequential quenched is always used).

All other keyword arguments are identical to get_C.

get_bands_gpu(H, Ncheb, ω_phys_vals; kwargs...)
    -> Matrix{Float64}  or  NamedTuple(Ak, ticks, labels)

GPU-accelerated version of get_bands.

GPU handles: the full Chebyshev MPO recurrence (the dominant cost) and the QFT sandwich applied to each Chebyshev moment. CPU handles: k-group setup, KPM weight matrix, final scalar accumulation.

Use type=ComplexF32 or type=ComplexF64 to choose the GPU tensor datatype. dtype=... is accepted as an alias for consistency with the non-Hermitian GPU entry points. ComplexF32 is faster, while ComplexF64 is safer at tight cutoffs on large systems.

All keyword arguments are identical to the TBHamiltonian overload of get_bands. The return value is also identical: a plain Matrix{Float64} when no kpath is given, or a NamedTuple(Ak, ticks, labels) when a high-symmetry path is requested.

Usage:

using CUDA
res = TensorBinding.get_bands_gpu(H, 500, omega;
        kpath=[:G, :M, :Kp, :G], kpath_lattice=:honeycomb,
        num_x=50, maxdim=200, type=ComplexF64, printinfo=true)
heatmap(1:size(res.Ak,2), omega, res.Ak; xticks=(res.ticks, res.labels))

get_dos_stochastic_gpu(H, Ncheb, ω_phys_vals; kwargs...)
    -> Vector{Float64}   length Nω

GPU-accelerated stochastic density of states via MPS Chebyshev KPM.

For each random sample the scaled Hamiltonian MPO lives on GPU and the product- state MPS is transferred to GPU once before the recursion starts. The Chebyshev moments ⟨ψ₀|T_n(H̃)|ψ₀⟩ are scalars pulled to CPU at each step.

Signature and optional kwargs are identical to get_dos_stochastic (CPU). N_bound (exciton bound-sector enrichment) is supported. Use dos_weighting=:sample to return the unweighted sampled signal avg_full + avg_bound, which is useful when visualising exciton peaks that are otherwise hidden by continuum phase-space factors in the trace DOS. For exciton Hamiltonians, continuum_only=true samples ordered electron-hole product states with x_e != x_h for the N_sample branch.

kernel=:hodc selects the Higher-Order Delta Chebyshev reconstruction (eta, m_order control the contour); its weights already carry the full KPM normalisation, so no √(1−ω²) denominator is applied. eta=0 falls back to 1/(Ncheb+1). Otherwise kernel is a convolution kernel (:jackson default, :lorentz with lambda, :fejer, :dirichlet).

get_exciton_cheb_convergence_gpu(H, X, Ncheb_max;
                                  maxdim_test, maxdim_ref, cutoff, printinfo)
    -> NamedTuple

Run two parallel GPU Chebyshev KPM recursions starting from |X,X⟩ = mpsexciton(X, H.sites): a reference recursion at maxdim_ref and a test recursion at maxdim_test. At each step n the two Chebyshev vectors φref^n and φtest^n are compared to yield:

errfidelityn = 1 − |⟨φref^n | φtest^n⟩|² / (‖φref^n‖² ‖φtest^n‖²)

which grows from 0 (perfect agreement) towards 1 as truncation errors accumulate.

Returns a NamedTuple with Float64 / Int vectors of length Nchebmax: n — Chebyshev index (1-based) muref — KPM moment ⟨X,X|Tn(H̃)|X,X⟩ from reference recursion mutest — KPM moment from test recursion deltamu — |muref − mutest| (moment discrepancy) errfidelity — infidelity 1 − fidelity as above mdimref — maxlinkdim of φref^n mdimtest — maxlinkdim of φtest^n normref — ‖φref^n‖ (should stay ≤ 1; growth indicates instability) normtest — ‖φtest^n‖

The Hamiltonian is rescaled internally: H̃ = (H − center·I) / scale, same as in getexcitonldosspatialgpu. All MPS live on GPU in ComplexF32 throughout.

Use this to find the critical Ncheb beyond which maxdim_test is too small for a given system size — the threshold is where err_fidelity departs significantly from 0 (e.g. > 0.01 for a tight criterion, > 0.1 for a loose one).

get_exciton_ldos_spatial_gpu(H, Ncheb, ω_phys_vals;
                             Lx, num_y, reduce,
                             X_list, X_groups, num_x, num_avg, x_start, x_end,
                             kernel, lambda, eta, m_order, maxdim, cutoff,
                             verbose, printinfo)
    -> Matrix{Float64}   (Nω × n_cols)

GPU-accelerated spatial exciton LDOS A(X,ω) = ⟨X,X|δ(ω−H)|X,X⟩ via MPS Chebyshev KPM. One GPU Chebyshev recursion runs per probe position X (electron = hole = X, 1-indexed in 1:H.N) from |X,X⟩ = mpsexciton(X, H.sites); moments are scalars pulled to CPU.

1D sampling (default, Lx=nothing): num_x coarse positions over [x_start, x_end] with num_avg sub-positions per coarse cell averaged per output pixel.

2D grid (Lx provided): positions are 1-indexed on the (Lx+Ly)-qubit quantics grid, encoded as X = ix + iy·2^Lx + 1 (row-major, 0-indexed). num_x × num_y coarse cells are sampled via spatial_sampling_plan with num_avg sub-positions per cell. Output columns are row-major over coarse cells (iy outer, ix inner).

X_list / X_groups / x_groups bypass the automatic plan and pass positions directly.

kernel=:hodc selects HODC reconstruction (eta, m_order; eta=0 → 1/(Ncheb+1)). Other kernels: :jackson (default), :lorentz (lambda), :fejer, :dirichlet.

Use type=ComplexF32 (default, faster) or type=ComplexF64 (safer at tight cutoffs 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.

get_ldos_spatial_gpu(H, Ncheb, ω_phys_vals; kwargs...)
    -> Matrix{Float64}   shape (Nω × n_spatial_cols)

GPU-accelerated version of get_ldos_spatial (MPO mode only).

Sampling procedures (reduce) — see spatial_sampling_plan.

• :point (default) — read the LDOS at num_x[×num_y] cells / x_groups (optionally box-averaged). Coarse grids alias thin features.
• :block — integrate over num_x × num_y blocks (powers of two) by tracing out the within-block bits; gap-free, so thin in-gap edge channels on a large system cannot be missed. The scalable tool for large-scale edge-state maps.

Column layout

• No sublattice DOF: (Nω × ng), one column per pixel (group or block).
• Sublattice resolved: (Nω × ng×n_sub), interleaved [A₀, B₀, A₁, B₁, …].
• Sublattice averaged (large scale / :block): (Nω × ng), one value per pixel.

For :block, columns are row-major over coarse pixels (col = ixp + iyp·num_x + 1).

GPU/CPU split

GPU: entire Chebyshev MPO recurrence, aux projections, diagonal extraction, real-space scalar sampling (point eval or block integration). CPU: KPM weight matrix, output accumulation (scalars only).

Keyword arguments are identical to get_ldos_spatial (:mps mode is not available on GPU; only the single-pass MPO mode is implemented here). Pass type=ComplexF32 or type=ComplexF64 to choose the GPU tensor datatype consistently throughout the MPO recurrence, projections, and diagonal extraction.

Usage

using CUDA
# point map
ldos = TensorBinding.get_ldos_spatial_gpu(H, 200, ωlist;
    x_groups = [[uc] for uc in 1:H.N], maxdim=200, printinfo=true)
# block-integrated large-scale edge-state map (num_x, num_y powers of two)
ldos = TensorBinding.get_ldos_spatial_gpu(H, 200, ωlist;
    reduce=:block, num_x=128, num_y=128, sublattice=:average, maxdim=200)

get_nh_density_trajectory_gpu(H, rho0; nsteps, dt, sample_every, ...)
    -> (density, times, centers, groups, maxlinkdims)

GPU RK4 evolution of a density matrix under a static non-Hermitian Hamiltonian using d rho/dt = -i(H rho - rho Hdagger). H may be a TBHamiltonian or an MPO; rho0 may be a CPU or GPU MPO. The Hamiltonian and density matrix are uploaded once and the RK4 loop stays on GPU. At every sampled time, the diagonal of rho(t) is extracted on GPU and only scalar values at the requested groups are copied back to CPU.

Sampling follows the 1D spatial_sampling_plan convention: use num_x=0 to sample all sites, or set num_x to a smaller number for coarse production output. num_avg > 1 averages a few sub-points per sampled spatial bin in :point mode. With reduce=:block, num_x must be a power of two and each output value is the GPU block average over a contiguous interval of size 2^L / num_x.

get_nh_dos_grid_diag_trace_gpu(H, xlims, nx, ylims, ny, n; scale=nothing,
                               nh_scale_padding=1.05, ...)
    -> (xgrid, ygrid, Z)

Grid companion to get_nh_dos_points_diag_trace_gpu.

get_nh_dos_grid_gpu(H, xlims, nx, ylims, ny, n; scale=nothing,
                    nh_scale_padding=1.05, n_random, ...)
    -> (xgrid, ygrid, Z)

GPU stochastic non-Hermitian KPM spectral-weight grid. This mirrors nh_spectrum_grid(...; mode=:stochastic): for each complex point z = x + im*y, the non-Hermitian Hamiltonian is hermitized on CPU, then the dual-chain stochastic MPS recurrence runs on GPU. Only scalar moments are copied back to CPU.

The integer n follows the existing NH convention: the partial recurrence runs to order 2n.

get_nh_dos_points_diag_trace_gpu(H, z_points, n; scale=nothing,
                                 nh_scale_padding=1.05, point_ids, ...)
    -> Vector{Float64}

Deterministic GPU NH KPM at an explicit list of complex energies. For each z, the hermitized NH problem is built on CPU, then the online MPO-MPO recurrence runs on GPU and evaluates the total trace through diagonal extraction plus a GPU-resident all-sites sum. This avoids stochastic probes.

The integer n follows the NH convention used elsewhere in this file: the partial recurrence runs to order 2n.

get_nh_dos_points_gpu(H, z_points, n; scale=nothing,
                      nh_scale_padding=1.05, n_random, point_ids, ...)
    -> Vector{Float64}

GPU stochastic NH KPM at an explicit list of complex energies. This is the array-job companion to get_nh_dos_grid_gpu: each z_points[j] is independent, so production scripts can split a large grid over many GPUs and concatenate the long-form CSV outputs afterward.

When seed is an integer, each point uses a deterministic seed seed + seed_stride * point_id, where point_id defaults to the local point index. Supplying global flattened grid indices as point_ids makes stochastic samples reproducible independent of how the grid is tiled.

get_scf_bands_gpu(res, Ncheb, omega; kwargs...) -> (Ak, omega, ticks, labels)

Compute spin-summed mean-field bands from a converged magnetic SCF result. This is deliberately separate from scf_magnetic_hubbard_gpu: it initializes from the CPU res.H_up/res.H_dn, then each get_bands_gpu call uploads once and keeps the Chebyshev/QFT accumulation on GPU, extracting only scalars.

get_scf_magnetization_gpu(res; kwargs...) -> (values, centers, groups, n_up, n_dn)

Sample the converged magnetic SCF density matrices on GPU and extract only the final scalar values. If res carries GPU density MPOs from scf_magnetic_hubbard_gpu, they are reused directly; otherwise the CPU density MPOs are uploaded once. Each sampled point is evaluated in the same big-endian real-space convention as binary_to_MPS.

Set reduce=:block to average over every unit cell in each coarse block by tracing the within-block position bits on GPU. In block mode, num_x and num_y must be powers of two.

get_state_amplitude_trajectory_gpu(H, psi0; nsteps, dt, sample_every, ...)
    -> (amplitude, times, centers, groups, norms, maxlinkdims)

GPU TDVP evolution of a single-particle MPS state under the physical Hamiltonian H. ITensorMPS.tdvp(operator, t, init) computes exp(t * operator) * init (generator form), so the operator passed to tdvp is -im * H, which implements the Schrödinger evolution dψ/dt = -im * Hψ; therefore a loss term -im * Γ (Γ >= 0) damps the norm, matching evolve_with_tdvp(H::TBHamiltonian,...) on CPU and the NH RK4 convention dρ/dt = -i(Hρ - ρH†). H may be a TBHamiltonian or an MPO. The Hamiltonian and initial state are uploaded once, the TDVP loop stays on GPU, and only sampled scalar amplitudes are copied back to CPU.

The returned amplitude matrix has rows = sampled positions and columns = sampled times. By default it stores real(<x|psi(t)>), matching panel (d) of APSOS_NH_testing. Set component=:imag, :abs, or :probability if needed.

Sampling follows spatial_sampling_plan in 1D. reduce=:point samples representative positions or explicit groups. reduce=:block returns the block-averaged complex amplitude over contiguous intervals.

scf_magnetic_hubbard_gpu(H0, U; kwargs...) -> NamedTuple

GPU-accelerated two-channel collinear magnetic mean-field loop for the on-site Hubbard model. The SCF iteration keeps the density profiles, Hartree MPOs, Hamiltonian MPOs, density matrices, RMS checks, and mixing on GPU; CPU objects are built only for initialization and for the compatibility fields returned at the end.

H_up = H0_up + U·diag(n_dn − background)
H_dn = H0_dn + U·diag(n_up   − background)

Only density_method=:mcweeny is supported here (grand-canonical at fermi); for particle-number-fixed SP2 use the CPU scf_magnetic_hubbard. A concrete purification scale is required so the GPU initial guess can be formed without estimating spectral bounds on CPU during the loop.

ComplexF32 eigen-decompositions can NaN at very tight cutoffs; a warning is emitted if cutoff < 1e-5, but the requested cutoff is used as-is.

Post-convergence observables are intentionally separate. Use get_scf_magnetization_gpu or get_scf_bands_gpu on the returned result when you want those GPU-accelerated diagnostics.