SpaceMismatch in correlation_length from CTMRG

[ogauthe]

Hi! I am trying to compute the correlation length of a wavefunction (say AKLT spin 2) from a CTMRG environment. I get a SpaceMismatch. Not sure whether I am doing something wrong. Here is my code:

using TensorKit: ←, ⊗, SU2Irrep, TensorMap, Vect, permute, truncdim, truncbelow
using MPSKit: correlation_length, leading_boundary
using PEPSKit: CTMRGEnv, InfinitePEPS

vD = Vect[SU2Irrep](1//2 => 1)
vd = Vect[SU2Irrep](2 => 1)
tAKLT_A = TensorMap(ones, vd ← vD ⊗ vD ⊗ vD ⊗ vD)
tAKLT_B = permute(adjoint(tAKLT_A), (5,), (1, 2, 3, 4))

ψ = InfinitePEPS([tAKLT_A tAKLT_B; tAKLT_B tAKLT_A])
trscheme = truncdim(100) & truncbelow(1e-12)
boundary_alg = (; tol=1e-10, trscheme=trscheme, maxiter=100);
env0 = CTMRGEnv(randn, Float64, ψ, oneunit(vD));
env, info_ctmrg = leading_boundary(env0, ψ; boundary_alg...);

ξ = correlation_length(ψ, env)  # SpaceMismatch

ERROR: SpaceMismatch("(Rep[SU₂](0=>6, 1=>13, 2=>8, 3=>2) ⊗ Rep[SU₂](1/2=>1) ⊗ Rep[SU₂](1/2=>1)') ≠ (Rep[SU₂](0=>6, 1=>13, 2=>8, 3=>2) ⊗ Rep[SU₂](1/2=>1)' ⊗ Rep[SU₂](1/2=>1))")
Stacktrace:
  [1] compose(W::TensorKit.TensorMapSpace{…}, V::TensorKit.TensorMapSpace{…})
    @ TensorKit ~/.julia/packages/TensorKit/hkxhv/src/spaces/homspace.jl:172
  [2] tensorcontract_structure(A::TensorMap{…}, pA::Tuple{…}, conjA::Bool, B::TensorMap{…}, pB::Tuple{…}, conjB::Bool, pAB::Tuple{…})
    @ TensorKit ~/.julia/packages/TensorKit/hkxhv/src/tensors/tensoroperations.jl:128
  [3] tensoralloc_contract
    @ ~/.julia/packages/TensorOperations/jVddt/src/implementation/allocator.jl:101 [inlined]
  [4] tensoralloc_contract
    @ ~/.julia/packages/TensorOperations/jVddt/src/implementation/allocator.jl:100 [inlined]
  [5] macro expansion
    @ ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfer.jl:0 [inlined]
  [6] transfer_left(v::TensorMap{…}, A::TensorMap{…}, Abar::TensorMap{…})
    @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfer.jl:18
  [7] transfer_left
    @ ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfer.jl:101 [inlined]
  [8] SingleTransferMatrix
    @ ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfermatrix.jl:50 [inlined]
  [9] #208
    @ ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfermatrix.jl:47 [inlined]
 [10] FlipArgs
    @ ./reduce.jl:214 [inlined]
 [11] BottomRF
    @ ./reduce.jl:86 [inlined]
 [12] _foldl_impl(op::Base.BottomRF{…}, init::TensorMap{…}, itr::Base.Iterators.Reverse{…})
    @ Base ./reduce.jl:58
 [13] foldl_impl
    @ ./reduce.jl:48 [inlined]
 [14] mapfoldr_impl
    @ ./reduce.jl:204 [inlined]
 [15] mapfoldr
    @ ./reduce.jl:223 [inlined]
 [16] foldr
    @ ./reduce.jl:242 [inlined]
 [17] ProductTransferMatrix
    @ ~/.julia/packages/MPSKit/bM5WR/src/transfermatrix/transfermatrix.jl:47 [inlined]
 [18] apply
    @ ~/.julia/packages/KrylovKit/jC5gU/src/apply.jl:2 [inlined]
 [19] initialize(iter::KrylovKit.ArnoldiIterator{…}; verbosity::Int64)
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/factorizations/arnoldi.jl:141
 [20] initialize
    @ ~/.julia/packages/KrylovKit/jC5gU/src/factorizations/arnoldi.jl:136 [inlined]
 [21] _schursolve(A::MPSKit.ProductTransferMatrix{…}, x₀::TensorMap{…}, howmany::Int64, which::Symbol, alg::KrylovKit.Arnoldi{…})
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/eigsolve/arnoldi.jl:360
 [22] eigsolve(A::MPSKit.ProductTransferMatrix{…}, x₀::TensorMap{…}, howmany::Int64, which::Symbol, alg::KrylovKit.Arnoldi{…}; alg_rrule::KrylovKit.Arnoldi{…})
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/eigsolve/arnoldi.jl:148
 [23] eigsolve
    @ ~/.julia/packages/KrylovKit/jC5gU/src/eigsolve/arnoldi.jl:147 [inlined]
 [24] eigsolve(f::MPSKit.ProductTransferMatrix{…}, x₀::TensorMap{…}, howmany::Int64, which::Symbol; kwargs::@Kwargs{…})
    @ KrylovKit ~/.julia/packages/KrylovKit/jC5gU/src/eigsolve/eigsolve.jl:223
 [25] eigsolve
    @ ~/.julia/packages/KrylovKit/jC5gU/src/eigsolve/eigsolve.jl:197 [inlined]
 [26] transfer_spectrum(above::MPSKit.InfiniteMPS{…}; below::MPSKit.InfiniteMPS{…}, tol::Float64, num_vals::Int64, sector::SU2Irrep)
    @ MPSKit ~/.julia/packages/MPSKit/bM5WR/src/algorithms/toolbox.jl:89
 [27] transfer_spectrum
    @ ~/.julia/packages/MPSKit/bM5WR/src/algorithms/toolbox.jl:82 [inlined]
 [28] (::PEPSKit.var"#507#511"{Int64, @Kwargs{}, CTMRGEnv{TensorMap{…}, TensorMap{…}}})(r::Int64)
    @ PEPSKit ~/.julia/packages/PEPSKit/2FBaz/src/algorithms/toolbox.jl:289
 [29] iterate
    @ ./generator.jl:48 [inlined]
 [30] _collect(c::UnitRange{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base ./array.jl:811
 [31] collect_similar
    @ ./array.jl:720 [inlined]
 [32] map
    @ ./abstractarray.jl:3371 [inlined]
 [33] _correlation_length(env::CTMRGEnv{TensorMap{…}, TensorMap{…}}; num_vals::Int64, kwargs::@Kwargs{})
    @ PEPSKit ~/.julia/packages/PEPSKit/2FBaz/src/algorithms/toolbox.jl:286
 [34] _correlation_length
    @ ~/.julia/packages/PEPSKit/2FBaz/src/algorithms/toolbox.jl:278 [inlined]
 [35] correlation_length(state::InfinitePEPS{TensorKit.AbstractTensorMap{…}}, env::CTMRGEnv{TensorMap{…}, TensorMap{…}})
    @ PEPSKit ~/.julia/packages/PEPSKit/2FBaz/src/algorithms/toolbox.jl:275
 [36] top-level scope
    @ REPL[57]:1
Some type information was truncated. Use `show(err)` to see complete types.

I have been able to get a value from a lower level approach

using MPSKit: InfiniteMPS, transfer_spectrum
above = InfiniteMPS([env.edges[1, 1, 1], env.edges[1, 1, 2]]);
below = InfiniteMPS([env.edges[3, 1, 1], env.edges[3, 1, 2]]);
tsp = transfer_spectrum(above; below=below, num_vals=2)
ξ2 = -length(above) / log(abs(real(tsp[2])))  # around 1

But the result is wrong from a factor 2 (AKLT spin 2 wavefunction has correlation length around 2 according to https://link.aps.org/doi/10.1103/PhysRevB.96.121115)

[lkdvos]

I think the wrong factor is fine, I would have to double check but probably it is because the correlation length is in units of unit cell, not actual sites.

The space mismatch is interesting, it seems to indicate we are misaligning something since your example has the interesting case where the physical spaces are actually different on each site: half of them are dual and half of them aren't.

Edit: thinking about this slightly more, probably we are not using the second row for the bottom of the transfer matrix , ie

PEPSKit.jl/src/algorithms/toolbox.jl

Line 288 in bb53c65

below = InfiniteMPS(_dag.(env.edges[SOUTH, r, :]))

has to be row plus one.