Modify mixed_fixpoints, add correlation_length

src/PEPSKit.jl

@@ -71,7 +71,7 @@ module Defaults
 end
 
 export SVDAdjoint, IterSVD, NonTruncSVDAdjoint
-export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv
+export FixedSpaceTruncation, ProjectorAlg, CTMRG, CTMRGEnv, correlation_length
 export LocalOperator
 export expectation_value, costfun
 export leading_boundary

src/algorithms/ctmrg.jl

@@ -328,3 +328,35 @@ function LinearAlgebra.norm(peps::InfinitePEPS, env::CTMRGEnv)
 
     return total
 end
+
+"""
+    correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2)
+
+Compute the PEPS correlation length based on the horizontal and vertical
+transfer matrices. Additionally the (normalized) eigenvalue spectrum is
+returned. Specify the number of computed eigenvalues with `howmany`.
+"""
+function MPSKit.correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2)
+    ξ_h = Vector{Float64}(undef, size(peps, 1))
+    ξ_v = Vector{Float64}(undef, size(peps, 2))
+    λ_h = Matrix{ComplexF64}(undef, size(peps, 1), howmany)
+    λ_v = Matrix{ComplexF64}(undef, size(peps, 2), howmany)
+
+    # Horizontal
+    above_h = MPSMultiline(map(r -> InfiniteMPS(env.edges[1, r, :]), 1:size(peps, 1)))
+    below_h = MPSMultiline(map(r -> InfiniteMPS(env.edges[3, r, :]), 1:size(peps, 1)))
+    transfer_peps_h = TransferPEPSMultiline(peps, NORTH)
+    _, (_, vals_h) = MPSKit.mixed_fixpoints(above_h, transfer_peps_h, below_h; howmany)
+    λ_h = map(λ_row -> λ_row / abs(λ_row[1]), vals_h)  # Normalize largest eigenvalue
+    ξ_h = map(λ_row -> -1 / log(abs(λ_row[2])), λ_h)
+
+    # Vertical
+    above_v = MPSMultiline(map(c -> InfiniteMPS(env.edges[2, :, c]), 1:size(peps, 2)))
+    below_v = MPSMultiline(map(c -> InfiniteMPS(env.edges[4, :, c]), 1:size(peps, 2)))
+    transfer_peps_v = TransferPEPSMultiline(peps, EAST)
+    _, (_, vals_v) = MPSKit.mixed_fixpoints(above_v, transfer_peps_v, below_v; howmany)
+    λ_v = map(λ_row -> λ_row / abs(λ_row[1]), vals_v)  # Normalize largest eigenvalue
+    ξ_v = map(λ_row -> -1 / log(abs(λ_row[2])), λ_v)
+
+    return ξ_h, ξ_v, λ_h, λ_v
+end

src/mpskit_glue/transferpeps_environments.jl

@@ -20,6 +20,7 @@ function MPSKit.mixed_fixpoints(
     below::MPSMultiline,
     init=gen_init_fps(above, O, below);
     solver=MPSKit.Defaults.eigsolver,
+    howmany=1
 )
     T = eltype(above)
 
@@ -30,6 +31,8 @@ function MPSKit.mixed_fixpoints(
     envtype = eltype(init[1])
     lefties = PeriodicArray{envtype,2}(undef, numrows, numcols)
     righties = PeriodicArray{envtype,2}(undef, numrows, numcols)
+    lvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
+    rvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
 
     @threads for cr in 1:numrows
         c_above = above[cr]
@@ -40,18 +43,20 @@ function MPSKit.mixed_fixpoints(
         @sync begin
             Threads.@spawn begin
                 E_LL = TransferMatrix($c_above.AL, $O[cr], $c_below.AL)
-                (_, Ls, convhist) = eigsolve(flip(E_LL), $L0, 1, :LM, $solver)
-                convhist.converged < 1 &&
+                (Lvals, Ls, convhist) = eigsolve(flip(E_LL), $L0, howmany, :LM, $solver)
+                convhist.converged < howmany &&
                     @info "left eigenvalue failed to converge $(convhist.normres)"
                 L0 = first(Ls)
+                lvals[cr] = Lvals[1:howmany]
             end
 
             Threads.@spawn begin
                 E_RR = TransferMatrix($c_above.AR, $O[cr], $c_below.AR)
-                (_, Rs, convhist) = eigsolve(E_RR, $R0, 1, :LM, $solver)
-                convhist.converged < 1 &&
+                (Rvals, Rs, convhist) = eigsolve(E_RR, $R0, howmany, :LM, $solver)
+                convhist.converged < howmany &&
                     @info "right eigenvalue failed to converge $(convhist.normres)"
                 R0 = first(Rs)
+                rvals[cr] = Rvals[1:howmany]
             end
         end
 
@@ -81,7 +86,7 @@ function MPSKit.mixed_fixpoints(
         end
     end
 
-    return (lefties, righties)
+    return (lefties, righties), (lvals, rvals)
 end
 
 function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSMultiline)