Change sector convention for excitations

src/algorithms/ED.jl

@@ -31,47 +31,51 @@ equivalent to dense eigenvectors.
 
 """
 function exact_diagonalization(H::FiniteMPOHamiltonian;
-                               sector=first(sectors(oneunit(physicalspace(H, 1)))),
-                               len::Int=length(H), num::Int=1, which::Symbol=:SR,
+                               sector=one(sectortype(H)),
+                               num::Int=1, which::Symbol=:SR,
                                alg=Defaults.alg_eigsolve(; dynamic_tols=false))
-    left = ℂ[typeof(sector)](sector => 1)
-    right = oneunit(left)
+    L = length(H)
+    @assert L > 1 "FiniteMPOHamiltonian must have length > 1"
+    middle_site = (L >> 1) + 1
 
-    middle_site = Int(round(len / 2))
-
-    Ot = eltype(H[1])
-
-    mpst_type = tensormaptype(spacetype(Ot), 2, 1, storagetype(Ot))
-    mpsb_type = tensormaptype(spacetype(Ot), 1, 1, storagetype(Ot))
-    Cs = Vector{Union{Missing,mpsb_type}}(missing, len + 1)
-    ALs = Vector{Union{Missing,mpst_type}}(missing, len)
-    ARs = Vector{Union{Missing,mpst_type}}(missing, len)
-    ACs = Vector{Union{Missing,mpst_type}}(missing, len)
+    T = storagetype(eltype(H))
+    TA = tensormaptype(spacetype(H), 2, 1, T)
 
+    # fuse from left to right
+    ALs = Vector{Union{Missing,TA}}(missing, L)
+    left = oneunit(spacetype(H))
     for i in 1:(middle_site - 1)
-        ALs[i] = isomorphism(storagetype(Ot), left * physicalspace(H, i),
-                             fuse(left * physicalspace(H, i)))
-        left = _lastspace(ALs[i])'
+        P = physicalspace(H, i)
+        ALs[i] = isomorphism(T, left ⊗ P ← fuse(left ⊗ P))
+        left = right_virtualspace(ALs[i])
     end
-    for i in len:-1:(middle_site + 1)
-        ARs[i] = _transpose_front(isomorphism(storagetype(Ot),
-                                              fuse(right * physicalspace(H, i)'),
-                                              right * physicalspace(H, i)'))
-        right = _firstspace(ARs[i])
+
+    # fuse from right to left
+    ARs = Vector{Union{Missing,TA}}(missing, L)
+    right = spacetype(H)(sector => 1)
+    for i in reverse((middle_site + 1):L)
+        P = physicalspace(H, i)
+        ARs[i] = _transpose_front(isomorphism(T, fuse(right ⊗ P') ← right ⊗ P'))
+        right = left_virtualspace(ARs[i])
     end
-    ACs[middle_site] = randomize!(similar(H[1][1, 1, 1, 1],
-                                          left * physicalspace(H, middle_site) ← right))
-    norm(ACs[middle_site]) == 0 && throw(ArgumentError("invalid sector"))
-    normalize!(ACs[middle_site])
 
-    #construct the largest possible finite mps of that length
+    # center
+    ACs = Vector{Union{Missing,TA}}(missing, L)
+    P = physicalspace(H, middle_site)
+    ACs[middle_site] = rand!(similar(ALs[1], left ⊗ P ← right))
+
+    TB = tensormaptype(spacetype(H), 1, 1, T)
+    Cs = Vector{Union{Missing,TB}}(missing, L + 1)
     state = FiniteMPS(ALs, ARs, ACs, Cs)
     envs = environments(state, H)
 
-    #optimize the middle site. Because there is no truncation, this single site captures the entire possible hilbert space
-    H_ac = ∂∂AC(middle_site, state, H, envs) # this linear operator is now the actual full hamiltonian!
-    (vals, vecs, convhist) = eigsolve(H_ac, state.AC[middle_site], num, which, alg)
+    # optimize the middle site
+    # Because the MPS is full rank - this is equivalent to the full Hamiltonian
+    AC₀ = state.AC[middle_site]
+    H_ac = ∂∂AC(middle_site, state, H, envs)
+    vals, vecs, convhist = eigsolve(H_ac, AC₀, num, which, alg)
 
+    # repack eigenstates
     state_vecs = map(vecs) do v
         cs = copy(state)
         cs.AC[middle_site] = v

src/operators/abstractmpo.jl

@@ -261,7 +261,7 @@ end
 
 function add_physical_charge(O::MPOTensor, charge::Sector)
     sectortype(O) === typeof(charge) || throw(SectorMismatch())
-    auxspace = Vect[typeof(charge)](charge => 1)
+    auxspace = Vect[typeof(charge)](charge => 1)'
     F = fuser(scalartype(O), physicalspace(O), auxspace)
     @plansor O_charged[-1 -2; -3 -4] := F[-2; 1 2] *
                                         O[-1 1; 4 3] *
@@ -270,16 +270,14 @@ function add_physical_charge(O::MPOTensor, charge::Sector)
 end
 function add_physical_charge(O::BraidingTensor, charge::Sector)
     sectortype(O) === typeof(charge) || throw(SectorMismatch())
-    auxspace = Vect[typeof(charge)](charge => 1)
+    auxspace = Vect[typeof(charge)](charge => 1)'
     V = left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
         fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)
     return BraidingTensor{scalartype(O)}(V)
 end
 function add_physical_charge(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}, charge::Sector)
     sectortype(O) == typeof(charge) || throw(SectorMismatch())
-
-    auxspace = Vect[typeof(charge)](charge => 1)
-
+    auxspace = Vect[typeof(charge)](charge => 1)'
     Odst = similar(O,
                    left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
                    fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O))

src/states/finitemps.jl

@@ -234,13 +234,13 @@
 
     Vspaces[1] = left
     for k in 2:N
-        Vspaces[k] = infimum(fuse(Vspaces[k - 1], fusedPspaces[k]), maxVspaces[k - 1])
+        Vspaces[k] = infimum(fuse(Vspaces[k - 1], fusedPspaces[k - 1]), maxVspaces[k - 1])
         dim(Vspaces[k]) > 0 || @warn "no fusion channels available at site $k"
     end
 
     Vspaces[end] = right
-    for k in N:-1:2
-        Vspaces[k] = infimum(maxVspaces[k - 1], fuse(Vspaces[k + 1], dual(fusedPspaces[k])))
+    for k in reverse(2:N)
+        Vspaces[k] = infimum(Vspaces[k], fuse(Vspaces[k + 1], dual(fusedPspaces[k])))
         dim(Vspaces[k]) > 0 || @warn "no fusion channels available at site $k"
     end
 
@@ -383,18 +383,36 @@
     return ProductSpace{S}(space.(Ref(t), Base.front(Base.tail(TensorKit.allind(t)))))
 end
 
+"""
+    max_virtualspaces(ψ::FiniteMPS)
+    max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=oneunit(S), right=oneunit(S))
+
+Compute the maximal virtual spaces of a given finite MPS or its physical spaces.
+"""
+function max_virtualspaces(Ps::Vector{<:Union{S,CompositeSpace{S}}}; left=oneunit(S),
+                           right=oneunit(S)) where {S<:ElementarySpace}
+    Vs = similar(Ps, length(Ps) + 1)
+    Vs[1] = left
+    Vs[end] = right
+    for k in 2:length(Ps)
+        Vs[k] = fuse(Vs[k - 1], fuse(Ps[k - 1]))
+    end
+    for k in reverse(2:length(Ps))
+        Vs[k] = infimum(Vs[k], fuse(Vs[k + 1], dual(fuse(Ps[k]))))
+    end
+    return Vs
+end
+function max_virtualspaces(ψ::FiniteMPS)
+    return max_virtualspaces(physicalspace(ψ); left=left_virtualspace(ψ, 1),
+                             right=right_virtualspace(ψ, length(ψ)))
+end
+
 """
     max_Ds(ψ::FiniteMPS) -> Vector{Float64}
 
 Compute the dimension of the maximal virtual space at a given site.
 """
-function max_Ds(ψ::FiniteMPS)
-    N = length(ψ)
-    physicaldims = dim.(space.(Ref(ψ)), 1:N)
-    D_left = cumprod(vcat(dim(left_virtualspace(ψ, 1)), physicaldims))
-    D_right = cumprod(vcat(dim(right_virtualspace(ψ, N)), reverse(physicaldims)))
-    return min.(D_left, D_right)
-end
+max_Ds(ψ::FiniteMPS) = dim.(max_virtualspaces(ψ))
 
 function Base.summary(io::IO, ψ::FiniteMPS)
     return print(io, "$(length(ψ))-site FiniteMPS ($(scalartype(ψ)), $(spacetype(ψ)))")

src/states/quasiparticle_state.jl

@@ -41,10 +41,11 @@ function LeftGaugedQP(datfun, left_gs, right_gs=left_gs;
     Xs = map(enumerate(VLs)) do (loc, vl)
         x = similar(vl,
                     right_virtualspace(vl) ←
-                    excitation_space' ⊗ right_virtualspace(right_gs, loc))
+                    excitation_space ⊗ right_virtualspace(right_gs, loc))
         fill_data!(x, datfun)
         return x
     end
+    sum(dim, Xs) == 0 && @warn "LeftGaugedQP: No possible fusion channels"
     left_gs isa InfiniteMPS ||
         momentum == zero(momentum) ||
         @warn "momentum is ignored for finite quasiparticles"
@@ -63,13 +64,15 @@ end
 
 function RightGaugedQP(datfun, left_gs, right_gs=left_gs;
                        sector=one(sectortype(left_gs)), momentum=0.0)
-    #find the left null spaces for the TNS
+    # find the left null spaces for the TNS
     excitation_space = Vect[typeof(sector)](sector => 1)
-    VRs = [adjoint(leftnull(adjoint(v))) for v in _transpose_tail.(right_gs.AR)]
-    Xs = [TensorMap{scalartype(left_gs)}(undef, left_virtualspace(right_gs, loc)',
-                                         excitation_space' * _firstspace(VRs[loc]))
-          for loc in 1:length(left_gs)]
-    fill_data!.(Xs, datfun)
+    VRs = convert(Vector, map(rightnull! ∘ _transpose_tail, right_gs.AR))
+    Xs = map(enumerate(VRs)) do (i, vr)
+        x = similar(vr,
+                    left_virtualspace(left_gs, i)' ←
+                    excitation_space ⊗ _firstspace(vr))
+        return fill_data!(x, datfun)
+    end
     left_gs isa InfiniteMPS ||
         momentum == zero(momentum) ||
         @warn "momentum is ignored for finite quasiparticles"
@@ -84,10 +87,10 @@ function Base.similar(v::RightGaugedQP, ::Type{T}=scalartype(v)) where {T<:Numbe
     return RightGaugedQP(v.left_gs, v.right_gs, similar.(v.Xs, T), v.VRs, v.momentum)
 end
 
-Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)];
+Base.getindex(v::LeftGaugedQP, i::Int) = v.VLs[mod1(i, end)] * v.Xs[mod1(i, end)]
 function Base.getindex(v::RightGaugedQP, i::Int)
     @plansor t[-1 -2; -3 -4] := v.Xs[mod1(i, end)][-1; -3 1] * v.VRs[mod1(i, end)][1; -4 -2]
-end;
+end
 
 function Base.setindex!(v::LeftGaugedQP, B, i::Int)
     v.Xs[mod1(i, end)] = v.VLs[mod1(i, end)]' * B

test/algorithms.jl

@@ -11,6 +11,7 @@ using MPSKit
 using MPSKit: fuse_mul_mpo
 using TensorKit
 using TensorKit: ℙ
+using LinearAlgebra: eigvals
 
 verbosity_full = 5
 verbosity_conv = 1
@@ -917,4 +918,50 @@ end
     @test Z_tdvp ≈ Z_dense_2 atol = 1e-2
 end
 
+@testset "Sector conventions" begin
+    L = 4
+    H = TestSetup.XY_model(U1Irrep; L)
+
+    H_dense = convert(TensorMap, H)
+    vals_dense = eigvals(H_dense)
+
+    tol = 1e-18 # tolerance required to separate degenerate eigenvalues
+    alg = MPSKit.Defaults.alg_eigsolve(; dynamic_tols=false, tol)
+
+    maxVspaces = MPSKit.max_virtualspaces(physicalspace(H))
+    gs, = find_groundstate(FiniteMPS(physicalspace(H), maxVspaces[2:(end - 1)]), H;
+                           verbosity=0)
+    E₀ = expectation_value(gs, H)
+    @test E₀ ≈ first(vals_dense[one(U1Irrep)])
+
+    for (sector, vals) in vals_dense
+        # ED tests
+        num = length(vals)
+        E₀s, ψ₀s, info = exact_diagonalization(H; num, sector, alg)
+        @test E₀s[1:num] ≈ vals[1:num]
+        # this is a trick to make the mps full-rank again, which is not guaranteed by ED
+        ψ₀ = changebonds(first(ψ₀s), SvdCut(; trscheme=notrunc()))
+        Vspaces = left_virtualspace.(Ref(ψ₀), 1:L)
+        push!(Vspaces, right_virtualspace(ψ₀, L))
+        @test all(splat(==), zip(Vspaces, MPSKit.max_virtualspaces(ψ₀)))
+
+        # Quasiparticle tests
+        Es, Bs = excitations(H, QuasiparticleAnsatz(; tol), gs; sector, num=1)
+        Es = Es .+ E₀
+        # first excited state is second eigenvalue if sector is trivial
+        @test Es[1] ≈ vals[isone(sector) ? 2 : 1] atol = 1e-8
+    end
+
+    # shifted charges tests
+    # targeting states with Sz = 1 => vals_shift_dense[0] == vals_dense[1]
+    # so effectively shifting the charges by -1
+    H_shift = MPSKit.add_physical_charge(H, U1Irrep.([1, 0, 0, 0]))
+    H_shift_dense = convert(TensorMap, H_shift)
+    vals_shift_dense = eigvals(H_shift_dense)
+    for (sector, vals) in vals_dense
+        sector′ = only(sector ⊗ U1Irrep(-1))
+        @test vals ≈ vals_shift_dense[sector′]
+    end
+end
+
 end

test/excis.jl

@@ -0,0 +1,100 @@
+using MPSKit, TensorKit
+import LinearAlgebra.eigvals
+using MPSKit: max_virtualspaces
+
+V = U1Space(i => 1 for i in 0:1)
+O = randn(V^2 ← V^2)
+O += O'  # Hermitian
+
+N = 3
+H = FiniteMPOHamiltonian(fill(V, N), (i, i + 1) => O for i in 1:(N - 1));
+
+h = convert(TensorMap, H);
+
+sectors = collect(blocksectors(h))
+sec = U1Irrep(1)
+@assert sec in sectors
+num = 10
+vals1 = eigvals(block(h, sec))
+vals2 = eigvals(h)[sec]
+@assert vals1 ≈ vals2
+vals3, vecs = exact_diagonalization(H; sector=(sec), num);
+vals1
+vals3
+
+right = U1Space(sec => 1)
+max_virtualspaces(physicalspace(H); right)
+
+psi_full = rand(oneunit(V) * V^N ← right);
+dim(space(psi_full))
+psi = MPSKit.decompose_localmps(psi_full);
+left_virtualspace.(psi)
+
+left_virtualspace.(vecs[1].AL)
+right_virtualspace.(vecs[1].AL)
+max_virtualspaces(physicalspace(H); right)
+psi1 = FiniteMPS(physicalspace(H), max_virtualspaces(physicalspace(H))[2:(end - 1)])
+psi1, = find_groundstate(psi1, H);
+psi2 = FiniteMPS(physicalspace(H), max_virtualspaces(physicalspace(H); right)[2:(end - 1)];
+                 right)
+left_virtualspace.(psi2.AL)
+right_virtualspace.(psi2.AL)
+
+psi2, = find_groundstate(psi2, H);
+expectation_value(psi2, H)
+
+Es, Bs = excitations(H, QuasiparticleAnsatz(), FiniteMPS(psi); sector=sec);
+@inferred excitations(H, QuasiparticleAnsatz(), FiniteMPS(psi); sector=sec);
+
+Es .+ expectation_value(psi1, H)
+
+vals1
+vals3
+
+psi2
+using TestEnv
+using Test
+TestEnv.activate()
+include("setup.jl")
+using .TestSetup
+using TensorKit, MPSKit
+using MPSKit: Multiline
+using KrylovKit
+
+H = repeat(TestSetup.sixvertex(), 2)
+ψ = InfiniteMPS([ℂ^2, ℂ^2], [ℂ^10, ℂ^10])
+ψ, envs, _ = leading_boundary(ψ, H,
+                              VUMPS(; maxiter=400, tol=1e-10))
+energies, ϕs = @inferred excitations(H, QuasiparticleAnsatz(),
+                                     [0.0, Float64(pi / 2)], ψ,
+                                     envs; verbosity=0)
+@test abs(energies[1]) > abs(energies[2]) # has a minimum at pi/2
+alg = QuasiparticleAnsatz()
+ps = [0.0, Float64(pi / 2)]
+excitations(H, alg, ps, ψ, envs; verbosity=0);
+using Cthulhu
+Hm = convert(MultilineMPO, H);
+psim = convert(MultilineMPS, ψ);
+envs = environments(psim, Hm);
+excitations(Hm, alg, ps[1], psim, envs, psim, envs; verbosity=0);
+
+@descend excitations(Hm, alg, ps[1], psim, envs, psim, envs);
+
+qp = Multiline([LeftGaugedQP(rand, psim[1], psim[1]; sector=one(sectortype(psim[1])),
+                             momentum=ps[1])]);
+excitations(Hm, alg, qp, envs, envs; num=1);
+
+@descend excitations(Hm, alg, qp, envs, envs; num=1);
+@code_warntype excitations(Hm, alg, qp, envs, envs; num=1);
+
+Heff = MPSKit.EffectiveExcitationHamiltonian(H, envs[1], envs[1], fill(1.0, length(H)));
+Heffs = Multiline([Heff]);
+@code_warntype KrylovKit.apply(Heff, qp[1])
+@descend KrylovKit.apply(Heff, qp[1])
+
+KrylovKit.apply(Multiline([Heff]), qp);
+@code_warntype KrylovKit.apply(Heffs, qp);
+
+@descend eigsolve(Heffs, qp, 1, :LM, Lanczos());
+@code_warntype eigsolve(Heffs, qp, 1, :LM);
+@descend eigsolve(Heffs, qp, 1, :LM);