Refactor InfiniteQPEnvironments

lkdvos · lkdvos · commit 0a0c3f93e088 · 2025-01-14T15:41:43.000-05:00
diff --git a/src/algorithms/excitation/quasiparticleexcitation.jl b/src/algorithms/excitation/quasiparticleexcitation.jl
@@ -168,7 +168,7 @@ end
 #                           Statmech Excitations                               #
 ################################################################################
 
-function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::Multiline{<:InfiniteQP},
+function excitations(H::MultilineMPO, alg::QuasiparticleAnsatz, ϕ₀::MultilineQP,
                      lenvs, renvs; num=1, kwargs...)
     qp_envs(ϕ) = environments(ϕ, H, lenvs, renvs; kwargs...)
     function H_eff(ϕ′)
@@ -214,13 +214,15 @@ function excitations(H::DenseMPO, alg::QuasiparticleAnsatz, momentum::Real,
                      renvs=lmps === rmps ? lenvs : environments(rmps, H);
                      sector=one(sectortype(lmps)), num=1, kwargs...)
     multiline_lmps = convert(MultilineMPS, lmps)
+    lenvs′ = Multiline([lenvs])
     if lmps === rmps
-        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs,
+        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs′,
                     multiline_lmps,
-                    lenvs; sector, num, kwargs...)
+                    lenvs′; sector, num, kwargs...)
     else
-        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs,
-                    convert(MultilineMPS, rmps), renvs; sector, num, kwargs...)
+        renvs′ = Multiline([renvs])
+        excitations(convert(MultilineMPO, H), alg, momentum, multiline_lmps, lenvs′,
+                    convert(MultilineMPS, rmps), renvs′; sector, num, kwargs...)
     end
 end
 
@@ -254,41 +256,10 @@ function effective_excitation_hamiltonian(H::Union{InfiniteMPOHamiltonian,
     return ϕ′
 end
 
-function effective_excitation_hamiltonian(H::MultilineMPO, ϕ::Multiline{<:InfiniteQP},
+function effective_excitation_hamiltonian(H::MultilineMPO, ϕ::MultilineQP,
                                           envs=environments(ϕ, H))
-    ϕ′ = Multiline(similar.(ϕ.data))
-    left_gs = ϕ.left_gs
-    right_gs = ϕ.right_gs
-
-    for row in 1:size(H, 1)
-        Bs = [ϕ[row][i] for i in 1:size(H, 2)]
-        for col in 1:size(H, 2)
-            en = @plansor conj(left_gs.AC[row, col][2 6; 4]) *
-                          leftenv(envs.lenvs, row, col, left_gs)[2 5; 3] *
-                          left_gs.AC[row + 1, col][3 7; 1] *
-                          H[row, col][5 6; 7 8] *
-                          rightenv(envs.lenvs, row, col, left_gs)[1 8; 4]
-
-            @plansor T[-1 -2; -3 -4] := leftenv(envs.lenvs, row, col, left_gs)[-1 5; 4] *
-                                        Bs[col][4 2; -3 1] *
-                                        H[row, col][5 -2; 2 3] *
-                                        rightenv(envs.renvs, row, col, right_gs)[1 3; -4]
-
-            @plansor T[-1 -2; -3 -4] += envs.lBs[row, col - 1][-1 4; -3 5] *
-                                        right_gs.AR[row, col][5 2; 1] *
-                                        H[row, col][4 -2; 2 3] *
-                                        rightenv(envs.renvs, row, col, right_gs)[1 3; -4]
-
-            @plansor T[-1 -2; -3 -4] += leftenv(envs.lenvs, row, col, left_gs)[-1 2; 1] *
-                                        left_gs.AL[row, col][1 3; 4] *
-                                        H[row, col][2 -2; 3 5] *
-                                        envs.rBs[row, col + 1][4 5; -3 -4]
-
-            ϕ′[row + 1][col] = T / en
-        end
-    end
-
-    return ϕ′
+    return Multiline(map(effective_excitation_hamiltonian,
+                         parent(H), parent(ϕ), parent(envs)))
 end
 
 function effective_excitation_hamiltonian(H::InfiniteMPO, ϕ::InfiniteQP,
diff --git a/src/environments/qp_envs.jl b/src/environments/qp_envs.jl
@@ -1,28 +1,27 @@
-#=
-nothing fancy - only used internally (and therefore cryptic) - stores some partially contracted things
-seperates out this bit of logic from effective_excitation_hamiltonian (now more readable)
-can also - potentially - partially reuse this in other algorithms
-=#
-struct QPEnv{A,B} <: AbstractMPSEnvironments
-    lBs::PeriodicArray{A,2}
-    rBs::PeriodicArray{A,2}
-
-    lenvs::B
-    renvs::B
-end
-
-struct QuasiparticleEnvironments{A,B} <: AbstractMPSEnvironments
+"""
+    InfiniteQPEnvironments <: AbstractMPSEnvironments
+
+Environments for an infinite QP-MPO-QP combination. These solve the corresponding fixedpoint equations:
+```math
+GLs[i] * T_BL[i] + GBLs[i] * T_RL[i] = GBLs[i + 1]
+T_BR[i] * GRs[i] + T_LR[i] * GBRs[i] = GBRs[i - 1]
+```
+where `T_BL`, `T_BR`, `T_RL` and `T_LR` are the (regularized) transfer matrix operators on a given site for `B-O-AL`, `B-O-AR`, `AR-O-AL` and `AL-O-AR` respectively.
+"""
+struct InfiniteQPEnvironments{A,B} <: AbstractMPSEnvironments
     leftBenvs::PeriodicVector{A}
     rightBenvs::PeriodicVector{A}
 
     leftenvs::B
     rightenvs::B
 end
 
-function leftenv(envs::QuasiparticleEnvironments, site::Int, state)
+Base.length(envs::InfiniteQPEnvironments) = length(envs.leftenvs)
+
+function leftenv(envs::InfiniteQPEnvironments, site::Int, state)
     return leftenv(envs.leftenvs, site, state)
 end
-function rightenv(envs::QuasiparticleEnvironments, site::Int, state)
+function rightenv(envs::InfiniteQPEnvironments, site::Int, state)
     return rightenv(envs.rightenvs, site, state)
 end
 
@@ -39,20 +38,26 @@ function environments(exci::Union{InfiniteQP,MultilineQP}, H, lenvs;
     return environments(exci, H, lenvs, renvs; kwargs...)
 end
 
+function environments(qp::MultilineQP, operator::MultilineMPO, lenvs, renvs; kwargs...)
+    (rows = size(qp, 1)) == size(operator, 1) || throw(ArgumentError("Incompatible sizes"))
+    envs = map(1:rows) do row
+        return environments(qp[row], operator[row], lenvs[row], renvs[row]; kwargs...)
+    end
+    return Multiline(PeriodicVector(envs))
+end
+
 function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian,
                       lenvs, renvs;
                       kwargs...)
     ids = findall(Base.Fix1(isidentitylevel, H), 2:(size(H[1], 1) - 1))
     solver = environment_alg(exci, H, exci; kwargs...)
-    # ids = collect(Iterators.filter(x -> isidentitylevel(H, x), 2:(H.odim - 1)))
 
     AL = exci.left_gs.AL
     AR = exci.right_gs.AR
 
     lBs = PeriodicVector([allocate_GBL(exci, H, exci, i) for i in 1:length(exci)])
     rBs = PeriodicVector([allocate_GBR(exci, H, exci, i) for i in 1:length(exci)])
 
-    # lBs, rBs = gen_exci_lw_rw(exci.left_gs, H, exci.right_gs, space(exci[1], 3))
     zerovector!(lBs[1])
     for pos in 1:length(exci)
         lBs[pos + 1] = lBs[pos] * TransferMatrix(AR[pos], H[pos], AL[pos]) /
@@ -139,7 +144,7 @@ function environments(exci::InfiniteQP, H::InfiniteMPOHamiltonian,
         rBs[i - 1] += rB_cur
     end
 
-    return QuasiparticleEnvironments(lBs, rBs, lenvs, renvs)
+    return InfiniteQPEnvironments(lBs, rBs, lenvs, renvs)
 end
 
 function environments(exci::FiniteQP,
@@ -168,149 +173,10 @@ function environments(exci::FiniteQP,
                         rightenv(renvs, pos, exci.right_gs)
     end
 
-    return QuasiparticleEnvironments(lBs, rBs, lenvs, renvs)
+    return InfiniteQPEnvironments(lBs, rBs, lenvs, renvs)
 end
 
-function environments(exci::Multiline{<:InfiniteQP},
-                      ham::MultilineMPO,
-                      lenvs,
-                      renvs;
-                      kwargs...)
-    exci.trivial ||
-        @warn "there is a phase ambiguity in topologically nontrivial statmech excitations"
-    solver = environment_alg(exci, ham, exci; kwargs...)
-
-    left_gs = exci.left_gs
-    right_gs = exci.right_gs
-
-    exci_space = space(exci[1][1], 3)
-
-    (numrows, numcols) = size(left_gs)
-
-    st = site_type(typeof(left_gs))
-    B_type = tensormaptype(spacetype(st), 2, 2, storagetype(st))
-
-    lBs = PeriodicArray{B_type,2}(undef, size(left_gs, 1), size(left_gs, 2))
-    rBs = PeriodicArray{B_type,2}(undef, size(left_gs, 1), size(left_gs, 2))
-
-    for row in 1:numrows
-        c_lenvs = broadcast(col -> leftenv(lenvs, col, left_gs)[row], 1:numcols)
-        c_renvs = broadcast(col -> rightenv(renvs, col, right_gs)[row], 1:numcols)
-
-        hamrow = ham[row, :]
-
-        left_above = left_gs[row]
-        left_below = left_gs[row + 1]
-        right_above = right_gs[row]
-        right_below = right_gs[row + 1]
-
-        left_renorms = fill(zero(scalartype(B_type)), numcols)
-        right_renorms = fill(zero(scalartype(B_type)), numcols)
-
-        for col in 1:numcols
-            lv = leftenv(lenvs, col, left_gs)[row]
-            rv = rightenv(lenvs, col, left_gs)[row]
-            left_renorms[col] = @plansor lv[1 2; 3] *
-                                         left_above.AC[col][3 4; 5] *
-                                         hamrow[col][2 6; 4 7] *
-                                         rv[5 7; 8] *
-                                         conj(left_below.AC[col][1 6; 8])
-
-            lv = leftenv(renvs, col, right_gs)[row]
-            rv = rightenv(renvs, col, right_gs)[row]
-            right_renorms[col] = @plansor lv[1 2; 3] *
-                                          right_above.AC[col][3 4; 5] *
-                                          hamrow[col][2 6; 4 7] *
-                                          rv[5 7; 8] *
-                                          conj(right_below.AC[col][1 6; 8])
-        end
-
-        left_renorms = left_renorms .^ -1
-        right_renorms = right_renorms .^ -1
-
-        lB_cur = fill_data!(similar(left_below.AL[1],
-                                    left_virtualspace(left_below, 0) *
-                                    _firstspace(hamrow[1])',
-                                    exci_space' * right_virtualspace(right_above, 0)),
-                            zero)
-        rB_cur = fill_data!(similar(left_below.AL[1],
-                                    left_virtualspace(left_below, 0) *
-                                    _firstspace(hamrow[1]),
-                                    exci_space' * right_virtualspace(right_above, 0)),
-                            zero)
-        for col in 1:numcols
-            lB_cur = lB_cur *
-                     TransferMatrix(right_above.AR[col], hamrow[col], left_below.AL[col])
-            lB_cur += c_lenvs[col] *
-                      TransferMatrix(exci[row][col], hamrow[col], left_below.AL[col])
-            lB_cur *= left_renorms[col] * exp(-1im * exci.momentum)
-            lBs[row, col] = lB_cur
-
-            col = numcols - col + 1
-
-            rB_cur = TransferMatrix(left_above.AL[col], hamrow[col], right_below.AR[col]) *
-                     rB_cur
-            rB_cur += TransferMatrix(exci[row][col], hamrow[col], right_below.AR[col]) *
-                      c_renvs[col]
-            rB_cur *= exp(1im * exci.momentum) * right_renorms[col]
-            rBs[row, col] = rB_cur
-        end
-
-        tm_RL = TransferMatrix(right_above.AR, hamrow, left_below.AL)
-        tm_LR = TransferMatrix(left_above.AL, hamrow, right_below.AR)
-
-        if exci.trivial
-            @plansor rvec[-1 -2; -3] := rightenv(lenvs, 0, left_gs)[row][-1 -2; 1] *
-                                        conj(left_below.C[0][-3; 1])
-            @plansor lvec[-1 -2; -3] := leftenv(lenvs, 1, left_gs)[row][-1 -2; 1] *
-                                        left_above.C[0][1; -3]
-
-            tm_RL = regularize(tm_RL, lvec, rvec)
-
-            @plansor rvec[-1 -2; -3] := rightenv(renvs, 0, right_gs)[row][1 -2; -3] *
-                                        right_above.C[0][-1; 1]
-            @plansor lvec[-1 -2; -3] := conj(right_below.C[0][-3; 1]) *
-                                        leftenv(renvs, 1, right_gs)[row][-1 -2; 1]
-
-            tm_LR = regularize(tm_LR, lvec, rvec)
-        end
-
-        lBs[row, end], convhist = linsolve(flip(tm_RL), lB_cur, lB_cur, solver, 1,
-                                           -exp(-1im * numcols * exci.momentum) *
-                                           prod(left_renorms))
-        convhist.converged == 0 &&
-            @warn "GBL[$row] failed to converge: normres = $(convhist.normres)"
-
-        rBs[row, 1], convhist = linsolve(tm_LR, rB_cur, rB_cur, GMRES(), 1,
-                                         -exp(1im * numcols * exci.momentum) *
-                                         prod(right_renorms))
-        convhist.converged == 0 &&
-            @warn "GBR[$row] failed to converge: normres = $(convhist.normres)"
-
-        left_cur = lBs[row, end]
-        right_cur = rBs[row, 1]
-        for col in 1:(numcols - 1)
-            left_cur = left_renorms[col] * left_cur *
-                       TransferMatrix(right_above.AR[col], hamrow[col],
-                                      left_below.AL[col]) * exp(-1im * exci.momentum)
-            lBs[row, col] += left_cur
-
-            col = numcols - col + 1
-            right_cur = TransferMatrix(left_above.AL[col], hamrow[col],
-                                       right_below.AR[col]) * right_cur *
-                        exp(1im * exci.momentum) * right_renorms[col]
-            rBs[row, col] += right_cur
-        end
-    end
-
-    return QPEnv(lBs, rBs, lenvs, renvs)
-end
-
-function environments(exci::InfiniteQP,
-                      O::InfiniteMPO,
-                      lenvs,
-                      renvs;
-                      kwargs...)
+function environments(exci::InfiniteQP, O::InfiniteMPO, lenvs, renvs; kwargs...)
     exci.trivial ||
         @warn "there is a phase ambiguity in topologically nontrivial statmech excitations"
     solver = environment_alg(exci, O, exci; kwargs...)
@@ -397,5 +263,5 @@ function environments(exci::InfiniteQP,
         GBR[col] += right_cur
     end
 
-    return QuasiparticleEnvironments(GBL, GBR, lenvs, renvs)
+    return InfiniteQPEnvironments(GBL, GBR, lenvs, renvs)
 end
diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl
@@ -351,7 +351,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S<:FiniteMPS}
     return FiniteMPS(Ls + Rs + Bs; normalize=false)
 end
 
-function Base.getproperty(exci::Multiline{<:InfiniteQP}, s::Symbol)
+function Base.getproperty(exci::MultilineQP, s::Symbol)
     if s == :momentum
         return first(exci.data).momentum
     elseif s == :left_gs
