Refactor quasiparticle regularization (#318)

* Fix off-by-one error

* centralize regularization calls

* avoid braiding in regularize

* remove unused code

* Add `istrivial` and `istopological`

* add Z2 Ising model

* test sectors in quasiparticle ansatz

* Merge transverse field ising implementations

* small tests alterations

---------

Co-authored-by: Boris De Vos <[email protected]>

Author: lkdvos

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
@@ -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.
@@ -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(
@@ -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,
@@ -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.
@@ -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(
@@ -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)"

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

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
 

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
 
@@ -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
@@ -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
@@ -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
@@ -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
 
@@ -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
 
@@ -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
@@ -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...)
 
@@ -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] *