Move correlation_length out of PR (add separate PR later)

src/PEPSKit.jl

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

src/algorithms/ctmrg.jl

@@ -353,45 +353,3 @@ 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`.
-"""
-# TODO: Rewrite this similar to gauge_fix using transfermatrix_fixedpoint
-function MPSKit.correlation_length(peps::InfinitePEPS, env::CTMRGEnv; howmany=2)
-    ξ = Array{Float64,3}(undef, (2, size(peps)...))  # First index picks horizontal or vertical direction
-    λ = Array{ComplexF64,4}(undef, (2, howmany, size(peps)...))
-    for r in 1:size(peps, 1), c in 1:size(peps, 2)
-        @autoopt @tensor transferh[χ_LT D_Lab D_Lbe χ_LB; χ_RT D_Rab D_Rbe χ_RB] :=
-            env.edges[NORTH, _prev(r, end), c][χ_LT D1 D2; χ_RT] *
-            peps[r, c][d; D1 D_Rab D3 D_Lab] *
-            conj(peps[r, c][d; D2 D_Rbe D4 D_Lbe]) *
-            env.edges[SOUTH, _next(r, end), c][χ_RB D3 D4; χ_LB]
-        @autoopt @tensor transferv[χ_TL D_Tab D_Tbe χ_TL; χ_BL D_Bab D_Bbe χ_BR] :=
-            env.edges[EAST, r, _next(c, end)][χ_TR D1 D2; χ_BR] *
-            peps[r, c][d; D_Tab D1 D_Bab D3] *
-            conj(peps[r, c][d; D_Tbe D2 D_Bbe D4]) *
-            env.edges[WEST, r, _prev(c, end)][χ_BL D3 D4; χ_TL]
-
-        function lintransfer(v, t)
-            @tensor v′[-1 -2 -3 -4] := t[-1 -2 -3 -4; 1 2 3 4] * v[1 2 3 4]
-            return v′
-        end
-
-        v₀h = Tensor(randn, scalartype(transferh), domain(transferh))
-        valsh, = eigsolve(v -> lintransfer(v, transferh), v₀h, howmany, :LM)
-        λ[1, :, r, c] = valsh[1:howmany] / abs(valsh[1])  # Normalize largest eigenvalue to 1
-        ξ[1, r, c] = -1 / log(abs(λ[1, 2, r, c]))
-
-        v₀v = Tensor(rand, scalartype(transferv), domain(transferv))
-        valsv, = eigsolve(v -> lintransfer(v, transferv), v₀v, howmany, :LM)
-        λ[2, :, r, c] = valsv[1:howmany] / abs(valsv[1])  # Normalize largest eigenvalue to 1
-        ξ[2, r, c] = -1 / log(abs(λ[2, 2, r, c]))
-    end
-
-    return ξ, λ
-end