Add MPSKit.correlation_length method for PEPS (#57)

* Add `MPSKit.transfer_spectrum` and `MPSKit.correlation_length` method for PEPS types

* Add small comments regarding inconsistent row/column indexing

* Add small `correlation_length` checks to Heisenberg and TF Ising tests

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,42 @@ 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; num_vals=2)
+    T = scalartype(peps)
+    ξ_h = Vector{real(T)}(undef, size(peps, 1))
+    ξ_v = Vector{real(T)}(undef, size(peps, 2))
+    λ_h = Vector{Vector{T}}(undef, size(peps, 1))
+    λ_v = Vector{Vector{T}}(undef, size(peps, 2))
+
+    # Horizontal
+    above_h = MPSMultiline(map(r -> InfiniteMPS(env.edges[1, r, :]), 1:size(peps, 1)))
+    respaced_edges_h = map(zip(space.(env.edges)[1, :, :], env.edges[3, :, :])) do (V1, T3)
+        return TensorMap(T3.data, V1)
+    end
+    below_h = MPSMultiline(map(r -> InfiniteMPS(respaced_edges_h[r, :]), 1:size(peps, 1)))
+    transfer_peps_h = TransferPEPSMultiline(peps, NORTH)
+    vals_h = MPSKit.transfer_spectrum(above_h, transfer_peps_h, below_h; num_vals)
+    λ_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)))
+    respaced_edges_v = map(zip(space.(env.edges)[2, :, :], env.edges[4, :, :])) do (V2, T4)
+        return TensorMap(T4.data, V2)
+    end
+    below_v = MPSMultiline(map(c -> InfiniteMPS(respaced_edges_v[:, c]), 1:size(peps, 2)))
+    transfer_peps_v = TransferPEPSMultiline(peps, EAST)
+    vals_v = MPSKit.transfer_spectrum(above_v, transfer_peps_v, below_v; num_vals)
+    λ_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/transferpepo_environments.jl

@@ -30,7 +30,7 @@ function MPSKit.mixed_fixpoints(
     righties = PeriodicArray{envtype,2}(undef, numrows, numcols)
 
     @threads for cr in 1:numrows
-        c_above = above[cr]
+        c_above = above[cr]  # TODO: Update index convention to above[cr - 1]
         c_below = below[cr + 1]
 
         (L0, R0) = init[cr]
@@ -90,7 +90,7 @@ function gen_init_fps(above::MPSMultiline, O::TransferPEPOMultiline, below::MPSM
             rand,
             scalartype(T),
             left_virtualspace(below, cr + 1, 0) * prod(adjoint.(west_spaces(O[cr], 1))),
-            left_virtualspace(above, cr, 0),
+            left_virtualspace(above, cr, 0),  # TODO: Update index convention to above[cr - 1]
         )
         R0::T = TensorMap(
             rand,

src/mpskit_glue/transferpeps_environments.jl

@@ -32,7 +32,7 @@ function MPSKit.mixed_fixpoints(
     righties = PeriodicArray{envtype,2}(undef, numrows, numcols)
 
     @threads for cr in 1:numrows
-        c_above = above[cr]
+        c_above = above[cr]  # TODO: Update index convention to above[cr - 1]
         c_below = below[cr + 1]
 
         (L0, R0) = init[cr]
@@ -94,7 +94,7 @@ function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSM
             left_virtualspace(below, cr + 1, 0) *
             space(O[cr].top[1], 5)' *
             space(O[cr].bot[1], 5),
-            left_virtualspace(above, cr, 0),
+            left_virtualspace(above, cr, 0),  # TODO: Update index convention to above[cr - 1]
         )
         R0::T = TensorMap(
             rand,
@@ -107,3 +107,31 @@ function gen_init_fps(above::MPSMultiline, O::TransferPEPSMultiline, below::MPSM
         (L0, R0)
     end
 end
+
+function MPSKit.transfer_spectrum(
+    above::MPSMultiline,
+    O::TransferPEPSMultiline,
+    below::MPSMultiline,
+    init=gen_init_fps(above, O, below);
+    num_vals=2,
+    solver=MPSKit.Defaults.eigsolver,
+)
+    @assert size(above) == size(O)
+    @assert size(below) == size(O)
+
+    numrows = size(above, 2)
+    envtype = eltype(init[1])
+    eigenvals = Vector{Vector{scalartype(envtype)}}(undef, numrows)
+
+    @threads for cr in 1:numrows
+        L0, = init[cr]
+
+        E_LL = TransferMatrix(above[cr - 1].AL, O[cr], below[cr + 1].AL)  # Note that this index convention is different from above!
+        λ, _, convhist = eigsolve(flip(E_LL), L0, num_vals, :LM, solver)
+        convhist.converged < num_vals &&
+            @warn "correlation length failed to converge: normres = $(convhist.normres)"
+        eigenvals[cr] = λ
+    end
+
+    return eigenvals
+end

test/heisenberg.jl

@@ -31,8 +31,10 @@ env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, c
 
 # find fixedpoint
 result = fixedpoint(psi_init, H, opt_alg, env_init)
+λ_h, λ_v, = correlation_length(result.peps, result.env)
 
 @test result.E ≈ -0.6694421 atol = 1e-2
+@test all(@. λ_h > 0 && λ_v > 0)
 
 # same test but for 2x2 unit cell
 H_2x2 = square_lattice_heisenberg(; unitcell=(2, 2))
@@ -41,5 +43,7 @@ env_init_2x2 = leading_boundary(
     CTMRGEnv(psi_init_2x2, ComplexSpace(χenv)), psi_init_2x2, ctm_alg
 )
 result_2x2 = fixedpoint(psi_init_2x2, H_2x2, opt_alg, env_init_2x2)
+λ_h_2x2, λ_v_2x2, = correlation_length(result_2x2.peps, result_2x2.env)
 
 @test result_2x2.E ≈ 4 * result.E atol = 1e-2
+@test all(@. λ_h_2x2 > 0 && λ_v_2x2 > 0)

test/tf_ising.jl

@@ -41,6 +41,7 @@ env_init = leading_boundary(CTMRGEnv(psi_init, ComplexSpace(χenv)), psi_init, c
 
 # find fixedpoint
 result = fixedpoint(psi_init, H, opt_alg, env_init)
+λ_h, λ_v, = correlation_length(result.peps, result.env)
 
 # compute magnetization
 σx = TensorMap(scalartype(psi_init)[0 1; 1 0], ℂ^2, ℂ^2)
@@ -50,3 +51,10 @@ magn = expectation_value(result.peps, M, result.env)
 @test result.E ≈ e atol = 1e-2
 @test imag(magn) ≈ 0 atol = 1e-6
 @test abs(magn) ≈ mˣ atol = 5e-2
+
+# find fixedpoint in polarized phase and compute correlations lengths
+H_polar = square_lattice_tf_ising(; h=4.5)
+result_polar = fixedpoint(psi_init, H_polar, opt_alg, env_init)
+λ_h_polar, λ_v_polar, = correlation_length(result_polar.peps, result_polar.env)
+@test λ_h_polar < λ_h
+@test λ_v_polar < λ_v