update left and right_virtualspace (#205)

Author: lkdvos

src/MPSKit.jl

@@ -60,8 +60,6 @@ export exact_diagonalization
 export TransferMatrix
 export transfer_left, transfer_right
 
-@deprecate virtualspace left_virtualspace # there is a possible ambiguity when C isn't square, necessitating specifying left or right virtualspace
-
 # Abstract type defs
 abstract type Algorithm end
 

src/algorithms/expval.jl

@@ -66,8 +66,8 @@ function expectation_value(ψ::AbstractMPS, (inds, O)::Pair)
         # left side
         T = storagetype(site_type(ψ))
         @plansor Vl[-1 -2; -3] := isomorphism(T,
-                                              left_virtualspace(ψ, sites[1] - 1),
-                                              left_virtualspace(ψ, sites[1] - 1))[-1; -3] *
+                                              left_virtualspace(ψ, sites[1]),
+                                              left_virtualspace(ψ, sites[1]))[-1; -3] *
                                   conj(Ut[-2])
 
         # middle
@@ -104,7 +104,7 @@ function expectation_value(ψ::InfiniteMPS, H::InfiniteMPOHamiltonian,
                            envs::AbstractMPSEnvironments=environments(ψ, H))
     return sum(1:length(ψ)) do i
         util = fill_data!(similar(ψ.AL[1], right_virtualspace(H, i)[end]), one)
-        @plansor GR[-1 -2; -3] := r_LL(ψ, i)[-1; -3] * conj(util[-2])
+        @plansor GR[-1 -2; -3] := r_LL(ψ, i)[-1; -3] * util[-2]
         return contract_mpo_expval(ψ.AL[i], leftenv(envs, i, ψ), H[i][:, 1, 1, end], GR)
     end
 end

src/algorithms/fidelity_susceptibility.jl

