From 0f50dc55995bfc50f28f123a93088acece5f724a Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 10:24:38 +0200 Subject: [PATCH 1/9] Fix off-by-one error --- src/environments/qp_envs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index acc9fa9f3..e07f20369 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -46,7 +46,7 @@ function environments(qp::MultilineQP, operator::MultilineMPO, lenvs, renvs; kwa end function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; kwargs...) - ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1)) + ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1)) .+ 1 solver = environment_alg(exci, H, exci; kwargs...) AL = exci.left_gs.AL From 3fe94d2021541046336820aa9de48d1a2e5597cf Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 10:38:59 +0200 Subject: [PATCH 2/9] centralize regularization calls --- .../excitation/exci_transfer_system.jl | 10 +++--- src/environments/qp_envs.jl | 35 +++++++++---------- 2 files changed, 22 insertions(+), 23 deletions(-) diff --git a/src/algorithms/excitation/exci_transfer_system.jl b/src/algorithms/excitation/exci_transfer_system.jl index bb0b90095..0e6428841 100644 --- a/src/algorithms/excitation/exci_transfer_system.jl +++ b/src/algorithms/excitation/exci_transfer_system.jl @@ -59,8 +59,9 @@ function left_excitation_transfer_system( T = TransferMatrix(exci.right_gs.AR, H_partial, exci.left_gs.AL) start = scale!(last(found[1:i] * T), cis(-mom * len)) if exci.trivial && isidentitylevel(H, i) - @plansor start[-1 -2; -3 -4] -= start[1 4; -3 2] * r_RL(exci.right_gs)[2; 3] * - τ[3 4; 5 1] * l_RL(exci.right_gs)[-1; 6] * τ[5 6; -4 -2] + ρ_left = l_RL(exci.right_gs) + ρ_right = r_RL(exci.right_gs) + regularize!(start, ρ_right, ρ_left) end found[i] = add!(start, GBL[i]) @@ -147,8 +148,9 @@ function right_excitation_transfer_system( T = TransferMatrix(exci.left_gs.AL, H_partial, exci.right_gs.AR) start = scale!(first(T * found[i:odim]), cis(mom * len)) if exci.trivial && isidentitylevel(H, i) - @plansor start[-1 -2; -3 -4] -= τ[6 2; 3 4] * start[3 4; -3 5] * - l_LR(exci.right_gs)[5; 2] * r_LR(exci.right_gs)[-1; 1] * τ[-2 -4; 1 6] + ρ_left = l_LR(exci.right_gs) + ρ_right = r_LR(exci.right_gs) + regularize!(start, ρ_left, ρ_right) end found[i] = add!(start, GBR[i]) diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index e07f20369..483f45252 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -62,11 +62,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; lBs[pos + 1] += leftenv(lenvs, pos, exci.left_gs) * TransferMatrix(exci[pos], H[pos], AL[pos]) / cis(exci.momentum) - if exci.trivial # regularization of trivial excitations + if exci.trivial && !isempty(ids) # regularization of trivial excitations + ρ_left = l_RL(exci.left_gs, pos + 1) + ρ_right = r_RL(exci.left_gs, pos) for i in ids - @plansor lBs[pos + 1][i][-1 -2; -3 -4] -= lBs[pos + 1][i][1 4; -3 2] * - r_RL(exci.left_gs, pos)[2; 3] * τ[3 4; 5 1] * - l_RL(exci.left_gs, pos + 1)[ -1; 6 ] * τ[5 6; -4 -2] + regularize!(lBs[pos + 1][i], ρ_right, ρ_left) end end end @@ -78,12 +78,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; rBs[pos - 1] += TransferMatrix(exci[pos], H[pos], AR[pos]) * rightenv(renvs, pos, exci.right_gs) * cis(exci.momentum) - if exci.trivial + if exci.trivial && !isempty(ids) + ρ_left = l_LR(exci.left_gs, pos) + ρ_right = r_LR(exci.left_gs, pos - 1) for i in ids - ρ_left = l_LR(exci.left_gs, pos) - ρ_right = r_LR(exci.left_gs, pos - 1) - @plansor rBs[pos - 1][i][-1 -2; -3 -4] -= τ[6 4; 1 3] * - rBs[pos - 1][i][1 3; -3 2] * ρ_left[2; 4] * ρ_right[-1; 5] * τ[-2 -4; 5 6] + regularize!(rBs[pos - 1][i], ρ_left, ρ_right) end end end @@ -101,12 +100,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; for i in 1:(length(exci) - 1) lB_cur = lB_cur * TransferMatrix(AR[i], H[i], AL[i]) / cis(exci.momentum) - if exci.trivial + if exci.trivial && !isempty(ids) + ρ_left = l_RL(exci.left_gs, i + 1) + ρ_right = r_RL(exci.left_gs, i) for k in ids - ρ_left = l_RL(exci.left_gs, i + 1) - ρ_right = r_RL(exci.left_gs, i) - @plansor lB_cur[k][-1 -2; -3 -4] -= lB_cur[k][1 4; -3 2] * - ρ_right[2; 3] * τ[3 4; 5 1] * ρ_left[-1; 6] * τ[5 6; -4 -2] + regularize!(lB_cur[k], ρ_right, ρ_left) end end @@ -117,12 +115,11 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; for i in length(exci):-1:2 rB_cur = TransferMatrix(AL[i], H[i], AR[i]) * rB_cur * cis(exci.momentum) - if exci.trivial + if exci.trivial && !isempty(ids) + ρ_left = l_LR(exci.left_gs, i) + ρ_right = r_LR(exci.left_gs, i - 1) for k in ids - ρ_left = l_LR(exci.left_gs, i) - ρ_right = r_LR(exci.left_gs, i - 1) - @plansor rB_cur[k][-1 -2; -3 -4] -= τ[6 4; 1 3] * - rB_cur[k][1 3; -3 2] * ρ_left[2; 4] * ρ_right[-1; 5] * τ[-2 -4; 5 6] + regularize!(rB_cur[k], ρ_left, ρ_right) end end From fe444baed23a5b33d288c128e5030e7a5873271b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 10:48:57 +0200 Subject: [PATCH 3/9] avoid braiding in regularize --- src/transfermatrix/transfermatrix.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transfermatrix/transfermatrix.jl b/src/transfermatrix/transfermatrix.jl index 88030e91d..31a33ab06 100644 --- a/src/transfermatrix/transfermatrix.jl +++ b/src/transfermatrix/transfermatrix.jl @@ -87,6 +87,6 @@ function regularize!(v::MPOTensor, lvec::MPSTensor, rvec::MPSTensor) end function regularize!(v::MPOTensor, lvec::MPSBondTensor, rvec::MPSBondTensor) - return @plansor v[-1 -2; -3 -4] -= τ[6 2; 3 4] * v[3 4; -3 5] * lvec[5; 2] * rvec[-1; 1] * - τ[-2 -4; 1 6] + λ = @plansor lvec[2; 1] * removeunit(removeunit(v, 3), 2)[1; 2] + return add!(v, insertleftunit(insertrightunit(rvec, 1; dual = isdual(space(v, 2))), 3), -λ) end From f99504aeb5eec96572c503ffca541dd7a56592e1 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 10:51:38 +0200 Subject: [PATCH 4/9] remove unused code --- .../excitation/exci_transfer_system.jl | 86 ------------------- 1 file changed, 86 deletions(-) diff --git a/src/algorithms/excitation/exci_transfer_system.jl b/src/algorithms/excitation/exci_transfer_system.jl index 0e6428841..8a6e68dd0 100644 --- a/src/algorithms/excitation/exci_transfer_system.jl +++ b/src/algorithms/excitation/exci_transfer_system.jl @@ -1,47 +1,3 @@ -function left_excitation_transfer_system( - lBs, H, exci; mom = exci.momentum, - solver = Defaults.linearsolver - ) - len = length(H) - found = zero.(lBs) - odim = length(lBs) - - for i in 1:odim - # this operation can in principle be even further optimized for larger unit cells - # as we only require the terms that end at level i. - # this would require to check the finite state machine, and discard non-connected - # terms. - H_partial = map(site -> H.data[site, 1:i, 1:i], 1:len) - T = TransferMatrix(exci.right_gs.AR, H_partial, exci.left_gs.AL) - start = scale!(last(found[1:i] * T), cis(-mom * len)) - if exci.trivial && isid(H, i) - @plansor start[-1 -2; -3 -4] -= start[1 4; -3 2] * r_RL(exci.right_gs)[2; 3] * - τ[3 4; 5 1] * l_RL(exci.right_gs)[-1; 6] * τ[5 6; -4 -2] - end - - found[i] = add!(start, lBs[i]) - - if reduce(&, contains.(H.data, i, i)) - if isid(H, i) - tm = TransferMatrix(exci.right_gs.AR, exci.left_gs.AL) - if exci.trivial - tm = regularize(tm, l_RL(exci.right_gs), r_RL(exci.right_gs)) - end - else - tm = TransferMatrix( - exci.right_gs.AR, getindex.(H.data, i, i), exci.left_gs.AL - ) - end - - found[i], convhist = linsolve( - flip(tm), found[i], found[i], solver, 1, -cis(-mom * len) - ) - convhist.converged == 0 && - @warn "GBL$i failed to converge: normres = $(convhist.normres)" - end - end - return found -end function left_excitation_transfer_system( GBL, H::InfiniteMPOHamiltonian, exci; mom = exci.momentum, solver = Defaults.linearsolver @@ -87,49 +43,7 @@ function left_excitation_transfer_system( end return found end -function right_excitation_transfer_system( - rBs, H, exci; mom = exci.momentum, solver = Defaults.linearsolver - ) - len = length(H) - found = zero.(rBs) - odim = length(rBs) - - for i in odim:-1:1 - # this operation can in principle be even further optimized for larger unit cells - # as we only require the terms that end at level i. - # this would require to check the finite state machine, and discard non-connected - # terms. - H_partial = map(site -> H.data[site, i:odim, i:odim], 1:len) - T = TransferMatrix(exci.left_gs.AL, H_partial, exci.right_gs.AR) - start = scale!(first(T * found[i:odim]), cis(mom * len)) - if exci.trivial && isid(H, i) - @plansor start[-1 -2; -3 -4] -= τ[6 2; 3 4] * start[3 4; -3 5] * - l_LR(exci.right_gs)[5; 2] * r_LR(exci.right_gs)[-1; 1] * τ[-2 -4; 1 6] - end - - found[i] = add!(start, rBs[i]) - - if reduce(&, contains.(H.data, i, i)) - if isid(H, i) - tm = TransferMatrix(exci.left_gs.AL, exci.right_gs.AR) - if exci.trivial - tm = regularize(tm, l_LR(exci.left_gs), r_LR(exci.right_gs)) - end - else - tm = TransferMatrix( - exci.left_gs.AL, getindex.(H.data, i, i), exci.right_gs.AR - ) - end - found[i], convhist = linsolve( - tm, found[i], found[i], solver, 1, -cis(mom * len) - ) - convhist.converged < 1 && - @warn "GBR$i failed to converge: normres = $(convhist.normres)" - end - end - return found -end function right_excitation_transfer_system( GBR, H::InfiniteMPOHamiltonian, exci; mom = exci.momentum, From 38d0206d9695f3d7dbc4c291c25abdeb4e6ca80e Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 28 Aug 2025 20:41:08 +0200 Subject: [PATCH 5/9] Add `istrivial` and `istopological` --- .../excitation/exci_transfer_system.jl | 29 +++++++++++-------- .../excitation/quasiparticleexcitation.jl | 8 ++--- src/algorithms/toolbox.jl | 2 +- src/environments/qp_envs.jl | 16 +++++----- src/states/quasiparticle_state.jl | 22 +++++++++----- 5 files changed, 45 insertions(+), 32 deletions(-) diff --git a/src/algorithms/excitation/exci_transfer_system.jl b/src/algorithms/excitation/exci_transfer_system.jl index 8a6e68dd0..c2b4b217d 100644 --- a/src/algorithms/excitation/exci_transfer_system.jl +++ b/src/algorithms/excitation/exci_transfer_system.jl @@ -6,6 +6,11 @@ function left_excitation_transfer_system( found = zerovector(GBL) odim = length(GBL) + if istrivial(exci) + ρ_left = l_RL(exci.right_gs) + ρ_right = r_RL(exci.right_gs) + end + for i in 1:odim # this operation can in principle be even further optimized for larger unit cells # as we only require the terms that end at level i. @@ -14,9 +19,7 @@ function left_excitation_transfer_system( H_partial = map(h -> getindex(h, 1:i, 1, 1, 1:i), parent(H)) T = TransferMatrix(exci.right_gs.AR, H_partial, exci.left_gs.AL) start = scale!(last(found[1:i] * T), cis(-mom * len)) - if exci.trivial && isidentitylevel(H, i) - ρ_left = l_RL(exci.right_gs) - ρ_right = r_RL(exci.right_gs) + if istrivial(exci) && isidentitylevel(H, i) regularize!(start, ρ_right, ρ_left) end @@ -25,8 +28,8 @@ function left_excitation_transfer_system( if !isemptylevel(H, i) if isidentitylevel(H, i) T = TransferMatrix(exci.right_gs.AR, exci.left_gs.AL) - if exci.trivial - T = regularize(T, l_RL(exci.right_gs), r_RL(exci.right_gs)) + if istrivial(exci) + T = regularize(T, ρ_left, ρ_right) end else T = TransferMatrix( @@ -53,6 +56,11 @@ function right_excitation_transfer_system( found = zerovector(GBR) odim = length(GBR) + if istrivial(exci) + ρ_left = l_LR(exci.right_gs) + ρ_right = r_LR(exci.right_gs) + end + for i in odim:-1:1 # this operation can in principle be even further optimized for larger unit cells # as we only require the terms that end at level i. @@ -61,9 +69,7 @@ function right_excitation_transfer_system( H_partial = map(h -> h[i:end, 1, 1, i:end], parent(H)) T = TransferMatrix(exci.left_gs.AL, H_partial, exci.right_gs.AR) start = scale!(first(T * found[i:odim]), cis(mom * len)) - if exci.trivial && isidentitylevel(H, i) - ρ_left = l_LR(exci.right_gs) - ρ_right = r_LR(exci.right_gs) + if istrivial(exci) && isidentitylevel(H, i) regularize!(start, ρ_left, ρ_right) end @@ -72,8 +78,8 @@ function right_excitation_transfer_system( if !isemptylevel(H, i) if isidentitylevel(H, i) tm = TransferMatrix(exci.left_gs.AL, exci.right_gs.AR) - if exci.trivial - tm = regularize(tm, l_LR(exci.left_gs), r_LR(exci.right_gs)) + if istrivial(exci) + tm = regularize(tm, ρ_left, ρ_right) end else tm = TransferMatrix( @@ -82,8 +88,7 @@ function right_excitation_transfer_system( end found[i], convhist = linsolve( - tm, found[i], found[i], solver, 1, - -cis(mom * len) + tm, found[i], found[i], solver, 1, -cis(mom * len) ) convhist.converged < 1 && @warn "GBR$i failed to converge: normres = $(convhist.normres)" diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl index a2fc2759c..0eb400258 100644 --- a/src/algorithms/excitation/quasiparticleexcitation.jl +++ b/src/algorithms/excitation/quasiparticleexcitation.jl @@ -57,7 +57,7 @@ function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renv end function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs; num = 1, kwargs...) # Infer `renvs` in function body as it depends on `solver`. - renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) return excitations(H, alg, ϕ₀, lenvs, renvs; num, kwargs...) end function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP; num = 1, kwargs...) @@ -136,7 +136,7 @@ end function excitations( H, alg::QuasiparticleAnsatz, ϕ₀::FiniteQP, lenvs = environments(ϕ₀.left_gs, H), - renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H); + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H); num = 1 ) E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs) @@ -223,7 +223,7 @@ function excitations( kwargs... ) # Infer `renvs` in function body as it depends on `solver`. - renvs = ϕ₀.trivial ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) + renvs = !istopological(ϕ₀) ? lenvs : environments(ϕ₀.right_gs, H; kwargs...) return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...) end function excitations( @@ -372,7 +372,7 @@ function effective_excitation_renormalization_energy(H, ϕ, lenvs, renvs) E[i] = contract_mpo_expval( ψ_left.AC[i], leftenv(lenvs, i, ψ_left), H[i], rightenv(lenvs, i, ψ_left) ) - if !ϕ.trivial + if istopological(ϕ) E[i] += contract_mpo_expval( ψ_right.AC[i], leftenv(renvs, i, ψ_right), H[i], rightenv(renvs, i, ψ_right) ) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 9bf7c180e..150f276fb 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -225,7 +225,7 @@ end function variance(state::InfiniteQP, H::InfiniteMPOHamiltonian, envs = environments(state, H)) # I remember there being an issue here @gertian? - state.trivial || + istopological(state) && throw(ArgumentError("variance of domain wall excitations is not implemented")) gs = state.left_gs diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl index 483f45252..4aacf0dae 100644 --- a/src/environments/qp_envs.jl +++ b/src/environments/qp_envs.jl @@ -33,7 +33,7 @@ function environments(exci::Union{InfiniteQP, MultilineQP}, H; kwargs...) return environments(exci, H, lenvs; kwargs...) end function environments(exci::Union{InfiniteQP, MultilineQP}, H, lenvs; kwargs...) - renvs = exci.trivial ? lenvs : environments(exci.right_gs, H; kwargs...) + renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H; kwargs...) return environments(exci, H, lenvs, renvs; kwargs...) end @@ -62,7 +62,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; lBs[pos + 1] += leftenv(lenvs, pos, exci.left_gs) * TransferMatrix(exci[pos], H[pos], AL[pos]) / cis(exci.momentum) - if exci.trivial && !isempty(ids) # regularization of trivial excitations + if istrivial(exci) && !isempty(ids) # regularization of trivial excitations ρ_left = l_RL(exci.left_gs, pos + 1) ρ_right = r_RL(exci.left_gs, pos) for i in ids @@ -78,7 +78,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; rBs[pos - 1] += TransferMatrix(exci[pos], H[pos], AR[pos]) * rightenv(renvs, pos, exci.right_gs) * cis(exci.momentum) - if exci.trivial && !isempty(ids) + if istrivial(exci) && !isempty(ids) ρ_left = l_LR(exci.left_gs, pos) ρ_right = r_LR(exci.left_gs, pos - 1) for i in ids @@ -100,7 +100,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; for i in 1:(length(exci) - 1) lB_cur = lB_cur * TransferMatrix(AR[i], H[i], AL[i]) / cis(exci.momentum) - if exci.trivial && !isempty(ids) + if istrivial(exci) && !isempty(ids) ρ_left = l_RL(exci.left_gs, i + 1) ρ_right = r_RL(exci.left_gs, i) for k in ids @@ -115,7 +115,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian, lenvs, renvs; for i in length(exci):-1:2 rB_cur = TransferMatrix(AL[i], H[i], AR[i]) * rB_cur * cis(exci.momentum) - if exci.trivial && !isempty(ids) + if istrivial(exci) && !isempty(ids) ρ_left = l_LR(exci.left_gs, i) ρ_right = r_LR(exci.left_gs, i - 1) for k in ids @@ -132,7 +132,7 @@ end function environments( exci::FiniteQP, H::FiniteMPOHamiltonian, lenvs = environments(exci.left_gs, H), - renvs = exci.trivial ? lenvs : environments(exci.right_gs, H); + renvs = !istopological(exci) ? lenvs : environments(exci.right_gs, H); kwargs... ) AL = exci.left_gs.AL @@ -161,7 +161,7 @@ function environments( end function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...) - exci.trivial || + istopological(exci) && @warn "there is a phase ambiguity in topologically nontrivial statmech excitations" solver = environment_alg(exci, O, exci; kwargs...) @@ -203,7 +203,7 @@ function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...) T_RL = TransferMatrix(right_gs.AR, O, left_gs.AL) T_LR = TransferMatrix(left_gs.AL, O, right_gs.AR) - if exci.trivial + if istrivial(exci) @plansor rvec[-1 -2; -3] := rightenv(lenvs, 0, left_gs)[-1 -2; 1] * conj(left_gs.C[0][-3; 1]) @plansor lvec[-1 -2; -3] := leftenv(lenvs, 1, left_gs)[-1 -2; 1] * diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl index bd1f827b2..cf476f914 100644 --- a/src/states/quasiparticle_state.jl +++ b/src/states/quasiparticle_state.jl @@ -133,7 +133,7 @@ function Base.convert( tm = TransferMatrix(input.left_gs.AL, input.right_gs.AR) - if input.trivial + if istrivial(input) tm = regularize(tm, l_LR(input.right_gs), r_LR(input.right_gs)) end @@ -180,7 +180,7 @@ function Base.convert( end tm = TransferMatrix(input.right_gs.AR, input.left_gs.AL) - if input.trivial + if istrivial(input) tm = regularize(tm, l_RL(input.right_gs), r_RL(input.right_gs)) end @@ -222,9 +222,12 @@ left_virtualspace(state::QP) = map(Base.Fix1(left_virtualspace, state), eachsite right_virtualspace(state::QP, i::Int) = right_virtualspace(state.right_gs, i) right_virtualspace(state::QP) = map(Base.Fix1(right_virtualspace, state), eachsite(state)) auxiliaryspace(state::QP) = space(state.Xs[1], 2) - +auxiliarysector(state::QP) = only(sectors(auxiliaryspace(state))) eachsite(state::QP) = eachsite(state.left_gs) +istopological(qp::QP) = qp.left_gs !== qp.right_gs +istrivial(qp::QP) = !istopological(qp) && isone(auxiliarysector(qp)) + Base.copy(a::QP) = copy!(similar(a), a) Base.copyto!(a::QP, b::QP) = copy!(a, b) function Base.copy!(a::T, b::T) where {T <: QP} @@ -233,11 +236,12 @@ function Base.copy!(a::T, b::T) where {T <: QP} end return a end -function Base.getproperty(v::QP, s::Symbol) +Base.@constprop :aggressive function Base.getproperty(qp::QP, s::Symbol) if s == :trivial - return v.left_gs === v.right_gs + Base.depwarn("`qp.trivial` is deprecated in favor of `istrivial` and `istopological`", :trivial) + return !istopological(qp) else - return getfield(v, s) + return getfield(qp, s) end end @@ -402,7 +406,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S <: FiniteMPS} return FiniteMPS(Ls + Rs + Bs; normalize = false) end -function Base.getproperty(exci::MultilineQP, s::Symbol) +Base.@constprop :aggressive function Base.getproperty(exci::MultilineQP, s::Symbol) if s == :momentum return first(exci.data).momentum elseif s == :left_gs @@ -416,6 +420,10 @@ function Base.getproperty(exci::MultilineQP, s::Symbol) end end +# These should really all be the same, so it might make sense to take the first instead +istrivial(exci::MultilineQP) = all(istrivial, exci.data) +istopological(exci::MultilineQP) = all(istopological, exci.data) + # VectorInterface # --------------- From 1b8d76ccdf97fb243b62cdccf14eb09f5c2b0f6d Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 29 Aug 2025 15:17:56 +0200 Subject: [PATCH 6/9] add Z2 Ising model --- test/setup.jl | 38 +++++++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/test/setup.jl b/test/setup.jl index 09a18a83f..95c544934 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -94,6 +94,14 @@ function S_x(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) w throw(ArgumentError("spin $spin not supported")) end end +function S_x(::Type{Z2Irrep}, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number} + spin == 1 // 2 || throw(ArgumentError("spin $spin not supported")) + pspace = Z2Space(0 => 1, 1 => 1) + X = zeros(T, pspace, pspace) + block(X, Z2Irrep(0)) .= one(T) # no times 2 + block(X, Z2Irrep(1)) .= -one(T) + return X +end function S_y(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number} return if spin == 1 // 2 TensorMap(T[0 -im; im 0], ℂ^2 ← ℂ^2) @@ -121,8 +129,20 @@ end function S_zz(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number} return S_z(Trivial, T; spin) ⊗ S_z(Trivial, T; spin) end +function S_zz(::Type{Z2Irrep}, ::Type{T} = ComplexF64; spin = 1 // 2) where {T <: Number} + spin == 1 // 2 || throw(ArgumentError("spin $spin not supported")) + P = Z2Space(0 => 1, 1 => 1) + ZZ = zeros(ComplexF64, P ⊗ P ← P ⊗ P) + flip_charge(charge::Z2Irrep) = only(charge ⊗ Z2Irrep(1)) + for (s, f) in fusiontrees(ZZ) + if s.uncoupled == map(flip_charge, f.uncoupled) + ZZ[s, f] .= 1 # no divide by 4 + end + end + return ZZ +end -function transverse_field_ising(::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} +function transverse_field_ising(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} X = S_x(Trivial, T; spin = 1 // 2) ZZ = S_zz(Trivial, T; spin = 1 // 2) @@ -138,6 +158,22 @@ function transverse_field_ising(::Type{T} = ComplexF64; g = 1.0, L = Inf) where return H₁ + H₂ end +function transverse_field_ising(::Type{Z2Irrep}, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} + X = S_x(Z2Irrep, T; spin = 1 // 2) + ZZ = S_zz(Z2Irrep, T; spin = 1 // 2) + + if L == Inf + lattice = PeriodicArray([space(X, 1)]) + H₁ = InfiniteMPOHamiltonian(lattice, i => -g * X for i in 1:1) + H₂ = InfiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:1) + else + lattice = fill(space(X, 1), L) + H₁ = FiniteMPOHamiltonian(lattice, i => -g * X for i in 1:L) + H₂ = FiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:(L - 1)) + end + return H₁ + H₂ +end + function heisenberg_XXX(::Type{SU2Irrep}; spin = 1, L = Inf) h = ones(ComplexF64, SU2Space(spin => 1)^2 ← SU2Space(spin => 1)^2) for (c, b) in blocks(h) From b26f3adf5ba51abafdcddfa741dc50d7324fdbd0 Mon Sep 17 00:00:00 2001 From: Boris De Vos Date: Fri, 29 Aug 2025 15:18:27 +0200 Subject: [PATCH 7/9] test sectors in quasiparticle ansatz --- test/algorithms.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/test/algorithms.jl b/test/algorithms.jl index 469506e39..3bd38a755 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -475,6 +475,21 @@ module TestAlgorithms @test energies[1] ≈ 0.41047925 atol = 1.0e-4 @test variance(ϕs[1], H) < 1.0e-8 end + @testset "infinite sector convention" begin + g = 4 + H = repeat(transverse_field_ising(Z2Irrep; g = g, L = Inf), 2) + V = Z2Space(0 => 24, 1 => 24) + P = Z2Space(0 => 1, 1 => 1) + ψ = InfiniteMPS([P, P], [V, V]) + ψ, envs = find_groundstate(ψ, H, VUMPS(; tol = 1.0e-10, maxiter = 400)) + + # testing istrivial and istopological + momentum = 0 + exc0, qp0 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(0)) + exc1, qp1 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(1)) + @test isapprox(first(exc1), abs(2 * (g - 1)); atol = 1.0e-6) # charged excitation lower in energy + @test variance(qp1[1], H) < 1.0e-8 + end @testset "infinite (mpo)" begin H = repeat(sixvertex(), 2) ψ = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^10, ℂ^10]) From f4de7362d737e433e762431a234f6f69fd7c30a9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 31 Aug 2025 13:04:31 +0200 Subject: [PATCH 8/9] Merge transverse field ising implementations --- test/setup.jl | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/test/setup.jl b/test/setup.jl index 95c544934..a6b07bbc0 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -142,25 +142,9 @@ function S_zz(::Type{Z2Irrep}, ::Type{T} = ComplexF64; spin = 1 // 2) where {T < return ZZ end -function transverse_field_ising(::Type{Trivial} = Trivial, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} - X = S_x(Trivial, T; spin = 1 // 2) - ZZ = S_zz(Trivial, T; spin = 1 // 2) - - if L == Inf - lattice = PeriodicArray([ℂ^2]) - H₁ = InfiniteMPOHamiltonian(lattice, i => -g * X for i in 1:1) - H₂ = InfiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:1) - else - lattice = fill(ℂ^2, L) - H₁ = FiniteMPOHamiltonian(lattice, i => -g * X for i in 1:L) - H₂ = FiniteMPOHamiltonian(lattice, (i, i + 1) => -ZZ for i in 1:(L - 1)) - end - return H₁ + H₂ -end - -function transverse_field_ising(::Type{Z2Irrep}, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} - X = S_x(Z2Irrep, T; spin = 1 // 2) - ZZ = S_zz(Z2Irrep, T; spin = 1 // 2) +function transverse_field_ising(sym::Type{<:Sector} = Trivial, ::Type{T} = ComplexF64; g = 1.0, L = Inf) where {T <: Number} + X = S_x(sym, T; spin = 1 // 2) + ZZ = S_zz(sym, T; spin = 1 // 2) if L == Inf lattice = PeriodicArray([space(X, 1)]) From de2b0217fc8cb33743ed9763b4da44d8b40a9e42 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Sun, 31 Aug 2025 13:14:59 +0200 Subject: [PATCH 9/9] small tests alterations --- test/algorithms.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/algorithms.jl b/test/algorithms.jl index 3bd38a755..b009c121b 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -479,8 +479,7 @@ module TestAlgorithms g = 4 H = repeat(transverse_field_ising(Z2Irrep; g = g, L = Inf), 2) V = Z2Space(0 => 24, 1 => 24) - P = Z2Space(0 => 1, 1 => 1) - ψ = InfiniteMPS([P, P], [V, V]) + ψ = InfiniteMPS(physicalspace(H), [V, V]) ψ, envs = find_groundstate(ψ, H, VUMPS(; tol = 1.0e-10, maxiter = 400)) # testing istrivial and istopological @@ -488,6 +487,7 @@ module TestAlgorithms exc0, qp0 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(0)) exc1, qp1 = excitations(H, QuasiparticleAnsatz(), momentum, ψ; sector = Z2Irrep(1)) @test isapprox(first(exc1), abs(2 * (g - 1)); atol = 1.0e-6) # charged excitation lower in energy + @test first(exc0) > 2 * first(exc1) @test variance(qp1[1], H) < 1.0e-8 end @testset "infinite (mpo)" begin