Compatibility with multifusion categories (#297)

* add_util_leg spinoff for MPOHamiltonian

* apply add_util_mpoleg to correlator

* use add_util_mpoleg in FiniteMPO

* use add_util_mpoleg in test on correlator

* avoid `oneunit(S)` in one constructor of `InfiniteMPOHamiltonian`

* avoid `oneunit(S)` in one constructor of `FiniteMPOHamiltonian`

* remove type `LocalOperator` in `_oneunit`

* remove braidings in `left/right_excitation_transfer_system` of `InfiniteMPOHamiltonian`

* typo in `right_excitation_transfer_system`

* rewrite `(In)FiniteMPOHamiltonian` and the way `oneunit` was avoided

* avoid `oneunit(spacetype(H))` evaluation

* super important typo

* remove comments

* get `return_type` in `excitations` to evaluate the `sector` kwarg

* typo

* manipulate `oneunit`s in addition of `MPOHamiltonian`s

* format

* simplify expression

* add comments

* get `infinite_temperature_density_matrix` working

* generalise identifying space in `infinite_temperature_density_matrix`

* make util spaces generic in `left/right_excitation_transfer_system`

* remove braiding in `InfiniteQP` environments + rename variables for consistency

* Revert "remove braiding in `InfiniteQP` environments + rename variables for consistency"

This reverts commit 2114cdc.

* add multifusion tests

* actually use my variable everywhere

* add Z2 Ising model

* compare to Z2 Ising model + add injectivity check

* remove defaulting to `Z2Irrep` for new operators

* fix kitaev model term (even if it's not used)

* format

* actually run the multifusion test

* import TestSetup module, remove unused deps + add variance tests

* specify allowed spin for Z2 `S_zz` + don't use deprecated call

* add variance tests for `FiniteMPS` + allow it to be negative (?) + add `FiniteExcited` test

* Revert "Merge branch 'master' of https://github.com/QuantumKitHub/MPSKit.jl into boris-MTK-compat"

This reverts commit 2c96ab1, reversing
changes made to a17874a.

* Revert "Revert "Merge branch 'master' of https://github.com/QuantumKitHub/MPSKit.jl into boris-MTK-compat""

This reverts commit fdd95c7.

* Revert "Merge branch 'master' of https://github.com/QuantumKitHub/MPSKit.jl into boris-MTK-compat"

This reverts commit 2c96ab1, reversing
changes made to 8cf041f.

* Revert "Revert "Merge branch 'master' of https://github.com/QuantumKitHub/MPSKit.jl into boris-MTK-compat""

This reverts commit 2f1108e.

* fix the messed up merges hopefully

* format

* actually format everything

* typo

* deal with float dimension

* `left/rightunit` for MPS and MPO

* using `left/rightunit(space)`

* test updates

* some `unitspace`s

* some `rightunitspace`s

* docstring to ED

* remove ED error test

* more `rightunitspace`s in tests

* even more `rightunitspace`s

* some `unit`s, `leftunit`s and `isunit`s

* some code suggestions and fixes

* update `add_util_leg` and delete `add_util_mpoleg`

* format

* relax only to first

* check if hamiltonians are diagonal

* `leftunitspace` where it feels spatially correct and allowed

* Update src/algorithms/correlators.jl

Co-authored-by: Lukas Devos <ldevos98@gmail.com>

* deal with `unitspace` within `MPOHamiltonian` constructors

* get `left/right_canicalize!` working

* address some comments

* unit space insertions and removals

* bring back the wrapper

---------

Co-authored-by: Lukas Devos <ldevos98@gmail.com>

Author: borisdevos

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

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 optimize 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,23 @@ 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 +174,13 @@ Create and optimize 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 +269,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, leftunitspace(V_left), leftunitspace(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, leftunitspace(V_left), leftunitspace(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(
@@ -276,13 +276,14 @@ 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))
+        # kept as rightunitspace, but might need to change if we consider off-diagonal MPOs
+        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
         )
@@ -337,7 +338,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)
+            leftunitspace(V_wrap)
         else
             vs = Vector{S}(undef, chi_)
             for (k, v) in indmap
@@ -346,7 +347,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