Skip to content

Refactor quasiparticle regularization #318

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
lkdvos merged 9 commits into from
Sep 1, 2025
Merged

Refactor quasiparticle regularization #318

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0f50dc5
Fix off-by-one error
lkdvos Aug 28, 2025
3fe94d2
centralize regularization calls
lkdvos Aug 28, 2025
fe444ba
avoid braiding in regularize
lkdvos Aug 28, 2025
f99504a
remove unused code
lkdvos Aug 28, 2025
38d0206
Add `istrivial` and `istopological`
lkdvos Aug 28, 2025
1b8d76c
add Z2 Ising model
borisdevos Aug 29, 2025
b26f3ad
test sectors in quasiparticle ansatz
borisdevos Aug 29, 2025
f4de736
Merge transverse field ising implementations
lkdvos Aug 31, 2025
de2b021
small tests alterations
lkdvos Aug 31, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
117 changes: 19 additions & 98 deletions src/algorithms/excitation/exci_transfer_system.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -50,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.
Expand All @@ -58,18 +19,17 @@ 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)
@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]
if istrivial(exci) && isidentitylevel(H, i)
regularize!(start, ρ_right, ρ_left)
end

found[i] = add!(start, GBL[i])

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(
Expand All @@ -86,49 +46,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,
Expand All @@ -138,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.
Expand All @@ -146,18 +69,17 @@ 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)
@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]
if istrivial(exci) && isidentitylevel(H, i)
regularize!(start, ρ_left, ρ_right)
end

found[i] = add!(start, GBR[i])

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(
Expand All @@ -166,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)"
Expand Down
8 changes: 4 additions & 4 deletions src/algorithms/excitation/quasiparticleexcitation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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...)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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)
)
Expand Down
2 changes: 1 addition & 1 deletion src/algorithms/toolbox.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
45 changes: 21 additions & 24 deletions src/environments/qp_envs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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 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
@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
Expand All @@ -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 istrivial(exci) && !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
Expand All @@ -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 istrivial(exci) && !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

Expand All @@ -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 istrivial(exci) && !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

Expand All @@ -135,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
Expand Down Expand Up @@ -164,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...)

Expand Down Expand Up @@ -206,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] *
Expand Down
Loading
Loading