@@ -11,7 +11,7 @@ function fidelity_susceptibility(state::Union{FiniteMPS,InfiniteMPS}, H₀::T,
         Tos = LeftGaugedQP(rand, state)
         for (i, ac) in enumerate(state.AC)
             temp = ∂∂AC(i, state, V, venvs) * ac
-            help = fill_data!(similar(ac, utilleg(Tos)), one)
+            help = fill_data!(similar(ac, auxiliaryspace(Tos)), one)
             @plansor Tos[i][-1 -2; -3 -4] := temp[-1 -2; -4] * help[-3]
         end
 

src/algorithms/timestep/timeevmpo.jl

@@ -327,7 +327,7 @@ function make_time_mpo(H::InfiniteMPOHamiltonian{T}, dt, alg::WII) where {T}
         Vᵣ = right_virtualspace(H, i)[1:(end - 1)]
         P = physicalspace(H, i)
 
-        h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ')
+        h′ = similar(H[i], Vₗ ⊗ P ← P ⊗ Vᵣ)
         h′[2:end, 1, 1, 2:end] = WA[i]
         h′[2:end, 1, 1, 1] = WB[i]
         h′[1, 1, 1, 2:end] = WC[i]

src/algorithms/toolbox.jl

@@ -63,10 +63,10 @@ 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)))))
-    init = randomize!(similar(above.AL[1], left_virtualspace(below, 0),
-                              ℂ[typeof(sector)](sector => 1)' * left_virtualspace(above, 0)))
+    init = randomize!(similar(above.AL[1], left_virtualspace(below, 1),
+                              ℂ[typeof(sector)](sector => 1)' * left_virtualspace(above, 1)))
 
-    transferspace = fuse(left_virtualspace(above, 0) * left_virtualspace(below, 0)')
+    transferspace = fuse(left_virtualspace(above, 1) * left_virtualspace(below, 1)')
     num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values
     eigenvals, eigenvecs, convhist = eigsolve(flip(TransferMatrix(above.AL, below.AL)),
                                               init, num_vals, :LM; tol=tol)
@@ -239,26 +239,26 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O},
 
     # allocate output
     output = Vector{O}(undef, L)
-    V_wrap = left_virtualspace(mpo, 1)'
+    V_wrap = left_virtualspace(mpo, 1)
     ST = storagetype(O)
 
     util = fill!(similar(mpo[1], oneunit(V_wrap)), one(scalartype(O)))
-    @plansor cup[-1; -2 -3] := id(ST, V_wrap)[-3; -2] * util[-1]
+    @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_right = i == L ? oneunit(V_wrap) : fuse(V_wrap ⊗ right_virtualspace(mpo, i))
         output[i] = similar(mpo[i],
                             V_left * physicalspace(mpo, i) ←
                             physicalspace(mpo, i) * V_right)
         F_left = i == 1 ? cup : F_right
         F_right = i == L ? cup :
-                  isomorphism(ST, V_right ← V_wrap * right_virtualspace(mpo, i)')
-        @plansor output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
-                                           τ[-3 1; 4 3] *
-                                           mpo[i][2 -2; 3 5] *
-                                           conj(F_right[-4; 4 5])
+                  isomorphism(ST, V_right ← V_wrap' * right_virtualspace(mpo, i))
+        @plansor contractcheck = true output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
+                                                                τ[-3 1; 4 3] *
+                                                                mpo[i][2 -2; 3 5] *
+                                                                conj(F_right[-4; 4 5])
     end
 
     mpo isa SparseMPO && dropzeros!.(output) # the above process fills sparse mpos with zeros.

src/environments/abstract_envs.jl

@@ -14,12 +14,10 @@ Base.unlock(envs::AbstractMPSEnvironments) = unlock(envs.lock);
 # Allocating tensors
 # ------------------
 
-# TODO: fix the fucking left/right virtualspace bullshit
-# TODO: storagetype stuff
 function allocate_GL(bra::AbstractMPS, mpo::AbstractMPO, ket::AbstractMPS, i::Int)
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
-    V = left_virtualspace(bra, i - 1) ⊗ left_virtualspace(mpo, i)' ←
-        left_virtualspace(ket, i - 1)
+    V = left_virtualspace(bra, i) ⊗ left_virtualspace(mpo, i)' ←
+        left_virtualspace(ket, i)
     if V isa BlockTensorKit.TensorMapSumSpace
         TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
     else
@@ -30,7 +28,7 @@ end
 
 function allocate_GR(bra::AbstractMPS, mpo::AbstractMPO, ket::AbstractMPS, i::Int)
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
-    V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i)' ←
+    V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ←
         right_virtualspace(bra, i)
     if V isa BlockTensorKit.TensorMapSumSpace
         TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
@@ -42,8 +40,8 @@ end
 
 function allocate_GBL(bra::QP, mpo::AbstractMPO, ket::QP, i::Int)
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
-    V = left_virtualspace(bra.left_gs, i - 1) ⊗ left_virtualspace(mpo, i)' ←
-        auxiliaryspace(ket)' ⊗ left_virtualspace(ket.left_gs, i - 1)
+    V = left_virtualspace(bra, i) ⊗ left_virtualspace(mpo, i)' ←
+        auxiliaryspace(ket)' ⊗ left_virtualspace(ket, i)
     if V isa BlockTensorKit.TensorMapSumSpace
         TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
     else
@@ -54,8 +52,8 @@ end
 
 function allocate_GBR(bra::QP, mpo::AbstractMPO, ket::QP, i::Int)
     T = Base.promote_type(scalartype(bra), scalartype(mpo), scalartype(ket))
-    V = right_virtualspace(ket.right_gs, i) ⊗ right_virtualspace(mpo, i)' ←
-        auxiliaryspace(ket)' ⊗ right_virtualspace(bra.right_gs, i)
+    V = right_virtualspace(ket, i) ⊗ right_virtualspace(mpo, i) ←
+        auxiliaryspace(ket)' ⊗ right_virtualspace(bra, i)
     if V isa BlockTensorKit.TensorMapSumSpace
         TT = blocktensormaptype(spacetype(bra), numout(V), numin(V), T)
     else

src/environments/finite_envs.jl

@@ -34,30 +34,18 @@ function environments(below, operator, above, leftstart, rightstart)
                               rightenvs)
 end
 
-function environments(below::FiniteMPS{S}, O::DenseMPO, above=nothing) where {S}
-    N = length(below)
-    leftstart = isomorphism(storagetype(S),
-                            left_virtualspace(below, 0) ⊗ space(O[1], 1)' ←
-                            left_virtualspace(something(above, below), 0))
-    rightstart = isomorphism(storagetype(S),
-                             right_virtualspace(something(above, below), N) ⊗
-                             space(O[N], 4)' ←
-                             right_virtualspace(below, length(below)))
-    return environments(below, O, above, leftstart, rightstart)
-end
-
 function environments(below::FiniteMPS{S}, O::Union{FiniteMPO,FiniteMPOHamiltonian},
                       above=nothing) where {S}
-    Vl_bot = left_virtualspace(below, 0)
+    Vl_bot = left_virtualspace(below, 1)
     Vl_mid = left_virtualspace(O, 1)
-    Vl_top = isnothing(above) ? left_virtualspace(below, 0) : left_virtualspace(above, 0)
+    Vl_top = isnothing(above) ? left_virtualspace(below, 1) : left_virtualspace(above, 1)
     leftstart = isomorphism(storagetype(S), Vl_bot ⊗ Vl_mid' ← Vl_top)
 
     N = length(below)
     Vr_bot = right_virtualspace(below, N)
     Vr_mid = right_virtualspace(O, N)
     Vr_top = isnothing(above) ? right_virtualspace(below, N) : right_virtualspace(above, N)
-    rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid' ← Vr_bot)
+    rightstart = isomorphism(storagetype(S), Vr_top ⊗ Vr_mid ← Vr_bot)
 
     return environments(below, O, above, leftstart, rightstart)
 end

src/operators/abstractmpo.jl

@@ -156,7 +156,7 @@ Compute the mpo tensor that arises from multiplying MPOs.
 function fuse_mul_mpo(O1::MPOTensor, O2::MPOTensor)
     T = promote_type(scalartype(O1), scalartype(O2))
     F_left = fuser(T, left_virtualspace(O2), left_virtualspace(O1))
-    F_right = fuser(T, right_virtualspace(O2)', right_virtualspace(O1)')
+    F_right = fuser(T, right_virtualspace(O2), right_virtualspace(O1))
     @plansor O[-1 -2; -3 -4] := F_left[-1; 1 2] *
                                 O2[1 5; -3 3] *
                                 O1[2 -2; 5 4] *
@@ -203,7 +203,7 @@ function add_physical_charge(O::BraidingTensor, charge::Sector)
     sectortype(O) === typeof(charge) || throw(SectorMismatch())
     auxspace = Vect[typeof(charge)](charge => 1)
     V = left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
-        fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)'
+        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)
@@ -213,7 +213,7 @@ function add_physical_charge(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}, charge:
 
     Odst = similar(O,
                    left_virtualspace(O) ⊗ fuse(physicalspace(O), auxspace) ←
-                   fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O)')
+                   fuse(physicalspace(O), auxspace) ⊗ right_virtualspace(O))
     for (I, v) in nonzero_pairs(O)
         Odst[I] = add_physical_charge(v, charge)
     end

src/operators/mpo.jl

@@ -19,7 +19,7 @@ const FiniteMPO{O<:MPOTensor} = MPO{O,Vector{O}}
 
 function FiniteMPO(Os::AbstractVector{O}) where {O<:MPOTensor}
     for i in eachindex(Os)[1:(end - 1)]
-        dual(right_virtualspace(Os[i])) == left_virtualspace(Os[i + 1]) ||
+        right_virtualspace(Os[i]) == left_virtualspace(Os[i + 1]) ||
             throw(SpaceMismatch("unmatching virtual spaces at site $i"))
     end
     return FiniteMPO{O}(Os)
@@ -38,7 +38,7 @@ const InfiniteMPO{O<:MPOTensor} = MPO{O,PeriodicVector{O}}
 
 function InfiniteMPO(Os::AbstractVector{O}) where {O<:MPOTensor}
     for i in eachindex(Os)
-        dual(right_virtualspace(Os[i])) == left_virtualspace(Os[mod1(i + 1, end)]) ||
+        right_virtualspace(Os[i]) == left_virtualspace(Os[mod1(i + 1, end)]) ||
             throw(SpaceMismatch("umatching virtual spaces at site $i"))
     end
     return InfiniteMPO{O}(Os)
@@ -139,8 +139,8 @@ function Base.convert(::Type{TensorMap}, mpo::FiniteMPO)
     U_left = ones(scalartype(mpo), V_left)'
 
     V_right = right_virtualspace(mpo, length(mpo))
-    @assert V_right == oneunit(V_right)'
-    U_right = ones(scalartype(mpo), V_right')
+    @assert V_right == oneunit(V_right)
+    U_right = ones(scalartype(mpo), V_right)
 
     tensors = vcat(U_left, parent(mpo), U_right)
     indices = [[i, -i, -(2N - i + 1), i + 1] for i in 1:length(mpo)]
@@ -166,10 +166,10 @@ function Base.:+(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
     A = storagetype(TO)
 
     # left half
-    F₁ = isometry(A, (right_virtualspace(mpo1, 1) ⊕ right_virtualspace(mpo2, 1))',
-                  right_virtualspace(mpo1, 1)')
+    F₁ = isometry(A, (right_virtualspace(mpo1, 1) ⊕ right_virtualspace(mpo2, 1)),
+                  right_virtualspace(mpo1, 1))
     F₂ = leftnull(F₁)
-    @assert _lastspace(F₂) == right_virtualspace(mpo2, 1)
+    @assert _lastspace(F₂) == right_virtualspace(mpo2, 1)'
 
     @plansor O[-3 -1 -2; -4] := mpo1[1][-1 -2; -3 1] * conj(F₁[-4; 1]) +
                                 mpo2[1][-1 -2; -3 1] * conj(F₂[-4; 1])
@@ -184,10 +184,10 @@ function Base.:+(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
         @plansor O₂[-1 -2; -3 -4] := R[-1; 1] * F₂[1; 2] * mpo2[i][2 -2; -3 -4]
 
         # incorporate fusers from right side
-        F₁ = isometry(A, (right_virtualspace(mpo1, i) ⊕ right_virtualspace(mpo2, i))',
-                      right_virtualspace(mpo1, i)')
+        F₁ = isometry(A, (right_virtualspace(mpo1, i) ⊕ right_virtualspace(mpo2, i)),
+                      right_virtualspace(mpo1, i))
         F₂ = leftnull(F₁)
-        @assert _lastspace(F₂) == right_virtualspace(mpo2, i)
+        @assert _lastspace(F₂) == right_virtualspace(mpo2, i)'
         @plansor O[-3 -1 -2; -4] := O₁[-1 -2; -3 1] * conj(F₁[-4; 1]) +
                                     O₂[-1 -2; -3 1] * conj(F₂[-4; 1])
 
@@ -248,8 +248,8 @@ function Base.:*(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
     S = spacetype(TO)
     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))
+       (right_virtualspace(mpo1, N) != oneunit(S) ||
+        right_virtualspace(mpo2, N) != oneunit(S))
         @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
@@ -261,13 +261,8 @@ function Base.:*(mpo1::FiniteMPO{TO}, mpo2::FiniteMPO{TO}) where {TO}
     # note order of mpos: mpo1 * mpo2 * state -> mpo2 on top of mpo1
     local Fᵣ # trick to make Fᵣ defined in the loop
     for i in 1:N
-        Fₗ = i != 1 ? Fᵣ :
-             isomorphism(A, fuse(left_virtualspace(mpo2, i), left_virtualspace(mpo1, i)),
-                         left_virtualspace(mpo2, i) * left_virtualspace(mpo1, i))
-        Fᵣ = isomorphism(A,
-                         fuse(right_virtualspace(mpo2, i)', right_virtualspace(mpo1, i)'),
-                         right_virtualspace(mpo2, i)' * right_virtualspace(mpo1, i)')
-
+        Fₗ = i != 1 ? Fᵣ : fuser(A, left_virtualspace(mpo2, i), left_virtualspace(mpo1, i))
+        Fᵣ = fuser(A, right_virtualspace(mpo2, i), right_virtualspace(mpo1, i))
         @plansor O[i][-1 -2; -3 -4] := Fₗ[-1; 1 4] * mpo2[i][1 2; -3 3] *
                                        mpo1[i][4 -2; 2 5] *
                                        conj(Fᵣ[-4; 3 5])
@@ -284,11 +279,8 @@ function Base.:*(mpo::FiniteMPO, mps::FiniteMPS)
 
     local Fᵣ # trick to make Fᵣ defined in the loop
     for i in 1:length(mps)
-        Fₗ = i != 1 ? Fᵣ :
-             isomorphism(TT, fuse(left_virtualspace(A[i]), left_virtualspace(mpo, i)),
-                         left_virtualspace(A[i]) * left_virtualspace(mpo, i))
-        Fᵣ = isomorphism(TT, fuse(right_virtualspace(A[i])', right_virtualspace(mpo, i)'),
-                         right_virtualspace(A[i])' * right_virtualspace(mpo, i)')
+        Fₗ = i != 1 ? Fᵣ : fuser(TT, left_virtualspace(mps, i), left_virtualspace(mpo, i))
+        Fᵣ = fuser(TT, right_virtualspace(mps, i), right_virtualspace(mpo, i))
         A[i] = _fuse_mpo_mps(mpo[i], A[i], Fₗ, Fᵣ)
     end
 
@@ -298,12 +290,11 @@ function Base.:*(mpo::FiniteMPO, mps::FiniteMPS)
 end
 
 function Base.:*(mpo::InfiniteMPO, mps::InfiniteMPS)
-    check_length(mpo, mps)
+    L = check_length(mpo, mps)
     T = promote_type(scalartype(mpo), scalartype(mps))
-    fusers = PeriodicArray(map(mps.AL, mpo) do al, mp
-                               return fuser(T, _firstspace(al), _firstspace(mp))
-                           end)
-    As = map(1:length(mps)) do i
+    fusers = PeriodicArray(fuser.(T, left_virtualspace.(Ref(mps), 1:L),
+                                  left_virtualspace.(Ref(mpo), 1:L)))
+    As = map(1:L) do i
         return _fuse_mpo_mps(mpo[i], mps.AL[i], fusers[i], fusers[i + 1])
     end
     return changebonds(InfiniteMPS(As), SvdCut(; trscheme=notrunc()))
@@ -346,14 +337,14 @@ function TensorKit.dot(bra::FiniteMPS{T}, mpo::FiniteMPO, ket::FiniteMPS{T}) whe
     Nhalf = N ÷ 2
     # left half
     ρ_left = isomorphism(storagetype(T),
-                         left_virtualspace(bra, 0) ⊗ left_virtualspace(mpo, 1)',
-                         left_virtualspace(ket, 0))
+                         left_virtualspace(bra, 1) ⊗ left_virtualspace(mpo, 1)',
+                         left_virtualspace(ket, 1))
     T_left = TransferMatrix(ket.AL[1:Nhalf], mpo[1:Nhalf], bra.AL[1:Nhalf])
     ρ_left = ρ_left * T_left
 
     # right half
     ρ_right = isomorphism(storagetype(T),
-                          right_virtualspace(ket, N) ⊗ right_virtualspace(mpo, N)',
+                          right_virtualspace(ket, N) ⊗ right_virtualspace(mpo, N),
                           right_virtualspace(ket, length(ket)))
     T_right = TransferMatrix(ket.AR[(Nhalf + 1):end], mpo[(Nhalf + 1):end],
                              bra.AR[(Nhalf + 1):end])