Compatibility with multifusion categories

.github/workflows/Tests.yml

@@ -29,6 +29,7 @@ jobs:
           - states
           - operators
           - algorithms
+          - multifusion
           - other
         os:
           - ubuntu-latest

src/algorithms/ED.jl

@@ -1,6 +1,6 @@
 """
     exact_diagonalization(H::FiniteMPOHamiltonian;
-                          sector=first(sectors(oneunit(physicalspace(H, 1)))),
+                          sector=rightunit(H),
                           len::Int=length(H), num::Int=1, which::Symbol=:SR,
                           alg=Defaults.alg_eigsolve(; dynamic_tols=false))
                             -> vals, state_vecs, convhist
@@ -13,7 +13,7 @@ equivalent to dense eigenvectors.
 - `H::FiniteMPOHamiltonian`: the Hamiltonian to diagonalize.
 
 ### Keyword arguments
-- `sector=first(sectors(oneunit(physicalspace(H, 1))))`: the total charge of the
+- `sector=rightunit(H)`: the total charge of the
   eigenvectors, which is chosen trivial by default.
 - `len::Int=length(H)`: the length of the system.
 - `num::Int=1`: the number of eigenvectors to find.
@@ -32,7 +32,7 @@ equivalent to dense eigenvectors.
 """
 function exact_diagonalization(
         H::FiniteMPOHamiltonian;
-        sector = one(sectortype(H)), num::Int = 1, which::Symbol = :SR,
+        sector = rightunit(H), num::Int = 1, which::Symbol = :SR,
         alg = Defaults.alg_eigsolve(; dynamic_tols = false)
     )
     L = length(H)
@@ -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)(rightunit(sector) => 1)
     for i in 1:(middle_site - 1)
         P = physicalspace(H, i)
         ALs[i] = isomorphism(T, left ⊗ P ← fuse(left ⊗ P))

src/algorithms/correlators.jl

@@ -16,7 +16,7 @@ function correlator(
     )
     first(js) > i || @error "i should be smaller than j ($i, $(first(js)))"
     S₁ = _firstspace(O₁)
-    S₁ == oneunit(S₁) || throw(ArgumentError("O₁ should start with a trivial leg."))
+    isunitspace(S₁) || throw(ArgumentError("O₁ should start with a trivial leg."))
     S₂ = _lastspace(O₂)
     S₂ == S₁' || throw(ArgumentError("O₂ should end with a trivial leg."))
 
@@ -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/chepigaansatz.jl

@@ -35,10 +35,10 @@ end
 
 function excitations(
         H, alg::ChepigaAnsatz, ψ::FiniteMPS, envs = environments(ψ, H);
-        sector = one(sectortype(ψ)), num::Int = 1, pos::Int = length(ψ) ÷ 2
+        sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) || throw(ArgumentError("invalid position $pos"))
-    sector == one(sector) || error("not yet implemented for charged excitations")
+    isunit(sector) || error("not yet implemented for charged excitations")
 
     # add random offset to kickstart Krylov process:
     AC = ψ.AC[pos]
@@ -100,10 +100,10 @@ end
 
 function excitations(
         H, alg::ChepigaAnsatz2, ψ::FiniteMPS, envs = environments(ψ, H);
-        sector = one(sectortype(ψ)), num::Int = 1, pos::Int = length(ψ) ÷ 2
+        sector = leftunit(ψ), num::Int = 1, pos::Int = length(ψ) ÷ 2
     )
     1 ≤ pos ≤ length(ψ) - 1 || throw(ArgumentError("invalid position $pos"))
-    sector == one(sector) || error("not yet implemented for charged excitations")
+    isunit(sector) || error("not yet implemented for charged excitations")
 
     # add random offset to kickstart Krylov process:
     @plansor AC2[-1 -2; -3 -4] := ψ.AC[pos][-1 -2; 1] * ψ.AR[pos + 1][1 -4; -3]

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
@@ -86,14 +85,14 @@ Create and optimise infinite quasiparticle states.
 # Keywords
 - `num::Int`: number of excited states to compute
 - `solver`: algorithm for the linear solver of the quasiparticle environments
-- `sector=one(sectortype(left_ψ))`: charge of the quasiparticle state
+- `sector=leftunit(left_ψ)`: charge of the quasiparticle state
 - `parallel=true`: enable multi-threading over different momenta
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, momentum::Number, lmps::InfiniteMPS,
         lenvs = environments(lmps, H), rmps::InfiniteMPS = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), kwargs...
+        sector = leftunit(lmps), kwargs...
     )
     ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector, momentum)
     return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...)
@@ -103,15 +102,27 @@ function excitations(
         lenvs = environments(lmps, H), rmps = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
         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),
-        }
+        sector = leftunit(lmps), parallel = true, kwargs...
     )
+    # 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
@@ -167,13 +178,13 @@ Create and optimise finite quasiparticle states.
 
 # Keywords
 - `num::Int`: number of excited states to compute
-- `sector=one(sectortype(left_ψ))`: charge of the quasiparticle state
+- `sector=leftunit(lmps)`: charge of the quasiparticle state
 """
 function excitations(
         H, alg::QuasiparticleAnsatz, lmps::FiniteMPS,
         lenvs = environments(lmps, H), rmps::FiniteMPS = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), num = 1
+        sector = leftunit(lmps), num = 1
     )
     ϕ₀ = LeftGaugedQP(rand, lmps, rmps; sector)
     return excitations(H, alg, ϕ₀, lenvs, renvs; num)
@@ -262,7 +273,7 @@ function excitations(
         H::MultilineMPO, alg::QuasiparticleAnsatz, momentum::Real, lmps::MultilineMPS,
         lenvs = environments(lmps, H), rmps = lmps,
         renvs = lmps === rmps ? lenvs : environments(rmps, H);
-        sector = one(sectortype(lmps)), kwargs...
+        sector = leftunit(lmps), kwargs...
     )
     ϕ₀ = LeftGaugedQP(randn, lmps, rmps; sector, momentum)
     return excitations(H, alg, ϕ₀, lenvs, renvs; kwargs...)

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/timestep/taylorcluster.jl

@@ -186,14 +186,14 @@ function make_time_mpo(
     H′ = copy(parent(H))
 
     V_left = left_virtualspace(H[1])
-    V_left′ = ⊞(V_left, oneunit(V_left), oneunit(V_left))
+    V_left′ = ⊞(V_left, rightunitspace(V_left), rightunitspace(V_left))
     H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1]))
     for (I, v) in nonzero_pairs(H[1])
         H′[1][I] = v
     end
 
     V_right = right_virtualspace(H[end])
-    V_right′ = ⊞(oneunit(V_right), oneunit(V_right), V_right)
+    V_right′ = ⊞(rightunitspace(V_right), rightunitspace(V_right), V_right)
     H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′)
     for (I, v) in nonzero_pairs(H[end])
         H′[end][I[1], 1, 1, end] = v

src/algorithms/timestep/wii.jl

@@ -94,14 +94,14 @@ function make_time_mpo(
     H′ = copy(parent(H))
 
     V_left = left_virtualspace(H[1])
-    V_left′ = ⊞(V_left, oneunit(V_left), oneunit(V_left))
+    V_left′ = ⊞(V_left, rightunitspace(V_left), rightunitspace(V_left))
     H′[1] = similar(H[1], V_left′ ⊗ space(H[1], 2) ← domain(H[1]))
     for (I, v) in nonzero_pairs(H[1])
         H′[1][I] = v
     end
 
     V_right = right_virtualspace(H[end])
-    V_right′ = ⊞(oneunit(V_right), oneunit(V_right), V_right)
+    V_right′ = ⊞(rightunitspace(V_right), rightunitspace(V_right), V_right)
     H′[end] = similar(H[end], codomain(H[end]) ← space(H[end], 3)' ⊗ V_right′)
     for (I, v) in nonzero_pairs(H[end])
         H′[end][I[1], 1, 1, end] = v

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 = first(left_virtualspace(H[1]))
     W = map(1:length(H)) do site
         return BraidingTensor{scalartype(H)}(physicalspace(H, site), V)
     end
@@ -78,7 +78,7 @@ end
 
 """
     transfer_spectrum(above::InfiniteMPS; below=above, tol=Defaults.tol, num_vals=20,
-                           sector=first(sectors(oneunit(left_virtualspace(above, 1)))))
+                           sector=leftunit(above))
 
 Calculate the partial spectrum of the left mixed transfer matrix corresponding to the
 overlap of a given `above` state and a `below` state. The `sector` keyword argument can be
