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 4 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
96 changes: 6 additions & 90 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 @@ -59,8 +15,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])
Expand All @@ -86,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,
Expand All @@ -147,8 +62,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])
Expand Down
37 changes: 17 additions & 20 deletions src/environments/qp_envs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 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
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 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
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 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

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 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

Expand Down
4 changes: 2 additions & 2 deletions src/transfermatrix/transfermatrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading