Compatibility with multifusion categories

src/algorithms/ED.jl

@@ -44,7 +44,7 @@ function exact_diagonalization(
 
     # fuse from left to right
     ALs = Vector{Union{Missing, TA}}(missing, L)
-    left = oneunit(spacetype(H))
+    left = spacetype(H)(one(sector) => 1) # might need to be rightone, leave this for now
     for i in 1:(middle_site - 1)
         P = physicalspace(H, i)
         ALs[i] = isomorphism(T, left ⊗ P ← fuse(left ⊗ P))

src/algorithms/correlators.jl

@@ -39,8 +39,9 @@ function correlator(
 end
 
 function correlator(
-        state::AbstractMPS, O₁₂::AbstractTensorMap{<:Any, S, 2, 2}, i::Int, j
+        state::AbstractMPS, O₁₂::AbstractTensorMap{<:Any, S, 2, 2}, i::Int,
+        j
     ) where {S}
-    O₁, O₂ = decompose_localmpo(add_util_leg(O₁₂))
+    O₁, O₂ = decompose_localmpo(add_util_mpoleg(O₁₂))
     return correlator(state, O₁, O₂, i, j)
 end

src/algorithms/excitation/exci_transfer_system.jl

@@ -59,8 +59,14 @@ 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]
+            # not using braiding tensors here, leads to extra leg
+            util = similar(exci.left_gs.AL[1], space(parent(H)[1], 1)[1])
+            fill_data!(util, one)
+            @plansor start[-1 -2; -3 -4] -= start[2 1; -3 3] *
+                util[1] *
+                r_RL(exci.right_gs)[3; 2] *
+                l_RL(exci.right_gs)[-1; -4] *
+                conj(util[-2])
         end
 
         found[i] = add!(start, GBL[i])
@@ -69,7 +75,10 @@ function left_excitation_transfer_system(
             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))
