From 283a62c3a4889e6e8ae3b384caac49418e8c6fad Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 12 Dec 2024 06:49:10 -0500 Subject: [PATCH] update `left` and `right_virtualspace` (#205) --- src/MPSKit.jl | 2 - src/algorithms/expval.jl | 6 +-- src/algorithms/fidelity_susceptibility.jl | 2 +- src/algorithms/timestep/timeevmpo.jl | 2 +- src/algorithms/toolbox.jl | 22 ++++----- src/environments/abstract_envs.jl | 16 +++---- src/environments/finite_envs.jl | 18 ++------ src/operators/abstractmpo.jl | 6 +-- src/operators/mpo.jl | 55 ++++++++++------------- src/operators/mpohamiltonian.jl | 18 ++++---- src/states/abstractmps.jl | 32 +++++++------ src/states/finitemps.jl | 37 ++++----------- src/states/infinitemps.jl | 33 +++++--------- src/states/quasiparticle_state.jl | 16 +++---- src/states/windowmps.jl | 4 +- src/utility/utility.jl | 4 +- test/algorithms.jl | 4 +- 17 files changed, 113 insertions(+), 164 deletions(-) diff --git a/src/MPSKit.jl b/src/MPSKit.jl index 6a613fb9c..57fa0d3c5 100644 --- a/src/MPSKit.jl +++ b/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 diff --git a/src/algorithms/expval.jl b/src/algorithms/expval.jl index 8f239c890..0f4caf52d 100644 --- a/src/algorithms/expval.jl +++ b/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 diff --git a/src/algorithms/fidelity_susceptibility.jl b/src/algorithms/fidelity_susceptibility.jl index 8a58a6d2a..fb8c3b38c 100644 --- a/src/algorithms/fidelity_susceptibility.jl +++ b/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 diff --git a/src/algorithms/timestep/timeevmpo.jl b/src/algorithms/timestep/timeevmpo.jl index 0f537cd26..27c3756cf 100644 --- a/src/algorithms/timestep/timeevmpo.jl +++ b/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] diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 8ba3d589e..05a724e20 100644 --- a/src/algorithms/toolbox.jl +++ b/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. diff --git a/src/environments/abstract_envs.jl b/src/environments/abstract_envs.jl index bc2a95323..21965322f 100644 --- a/src/environments/abstract_envs.jl +++ b/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 diff --git a/src/environments/finite_envs.jl b/src/environments/finite_envs.jl index 0f1e6bb94..bf2c2d4fa 100644 --- a/src/environments/finite_envs.jl +++ b/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 diff --git a/src/operators/abstractmpo.jl b/src/operators/abstractmpo.jl index 3e5556180..f0218fb8e 100644 --- a/src/operators/abstractmpo.jl +++ b/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 diff --git a/src/operators/mpo.jl b/src/operators/mpo.jl index 2db031c08..d55b66e0b 100644 --- a/src/operators/mpo.jl +++ b/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]) diff --git a/src/operators/mpohamiltonian.jl b/src/operators/mpohamiltonian.jl index f6b985a21..b6a029446 100644 --- a/src/operators/mpohamiltonian.jl +++ b/src/operators/mpohamiltonian.jl @@ -36,7 +36,7 @@ const FiniteMPOHamiltonian{O<:MPOTensor} = MPOHamiltonian{O,Vector{O}} function FiniteMPOHamiltonian(Ws::AbstractVector{O}) where {O<:MPOTensor} for i in eachindex(Ws)[1:(end - 1)] - dual(right_virtualspace(Ws[i])) == left_virtualspace(Ws[i + 1]) || + right_virtualspace(Ws[i]) == left_virtualspace(Ws[i + 1]) || throw(ArgumentError("The virtual spaces of the MPO tensors at site $i do not match.")) end return FiniteMPOHamiltonian{O}(Ws) @@ -46,7 +46,7 @@ const InfiniteMPOHamiltonian{O<:MPOTensor} = MPOHamiltonian{O,PeriodicVector{O}} function InfiniteMPOHamiltonian(Ws::AbstractVector{O}) where {O<:MPOTensor} for i in eachindex(Ws) - dual(right_virtualspace(Ws[i])) == left_virtualspace(Ws[mod1(i + 1, end)]) || + right_virtualspace(Ws[i]) == left_virtualspace(Ws[mod1(i + 1, end)]) || throw(ArgumentError("The virtual spaces of the MPO tensors at site $i do not match.")) end return InfiniteMPOHamiltonian{O}(Ws) @@ -152,7 +152,7 @@ function FiniteMPOHamiltonian(lattice::AbstractArray{<:VectorSpace}, V[key_R == 0 ? end : key_R] = if O isa Number virtualspaces[i][key_L] else - right_virtualspace(O)' + right_virtualspace(O) end end end @@ -246,9 +246,9 @@ function InfiniteMPOHamiltonian(lattice′::AbstractArray{<:VectorSpace}, @assert virtualspaces[i - 1][key_L] == left_virtualspace(O) end if ismissing(virtualspaces[i][key_R]) - virtualspaces[i][key_R] = right_virtualspace(O)' + virtualspaces[i][key_R] = right_virtualspace(O) else - @assert virtualspaces[i][key_R] == right_virtualspace(O)' + @assert virtualspaces[i][key_R] == right_virtualspace(O) end end end @@ -368,8 +368,8 @@ function Base.convert(::Type{TensorMap}, H::FiniteMPOHamiltonian) U_left = ones(scalartype(H), V_left)' V_right = right_virtualspace(H, length(H)) - @assert V_right == oneunit(V_right)' - U_right = ones(scalartype(H), V_right') + @assert V_right == oneunit(V_right) + U_right = ones(scalartype(H), V_right) tensors = vcat(U_left, parent(H), U_right) indices = [[i, -i, -(i + N), i + 1] for i in 1:length(H)] @@ -402,7 +402,7 @@ function Base.:+(H₁::MPOH, H₂::MPOH) where {MPOH<:MPOHamiltonian} Vᵣ = (!isinf && i == length(H)) ? Vᵣ₁ : BlockTensorKit.oplus(Vᵣ₁[1:(end - 1)], Vᵣ₂[2:end]) - W = similar(eltype(H), Vₗ ⊗ physicalspace(H₁, i) ← physicalspace(H₁, i) ⊗ Vᵣ') + W = similar(eltype(H), Vₗ ⊗ physicalspace(H₁, i) ← physicalspace(H₁, i) ⊗ Vᵣ) #= (Direct) sum of two hamiltonians in Jordan form: (1 C₁ D₁) (1 C₂ D₂) (1 C₁ C₂ D₁+D₂) @@ -528,7 +528,7 @@ function Base.:*(H::FiniteMPOHamiltonian, mps::FiniteMPS) # right to middle U = ones(scalartype(H), right_virtualspace(H, length(H))) - @plansor a[-1 -2; -3 -4] := A[end][-1 2; -3] * H[end][-2 -4; 2 1] * conj(U[1]) + @plansor a[-1 -2; -3 -4] := A[end][-1 2; -3] * H[end][-2 -4; 2 1] * U[1] L, Q = rightorth!(a; alg=LQ()) A′[end] = transpose(convert(TensorMap, Q), ((1, 3), (2,))) diff --git a/src/states/abstractmps.jl b/src/states/abstractmps.jl index 248ea6eac..2d8221a8a 100644 --- a/src/states/abstractmps.jl +++ b/src/states/abstractmps.jl @@ -156,34 +156,38 @@ TensorKit.sectortype(ψ::AbstractMPS) = sectortype(typeof(ψ)) TensorKit.sectortype(ψtype::Type{<:AbstractMPS}) = sectortype(site_type(ψtype)) """ - left_virtualspace(ψ::AbstractMPS, i::Int) + left_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)]) -Return the left virtual space of the bond tensor to the right of site `i`. This is -equivalent to the left virtual space of the left-gauged site tensor at site `i + 1`. +Return the virtual space of the bond to the left of sites `pos`. + +!!! warning + In rare cases, the gauge tensor on the virtual space might not be square, and as a result it + cannot always be guaranteed that `right_virtualspace(ψ, i - 1) == left_virtualspace(ψ, i)` """ function left_virtualspace end left_virtualspace(A::GenericMPSTensor) = space(A, 1) left_virtualspace(O::MPOTensor) = space(O, 1) """ - right_virtualspace(ψ::AbstractMPS, i::Int) + right_virtualspace(ψ::AbstractMPS, [pos=1:length(ψ)]) + +Return the virtual space of the bond to the right of site(s) `pos`. -Return the right virtual space of the bond tensor to the right of site `i`. This is -equivalent to the right virtual space of the right-gauged site tensor at site `i`. +!!! warning + In rare cases, the gauge tensor on the virtual space might not be square, and as a result it + cannot always be guaranteed that `right_virtualspace(ψ, i - 1) == left_virtualspace(ψ, i)` """ function right_virtualspace end -right_virtualspace(A::GenericMPSTensor) = space(A, numind(A)) -right_virtualspace(O::MPOTensor) = space(O, 4) +right_virtualspace(A::GenericMPSTensor) = space(A, numind(A))' +right_virtualspace(O::MPOTensor) = space(O, 4)' """ - physicalspace(ψ::AbstractMPS, i::Int) + physicalspace(ψ::AbstractMPS, [pos=1:length(ψ)]) Return the physical space of the site tensor at site `i`. """ function physicalspace end +physicalspace(A::MPSTensor) = space(A, 2) physicalspace(A::GenericMPSTensor) = prod(x -> space(A, x), 2:(numind(A) - 1)) -function physicalspace(O::MPOTensor) - pspace = space(O, 2) - # Disallow SumSpace in physical space - return pspace isa SumSpace ? only(pspace) : pspace -end +physicalspace(O::MPOTensor) = space(O, 2) +physicalspace(O::AbstractBlockTensorMap{<:Any,<:Any,2,2}) = only(space(O, 2)) diff --git a/src/states/finitemps.jl b/src/states/finitemps.jl index 02e0eb783..fa29395d5 100644 --- a/src/states/finitemps.jl +++ b/src/states/finitemps.jl @@ -266,41 +266,22 @@ function TensorKit.storagetype(::Union{MPS,Type{MPS}}) where {A,MPS<:FiniteMPS{A end function left_virtualspace(ψ::FiniteMPS, n::Integer) - if n > 0 && !ismissing(ψ.ALs[n]) - dual(_lastspace(ψ.ALs[n])) - elseif n < length(ψ.ALs) && !ismissing(ψ.ALs[n + 1]) - _firstspace(ψ.ALs[n + 1]) - else - _firstspace(ψ.CR[n]) - end + checkbounds(ψ, n) + return !ismissing(ψ.ALs[n]) ? left_virtualspace(ψ.ALs[n]) : + !ismissing(ψ.ARs[n]) ? left_virtualspace(ψ.ARs[n]) : + dual(_lastspace(ψ.CR[n - 1])) end function right_virtualspace(ψ::FiniteMPS, n::Integer) - if n > 0 && !ismissing(ψ.ARs[n]) - dual(_lastspace(ψ.ARs[n])) - elseif n < length(ψ.ARs) && !ismissing(ψ.ARs[n + 1]) - _firstspace(ψ.ARs[n + 1]) - else - dual(_lastspace(ψ.CR[n])) - end + checkbounds(ψ, n) + return !ismissing(ψ.ARs[n]) ? right_virtualspace(ψ.ARs[n]) : + !ismissing(ψ.ALs[n]) ? right_virtualspace(ψ.ALs[n]) : + _firstspace(ψ.CR[n]) end physicalspace(ψ::FiniteMPS) = physicalspace.(Ref(ψ), 1:length(ψ)) function physicalspace(ψ::FiniteMPS{<:GenericMPSTensor{<:Any,N}}, n::Integer) where {N} N == 1 && return ProductSpace{spacetype(ψ)}() - A = if !ismissing(ψ.ALs[n]) - ψ.ALs[n] - elseif !ismissing(ψ.ARs[n]) - ψ.ARs[n] # should never reach last case? - else - ψ.AC[n] # should never reach last case? - end # should never reach last case? - - if N == 2 - return space(A, 2) - else - return ProductSpace{spacetype(ψ),N - 1}(space.(Ref(A), - Base.front(Base.tail(TensorKit.allind(A))))) - end + return physicalspace(coalesce(ψ.ALs[n], ψ.ARs[n], ψ.ACs[n])) end TensorKit.space(ψ::FiniteMPS{<:MPSTensor}, n::Integer) = space(ψ.AC[n], 2) diff --git a/src/states/infinitemps.jl b/src/states/infinitemps.jl index 694f236d6..38ff3299b 100644 --- a/src/states/infinitemps.jl +++ b/src/states/infinitemps.jl @@ -232,27 +232,17 @@ Base.checkbounds(::Type{Bool}, ψ::InfiniteMPS, i::Integer) = true site_type(::Type{<:InfiniteMPS{A}}) where {A} = A bond_type(::Type{<:InfiniteMPS{<:Any,B}}) where {B} = B -left_virtualspace(ψ::InfiniteMPS, n::Integer) = _firstspace(ψ.CR[n]) -right_virtualspace(ψ::InfiniteMPS, n::Integer) = dual(_lastspace(ψ.CR[n])) - -function physicalspace(ψ::InfiniteMPS{<:GenericMPSTensor{<:Any,N}}, n::Integer) where {N} - if N == 1 - return ProductSpace{spacetype(ψ)}() - elseif N == 2 - return space(ψ.AL[n], 2) - else - return ProductSpace{spacetype(ψ),N - 1}(space.(Ref(ψ.AL[n]), - Base.front(Base.tail(TensorKit.allind(ψ.AL[n]))))) - end -end +left_virtualspace(ψ::InfiniteMPS, n::Integer) = left_virtualspace(ψ.AL[n]) +right_virtualspace(ψ::InfiniteMPS, n::Integer) = right_virtualspace(ψ.AL[n]) +physicalspace(ψ::InfiniteMPS, n::Integer) = physicalspace(ψ.AL[n]) physicalspace(ψ::InfiniteMPS) = PeriodicArray(map(Base.Fix1(physicalspace, ψ), 1:length(ψ))) -TensorKit.space(ψ::InfiniteMPS{<:MPSTensor}, n::Integer) = space(ψ.AC[n], 2) -function TensorKit.space(ψ::InfiniteMPS{<:GenericMPSTensor}, n::Integer) - t = ψ.AC[n] - S = spacetype(t) - return ProductSpace{S}(space.(Ref(t), Base.front(Base.tail(TensorKit.allind(t))))) -end +# TensorKit.space(ψ::InfiniteMPS{<:MPSTensor}, n::Integer) = space(ψ.AC[n], 2) +# function TensorKit.space(ψ::InfiniteMPS{<:GenericMPSTensor}, n::Integer) +# t = ψ.AC[n] +# S = spacetype(t) +# return ProductSpace{S}(space.(Ref(t), Base.front(Base.tail(TensorKit.allind(t))))) +# end TensorKit.norm(ψ::InfiniteMPS) = norm(ψ.AC[1]) function TensorKit.normalize!(ψ::InfiniteMPS) @@ -333,7 +323,7 @@ l_LR(ψ::InfiniteMPS, loc::Int=1) = ψ.CR[loc - 1]' Left dominant eigenvector of the `AL`-`AL` transfermatrix. """ function l_LL(ψ::InfiniteMPS{A}, loc::Int=1) where {A} - return isomorphism(storagetype(A), space(ψ.AL[loc], 1), space(ψ.AL[loc], 1)) + return isomorphism(storagetype(A), left_virtualspace(ψ, loc), left_virtualspace(ψ, loc)) end """ @@ -342,7 +332,8 @@ end Right dominant eigenvector of the `AR`-`AR` transfermatrix. """ function r_RR(ψ::InfiniteMPS{A}, loc::Int=length(ψ)) where {A} - return isomorphism(storagetype(A), domain(ψ.AR[loc]), domain(ψ.AR[loc])) + return isomorphism(storagetype(A), right_virtualspace(ψ, loc), + right_virtualspace(ψ, loc)) end """ diff --git a/src/states/quasiparticle_state.jl b/src/states/quasiparticle_state.jl index 098656e93..01a203bc3 100644 --- a/src/states/quasiparticle_state.jl +++ b/src/states/quasiparticle_state.jl @@ -48,7 +48,7 @@ function RightGaugedQP(datfun, left_gs, right_gs=left_gs; #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 - 1)', + 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) @@ -87,7 +87,8 @@ end function Base.convert(::Type{RightGaugedQP}, input::LeftGaugedQP{S}) where {S<:InfiniteMPS} rg = RightGaugedQP(zero, input.left_gs, input.right_gs; - sector=first(sectors(utilleg(input))), momentum=input.momentum) + sector=first(sectors(auxiliaryspace(input))), + momentum=input.momentum) len = length(input) #construct environments @@ -131,7 +132,7 @@ end function Base.convert(::Type{LeftGaugedQP}, input::RightGaugedQP{S}) where {S<:InfiniteMPS} lg = LeftGaugedQP(zero, input.left_gs, input.right_gs; - sector=first(sectors(utilleg(input))), momentum=input.momentum) + sector=first(sectors(auxiliaryspace(input))), momentum=input.momentum) len = length(input) lBs = [@plansor t[-1; -2 -3] := input[1][1 2; -2 -3] * @@ -177,13 +178,10 @@ const InfiniteQP{S<:InfiniteMPS,T1,T2} = QP{S,T1,T2} TensorKit.spacetype(::Union{QP{S},Type{<:QP{S}}}) where {S} = spacetype(S) TensorKit.sectortype(::Union{QP{S},Type{<:QP{S}}}) where {S} = sectortype(S) -left_virtualspace(state::QP, i::Int) = space(state.Xs[mod1(i, end)], 1) -function right_virtualspace(state::QP, i::Int) - return space(state.Xs[mod1(i, end)], numind(state.Xs[mod1(i, end)])) -end +left_virtualspace(state::QP, i::Int) = left_virtualspace(state.left_gs, i) +right_virtualspace(state::QP, i::Int) = right_virtualspace(state.right_gs, i) auxiliaryspace(state::QP) = space(state.Xs[1], 2) -utilleg(v::QP) = space(v.Xs[1], 2) Base.copy(a::QP) = copy!(similar(a), a) Base.copyto!(a::QP, b::QP) = copy!(a, b) function Base.copy!(a::T, b::T) where {T<:QP} @@ -260,7 +258,7 @@ function Base.convert(::Type{<:FiniteMPS}, v::QP{S}) where {S<:FiniteMPS} elt = scalartype(v) - utl = utilleg(v) + utl = auxiliaryspace(v) ou = oneunit(utl) utsp = ou ⊕ ou upper = isometry(storagetype(site_type(v.left_gs)), utsp, ou) diff --git a/src/states/windowmps.jl b/src/states/windowmps.jl index 7225ddb21..37c8ba8bc 100644 --- a/src/states/windowmps.jl +++ b/src/states/windowmps.jl @@ -43,7 +43,7 @@ struct WindowMPS{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractFiniteMPS function WindowMPS(ψₗ::InfiniteMPS{A,B}, ψₘ::FiniteMPS{A,B}, ψᵣ::InfiniteMPS{A,B}=copy(ψₗ)) where {A<:GenericMPSTensor, B<:MPSBondTensor} - left_virtualspace(ψₗ, 0) == left_virtualspace(ψₘ, 0) && + left_virtualspace(ψₗ, 1) == left_virtualspace(ψₘ, 1) && right_virtualspace(ψₘ, length(ψₘ)) == right_virtualspace(ψᵣ, length(ψₘ)) || throw(SpaceMismatch("Mismatch between window and environment virtual spaces")) return new{A,B}(ψₗ, ψₘ, ψᵣ) @@ -62,7 +62,7 @@ end function WindowMPS(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S, ψₗ::InfiniteMPS, ψᵣ::InfiniteMPS=ψₗ) where {S<:ElementarySpace} - ψₘ = FiniteMPS(f, elt, physspaces, maxvirtspace; left=left_virtualspace(ψₗ, 0), + ψₘ = FiniteMPS(f, elt, physspaces, maxvirtspace; left=left_virtualspace(ψₗ, 1), right=right_virtualspace(ψᵣ, length(physspaces))) return WindowMPS(ψₗ, ψₘ, ψᵣ) end diff --git a/src/utility/utility.jl b/src/utility/utility.jl index 694eed205..6b82a5701 100644 --- a/src/utility/utility.jl +++ b/src/utility/utility.jl @@ -175,6 +175,6 @@ function check_length(a, b...) return L end -function fuser(::Type{T}, V1::S, V2::S) where {T<:Number,S<:IndexSpace} - return isomorphism(Vector{T}, fuse(V1 ⊗ V2), V1 ⊗ V2) +function fuser(::Type{T}, V1::S, V2::S) where {T,S<:IndexSpace} + return isomorphism(T, fuse(V1 ⊗ V2), V1 ⊗ V2) end diff --git a/test/algorithms.jl b/test/algorithms.jl index 98acc0044..6534fbea6 100644 --- a/test/algorithms.jl +++ b/test/algorithms.jl @@ -561,7 +561,7 @@ end state_tr = changebonds(state_oe, SvdCut(; trscheme=truncdim(dim(Dspace)))) - @test dim(left_virtualspace(state_tr, 5)) < dim(right_virtualspace(state_oe, 5)) + @test dim(left_virtualspace(state_tr, 5)) < dim(left_virtualspace(state_oe, 5)) end @testset "MPSMultiline" begin @@ -583,7 +583,7 @@ end state_tr = changebonds(state_oe, SvdCut(; trscheme=truncdim(dim(Dspace)))) - @test dim(right_virtualspace(state_tr, 1, 1)) < + @test dim(left_virtualspace(state_tr, 1, 1)) < dim(left_virtualspace(state_oe, 1, 1)) end end