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..585e9b5f 100644 --- a/src/algorithms/ctmrg.jl +++ b/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 diff --git a/src/mpskit_glue/transferpepo_environments.jl b/src/mpskit_glue/transferpepo_environments.jl index e6f00161..f79d641b 100644 --- a/src/mpskit_glue/transferpepo_environments.jl +++ b/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, diff --git a/src/mpskit_glue/transferpeps_environments.jl b/src/mpskit_glue/transferpeps_environments.jl index fb4c12ad..ccc6ff8b 100644 --- a/src/mpskit_glue/transferpeps_environments.jl +++ b/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 diff --git a/test/heisenberg.jl b/test/heisenberg.jl index 44a85b43..c95ab664 100644 --- a/test/heisenberg.jl +++ b/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) diff --git a/test/tf_ising.jl b/test/tf_ising.jl index 3b3284c3..f8094391 100644 --- a/test/tf_ising.jl +++ b/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