+                    # deal with extra leg
+                    @plansor lRL_util[-1 -2; -3] := l_RL(exci.right_gs)[-1;-3] * conj(util[-2])
+                    @plansor rRL_util[-1 -2; -3] := r_RL(exci.right_gs)[-1;-3] * util[-2]
+                    T = regularize(T, lRL_util, rRL_util)
                 end
             else
                 T = TransferMatrix(
@@ -129,6 +138,7 @@ function right_excitation_transfer_system(
     end
     return found
 end
+#TODO: check that the left/right discrepancies make sense
 function right_excitation_transfer_system(
         GBR, H::InfiniteMPOHamiltonian, exci;
         mom = exci.momentum,
@@ -147,8 +157,14 @@ 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]
+            # not using braiding tensors here, leads to extra leg
+            util = similar(exci.right_gs.AL[1], space(parent(H)[1], 1)[1])
+            fill_data!(util, one)
+            @plansor start[-1 -2; -3 -4] -= start[2 1; -3 3] *
+                conj(util[1]) *
+                l_LR(exci.right_gs)[3; 2] *
+                r_LR(exci.right_gs)[-1; -4] *
+                util[-2]
         end
 
         found[i] = add!(start, GBR[i])
@@ -157,7 +173,10 @@ function right_excitation_transfer_system(
             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))
+                    # deal with extra leg
+                    @plansor lLR_util[-1 -2; -3] := l_LR(exci.left_gs)[-1;-3] * conj(util[-2])
+                    @plansor rLR_util[-1 -2; -3] := r_LR(exci.right_gs)[-1;-3] * util[-2]
+                    tm = regularize(tm, lLR_util, rLR_util)
                 end
             else
                 tm = TransferMatrix(

src/algorithms/excitation/quasiparticleexcitation.jl

@@ -46,7 +46,6 @@ end
 function excitations(H, alg::QuasiparticleAnsatz, ϕ₀::InfiniteQP, lenvs, renvs; num::Int = 1)
     E = effective_excitation_renormalization_energy(H, ϕ₀, lenvs, renvs)
     H_eff = EffectiveExcitationHamiltonian(H, lenvs, renvs, E)
-
     Es, ϕs, convhist = eigsolve(ϕ₀, num, :SR, alg.alg) do ϕ
         return H_eff(ϕ; alg.alg_environments...)
     end
@@ -105,13 +104,25 @@ function excitations(
         verbosity = Defaults.verbosity, num = 1,
         sector = one(sectortype(lmps)), parallel = true, kwargs...
     )
-    Toutput = Core.Compiler.return_type(
-        excitations,
-        Tuple{
-            typeof(H), typeof(alg), eltype(momenta), typeof(lmps),
-            typeof(lenvs), typeof(rmps), typeof(renvs),
-        }
-    )
+    # wrapper to evaluate sector as positional argument
+    Toutput = let
+        function wrapper(H, alg, p, lmps, lenvs, rmps, renvs, sector; kwargs...)
+            return excitations(
+                H, alg, p, lmps, lenvs, rmps, renvs;
+                sector, kwargs...
+            )
+        end
+        Core.Compiler.return_type(
+            wrapper,
+            Tuple{
+                typeof(H), typeof(alg),
+                eltype(momenta), typeof(lmps),
+                typeof(lenvs),
+                typeof(rmps), typeof(renvs), typeof(sector),
+            }
+        )
+    end
+
     results = similar(momenta, Toutput)
     scheduler = parallel ? :greedy : :serial
     tmap!(results, momenta; scheduler) do momentum

src/algorithms/fixedpoint.jl

@@ -4,14 +4,14 @@
     fixedpoint(A, x₀, which::Symbol; kwargs...) -> val, vec
     fixedpoint(A, x₀, which::Symbol, alg) -> val, vec
 
-Compute the fixedpoint of a linear operator `A` using the specified eigensolver `alg`. The
+Compute the fixed point of a linear operator `A` using the specified eigensolver `alg`. The
 fixedpoint is assumed to be unique.
 """
 function fixedpoint(A, x₀, which::Symbol, alg::Lanczos)
     vals, vecs, info = eigsolve(A, x₀, 1, which, alg)
 
     if info.converged == 0
-        @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
+        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     end
 
     return vals[1], vecs[1]
@@ -21,10 +21,10 @@ function fixedpoint(A, x₀, which::Symbol, alg::Arnoldi)
     TT, vecs, vals, info = schursolve(A, x₀, 1, which, alg)
 
     if info.converged == 0
-        @warnv 1 "fixedpoint not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
+        @warnv 1 "fixed point not converged after $(info.numiter) iterations: normres = $(info.normres[1])"
     end
     if size(TT, 2) > 1 && TT[2, 1] != 0
-        @warnv 1 "non-unique fixedpoint detected"
+        @warnv 1 "non-unique fixed point detected"
     end
 
     return vals[1], vecs[1]

src/algorithms/groundstate/gradient_grassmann.jl

@@ -59,7 +59,7 @@ function find_groundstate(
         ψ::S, H, alg::GradientGrassmann, envs::P = environments(ψ, H)
     )::Tuple{S, P, Float64} where {S, P}
     !isa(ψ, FiniteMPS) || dim(ψ.C[end]) == 1 ||
-        @warn "This is not fully supported - split the mps up in a sum of mps's and optimize seperately"
+        @warn "This is not fully supported - split the mps up in a sum of mps's and optimize separately"
     normalize!(ψ)
 
     fg(x) = GrassmannMPS.fg(x, H, envs)

src/algorithms/toolbox.jl

@@ -28,7 +28,7 @@ Return the density matrix of the infinite temperature state for a given Hamilton
 This is the identity matrix in the physical space, and the identity in the auxiliary space.
 """
 function infinite_temperature_density_matrix(H::MPOHamiltonian)
-    V = oneunit(spacetype(H))
+    V = oneunit(space(first(H[1]), 1))
     W = map(1:length(H)) do site
         return BraidingTensor{scalartype(H)}(physicalspace(H, site), V)
     end

src/operators/mpo.jl

@@ -27,7 +27,7 @@ function FiniteMPO(Os::AbstractVector{O}) where {O}
 end
 
 function FiniteMPO(O::AbstractTensorMap{T, S, N, N}) where {T, S, N}
-    return FiniteMPO(decompose_localmpo(add_util_leg(O)))
+    return FiniteMPO(decompose_localmpo(add_util_mpoleg(O)))
 end
 
 """

src/operators/mpohamiltonian.jl

@@ -384,13 +384,19 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, local_opera
     E = scalartype(T)
     S = spacetype(T)
 
+    # avoid using one(S)
+    somempo = local_mpos[1].second[1]
+    sp_oneleg = space(somempo, 1)
+
+    _oneunit = oneunit(sp_oneleg) # should be rightoneunit, but MPOHamiltonians are always diagonal for now
+
     virtualsumspaces = Vector{SumSpace{S}}(undef, length(lattice) + 1)
-    virtualsumspaces[1] = SumSpace(fill(oneunit(S), 1))
-    virtualsumspaces[end] = SumSpace(fill(oneunit(S), 1))
+    virtualsumspaces[1] = SumSpace(fill(_oneunit, 1))
+    virtualsumspaces[end] = SumSpace(fill(_oneunit, 1))
 
     for i in 1:(length(lattice) - 1)
         n_channels = maximum(last, nonzero_keys[i]; init = 1) + 1
-        V = SumSpace(fill(oneunit(S), n_channels))
+        V = SumSpace(fill(_oneunit, n_channels))
         if n_channels > 2
             for ((key_L, key_R), O) in zip(nonzero_keys[i], nonzero_opps[i])
                 V[key_R == 0 ? end : key_R] = if O isa Number
@@ -468,9 +474,14 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, local_
     virtualspaces = PeriodicArray(
         [Vector{MissingS}(missing, operator_size) for _ in 1:length(nonzero_keys)]
     )
+    # avoid using one(S)
+    somempo = local_mpos[1].second[1]
+    sp_oneleg = space(somempo, 1)
+
+    _oneunit = oneunit(sp_oneleg)
     for V in virtualspaces
-        V[1] = oneunit(S)
-        V[end] = oneunit(S)
+        V[1] = _oneunit
+        V[end] = _oneunit
     end
 
     # start by filling in tensors -> space information available
@@ -516,7 +527,7 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, local_
         end
     end
 
-    foreach(Base.Fix2(replace!, missing => oneunit(S)), virtualspaces)
+    foreach(Base.Fix2(replace!, missing => _oneunit), virtualspaces)
     virtualsumspaces = map(virtualspaces) do V
         return SumSpace(collect(S, V))
     end
@@ -658,7 +669,7 @@ function Base.:+(
     ) where {O <: JordanMPOTensor}
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
-    Vtriv = oneunit(spacetype(H₁))
+    Vtriv = oneunit(space(first(H₁[1]), 1)) # should also be rightoneunit, but is currently diagonal
 
     for i in 1:N
         A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
@@ -682,7 +693,7 @@ function Base.:+(
     ) where {O <: JordanMPOTensor}
     N = check_length(H₁, H₂)
     H = similar(parent(H₁))
-    Vtriv = oneunit(spacetype(H₁))
+    Vtriv = oneunit(space(first(H₁[1]), 1))
     for i in 1:N
         A = cat(H₁[i].A, H₂[i].A; dims = (1, 4))
         B = cat(H₁[i].B, H₂[i].B; dims = 1)

src/utility/utility.jl

@@ -69,6 +69,19 @@ function add_util_leg(tensor::AbstractTensorMap{T, S, N1, N2}) where {T, S, N1,
     return util_front * tensor * util_back
 end
 
+function add_util_mpoleg(tensor::AbstractTensorMap{T, S, N1, N2}) where {T, S, N1, N2}
+    # separate function for mpo because add_util_leg is also used for mps from tensors
+    # and the additional legs there can be different depending on the fusion tree
+
+    #TODO?: might need a left and right utility space if MPO has module virtual space (doesn't happen for hamiltonians)
+    ou = S(one(first(blocksectors(tensor))) => 1)
+
+    util_front = isomorphism(storagetype(tensor), ou * codomain(tensor), codomain(tensor))
+    util_back = isomorphism(storagetype(tensor), domain(tensor), domain(tensor) * ou)
+
+    return util_front * tensor * util_back
+end
+
 function union_split(a::AbstractArray)
     T = reduce((a, b) -> Union{a, b}, typeof.(a))
     nA = similar(a, T)

test/algorithms.jl

@@ -710,7 +710,7 @@ module TestAlgorithms
 
         @test expectation_value(ψ, H) ≈
             expectation_value(ψ, 1 => -g * S_x()) + expectation_value(ψ, (1, 2) => -S_zz())
-        Z_mpo = MPSKit.add_util_leg(S_z())
+        Z_mpo = MPSKit.add_util_mpoleg(S_z())
         G = correlator(ψ, Z_mpo, Z_mpo, 1, 2:5)
         G2 = correlator(ψ, S_zz(), 1, 3:2:5)
         @test isapprox(G[2], G2[1], atol = 1.0e-2)