From ce32b745c175c8b0a2ea0008a781ce78cfb6ff47 Mon Sep 17 00:00:00 2001
From: Paul Brehmer <paul.brehmer@univie.ac.at>
Date: Tue, 16 Jul 2024 19:05:39 +0200
Subject: [PATCH] Modify mixed_fixpoints, add correlation_length

---
 src/PEPSKit.jl                               |  2 +-
 src/algorithms/ctmrg.jl                      | 32 ++++++++++++++++++++
 src/mpskit_glue/transferpeps_environments.jl | 15 ++++++---
 3 files changed, 43 insertions(+), 6 deletions(-)

diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl
index 1a3c941b..5c3695ba 100644
--- a/src/PEPSKit.jl
+++ b/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
diff --git a/src/algorithms/ctmrg.jl b/src/algorithms/ctmrg.jl
index 860257fb..924a91d2 100644
--- a/src/algorithms/ctmrg.jl
+++ b/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
diff --git a/src/mpskit_glue/transferpeps_environments.jl b/src/mpskit_glue/transferpeps_environments.jl
index fb4c12ad..9064c18d 100644
--- a/src/mpskit_glue/transferpeps_environments.jl
+++ b/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)