@@ -89,7 +89,7 @@ domain of each eigenvector. The `tol` and `num_vals` keyword arguments are passe
 """
 function transfer_spectrum(
         above::InfiniteMPS; below = above, tol = Defaults.tol, num_vals = 20,
-        sector = first(sectors(oneunit(left_virtualspace(above, 1))))
+        sector = leftunit(above)
     )
     init = randomize!(
         similar(
@@ -275,13 +275,13 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O}, L = length(mpo)) wher
     V_wrap = left_virtualspace(mpo, 1)
     ST = storagetype(O)
 
-    util = isometry(storagetype(O), oneunit(V_wrap) ← one(V_wrap))
+    util = isometry(storagetype(O), rightunitspace(V_wrap) ← one(V_wrap))
     @plansor cup[-1; -2 -3] := id(ST, V_wrap)[-2; -3] * util[-1]
 
     local F_right
     for i in 1:L
-        V_left = i == 1 ? oneunit(V_wrap) : fuse(V_wrap ⊗ left_virtualspace(mpo, i))
-        V_right = i == L ? oneunit(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
+        V_left = i == 1 ? rightunitspace(V_wrap) : fuse(V_wrap ⊗ left_virtualspace(mpo, i))
+        V_right = i == L ? rightunitspace(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
         output[i] = similar(
             mpo[i], V_left * physicalspace(mpo, i) ← physicalspace(mpo, i) * V_right
         )
@@ -336,7 +336,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L = length(H))
     output = Vector{O}(undef, L)
     for site in 1:L
         V_left = if site == 1
-            oneunit(V_wrap)
+            rightunitspace(V_wrap)
         else
             vs = Vector{S}(undef, chi_)
             for (k, v) in indmap
@@ -345,7 +345,7 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L = length(H))
             SumSpace(vs)
         end
         V_right = if site == L
-            oneunit(V_wrap)
+            rightunitspace(V_wrap)
         else
             vs = Vector{S}(undef, chi_)
             for (k, v) in indmap

src/environments/abstract_envs.jl

@@ -81,7 +81,7 @@ function environment_alg(
         tol = Defaults.tol, maxiter = Defaults.maxiter, krylovdim = Defaults.krylovdim,
         verbosity = Defaults.VERBOSE_NONE
     )
-    max_krylovdim = dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1))
+    max_krylovdim = ceil(Int, dim(left_virtualspace(above, 1)) * dim(left_virtualspace(below, 1)))
     return GMRES(; tol, maxiter, krylovdim = min(max_krylovdim, krylovdim), verbosity)
 end
 function environment_alg(

src/operators/abstractmpo.jl

@@ -31,6 +31,9 @@ right_virtualspace(mpo::AbstractMPO) = map(Base.Fix1(right_virtualspace, mpo), e
 physicalspace(mpo::AbstractMPO, site::Int) = physicalspace(mpo[site])
 physicalspace(mpo::AbstractMPO) = map(Base.Fix1(physicalspace, mpo), eachsite(mpo))
 
+TensorKit.leftunit(mpo::AbstractMPO) = only(sectors(leftunitspace(physicalspace(mpo, 1))))
+TensorKit.rightunit(mpo::AbstractMPO) = only(sectors(rightunitspace(physicalspace(mpo, 1))))
+
 for ftype in (:spacetype, :sectortype, :storagetype)
     @eval TensorKit.$ftype(mpo::AbstractMPO) = $ftype(typeof(mpo))
     @eval TensorKit.$ftype(::Type{MPO}) where {MPO <: AbstractMPO} = $ftype(eltype(MPO))

src/operators/mpo.jl

@@ -29,7 +29,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
 
 """
@@ -223,8 +223,9 @@ function Base.:*(mpo1::FiniteMPO{<:MPOTensor}, mpo2::FiniteMPO{<:MPOTensor})
     N = check_length(mpo1, mpo2)
     (S = spacetype(mpo1)) == spacetype(mpo2) || throw(SectorMismatch())
 
-    if (left_virtualspace(mpo1, 1) != oneunit(S) || left_virtualspace(mpo2, 1) != oneunit(S)) ||
-            (right_virtualspace(mpo1, N) != oneunit(S) || right_virtualspace(mpo2, N) != oneunit(S))
+
+    if (isunitspace(left_virtualspace(mpo1, 1)) || isunitspace(left_virtualspace(mpo2, 1))) ||
+            (isunitspace(right_virtualspace(mpo1, N)) || isunitspace(right_virtualspace(mpo2, N)))
         @warn "left/right virtual space is not trivial, fusion may not be unique"
         # this is a warning because technically any isomorphism that fuses the left/right
         # would work and for now I dont feel like figuring out if this is important