From 52a42ece34ef814e96dc753912480063e3096cf4 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Wed, 4 Jun 2025 15:23:54 +0200 Subject: [PATCH 01/10] add horizontal correlation function --- src/PEPSKit.jl | 2 + src/algorithms/correlators.jl | 154 ++++++++++++++++++++++++++++++++++ test/examples/tf_ising.jl | 25 +++++- 3 files changed, 177 insertions(+), 4 deletions(-) create mode 100644 src/algorithms/correlators.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 4259e4b55..e1f25be75 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -63,6 +63,7 @@ include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") include("algorithms/toolbox.jl") +include("algorithms/correlators.jl") include("utility/symmetrization.jl") @@ -78,6 +79,7 @@ export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector export LocalOperator, physicalspace export expectation_value, cost_function, product_peps, correlation_length, network_value +export correlator_horizontal export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl new file mode 100644 index 000000000..b5a8d0d20 --- /dev/null +++ b/src/algorithms/correlators.jl @@ -0,0 +1,154 @@ +function correlator_horizontal( + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, + O₁::AbstractTensorMap{T,S,1,2}, + O₂::AbstractTensorMap{T,S,2,1}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}; +) where {T,S} + @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." + (r, c₁) = Tuple(inds[1]) + (r₂, c₂) = Tuple(inds[2]) + @assert r == r₂ "Not a horizontal correlation function." + @assert c₁ < c₂ "The first column index must be less than the second." + + (Nr, Nc) = size(ket) + corr = T[] + + @autoopt @tensor left_side[χS DEt Dstring DEb; χN] := + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * + env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * + ket[mod1(r, Nr), mod1(c₁, Nc)][dt; DNt DEt DSt DWt] * + conj(bra[r, c₁][db; DNb DEb DSb DWb]) * + O₁[db; dt Dstring] * + env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * + env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] + @autoopt @tensor left_side_norm[χS DEt DEb; χN] := + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * + env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * + ket[mod1(r, Nr), mod1(c₁, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c₁, Nc)][d; DNb DEb DSb DWb]) * + env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * + env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] + for c in (c₁ + 1):c₂ + final = @autoopt @tensor left_side[χ6 DWt Dstring DWb; χ1] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * + env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * + ket[mod1(r, Nr), mod1(c, Nc)][dt; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][db; DNb DEb DSb DWb]) * + O₂[Dstring db; dt] * + env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * + env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] + final_norm = @autoopt @tensor left_side_norm[χ6 DWt DWb; χ1] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * + env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * + env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] + + push!(corr, final / final_norm) + if c ≠ c₂ + @autoopt @tensor left_side[χS DEt Dstring DEb; χN] = + left_side[χ1 DWt Dstring DWb; χ4] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] + @autoopt @tensor left_side_norm[χS DEt DEb; χN] = + left_side_norm[χ1 DWt DWb; χ4] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] + end + end + return corr +end + +function correlator_horizontal( + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, + O₁::AbstractTensorMap{T,S,1,1}, + O₂::AbstractTensorMap{T,S,1,1}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, +) where {T,S} + corr = T[] + (Nr, Nc) = size(ket) + @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." + (r, c₁) = Tuple(inds[1]) + (r₂, c₂) = Tuple(inds[2]) + @assert r == r₂ "Not a horizontal correlation function." + @assert c₁ < c₂ "The first column index must be less than the second." + + @autoopt @tensor left_side[χS DEt DEb; χN] := + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * + env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * + ket[mod1(r, Nr), mod1(c₁, Nc)][dt; DNt DEt DSt DWt] * + conj(bra[r, c₁][db; DNb DEb DSb DWb]) * + O₁[db; dt] * + env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * + env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] + @autoopt @tensor left_side_norm[χS DEt DEb; χN] := + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * + env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * + ket[mod1(r, Nr), mod1(c₁, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c₁, Nc)][d; DNb DEb DSb DWb]) * + env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * + env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] + for c in (c₁ + 1):c₂ + final = @autoopt @tensor left_side[χ6 DWt DWb; χ1] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * + env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * + ket[mod1(r, Nr), mod1(c, Nc)][dt; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][db; DNb DEb DSb DWb]) * + O₂[db; dt] * + env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * + env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] + final_norm = @autoopt @tensor left_side_norm[χ6 DWt DWb; χ1] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * + env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * + env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] + push!(corr, final / final_norm) + if c ≠ c₂ + @autoopt @tensor left_side[χS DEt DEb; χN] = + left_side[χ1 DWt DWb; χ4] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] + @autoopt @tensor left_side_norm[χS DEt DEb; χN] = + left_side_norm[χ1 DWt DWb; χ4] * + env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * + ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * + conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * + env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] + end + end + return corr +end + +function correlator_horizontal( + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, + O::AbstractTensorMap{T,S,2,2}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, +) where {T,S} + U, Σ, V = tsvd(O, ((1, 3), (2, 4))) + O₁ = permute(U * sqrt(Σ), ((1,), (2, 3))) + O₂ = permute(sqrt(Σ) * V, ((1, 2), (3,))) + return correlator_horizontal(ket, bra, env, O₁, O₂, inds) +end diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index cc6758db5..2255643f1 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -33,12 +33,29 @@ peps, env, E, = fixedpoint(H, peps₀, env₀; tol=gradtol) # compute magnetization σx = TensorMap(scalartype(peps₀)[0 1; 1 0], ℂ^2, ℂ^2) -M = LocalOperator(physicalspace(H), (CartesianIndex(1, 1),) => σx) -magn = expectation_value(peps, M, env) +σz = TensorMap(scalartype(peps₀)[1 0; 0 -1], ℂ^2, ℂ^2) +Mx = LocalOperator(physicalspace(H), (CartesianIndex(1, 1),) => σx) +Mz = LocalOperator(physicalspace(H), (CartesianIndex(1, 1),) => σz) +magnx = expectation_value(peps, Mx, env) +magnz = expectation_value(peps, Mz, env) @test E ≈ e atol = 1e-2 -@test imag(magn) ≈ 0 atol = 1e-6 -@test abs(magn) ≈ mˣ atol = 5e-2 +@test imag(magnx) ≈ 0 atol = 1e-6 +@test abs(magnx) ≈ mˣ atol = 5e-2 + +# compute horizontal connected correlation function +corr = + correlator_horizontal( + peps, peps, env, σz, σz, (CartesianIndex(1, 1), CartesianIndex(1, 21)) + ) .- magnz^2 +corr_2 = + correlator_horizontal( + peps, peps, env, σz ⊗ σz, (CartesianIndex(1, 1), CartesianIndex(1, 21)) + ) .- magnz^2 + +@test corr[end] ≈ 0.0 atol = 1e-5 +@test 1 / log(corr[18] / corr[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects +@test maximum(abs.(corr - corr_2)) < 1e-14 # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5) From 4933117a0ad00df2b104ddeb7f2fee6c7f822ded Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 5 Jun 2025 17:08:15 +0200 Subject: [PATCH 02/10] use TransferMatrix and add vertical correlators --- src/PEPSKit.jl | 3 +- src/algorithms/correlators.jl | 204 ++++++++++++++----------------- src/algorithms/transfermatrix.jl | 101 +++++++++++++++ test/examples/tf_ising.jl | 13 +- 4 files changed, 204 insertions(+), 117 deletions(-) create mode 100644 src/algorithms/transfermatrix.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index e1f25be75..a110884fa 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -63,6 +63,7 @@ include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") include("algorithms/toolbox.jl") +include("algorithms/transfermatrix.jl") include("algorithms/correlators.jl") include("utility/symmetrization.jl") @@ -79,7 +80,7 @@ export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector export LocalOperator, physicalspace export expectation_value, cost_function, product_peps, correlation_length, network_value -export correlator_horizontal +export correlator_horizontal, correlator_vertical export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index b5a8d0d20..6e5e272ac 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -2,9 +2,9 @@ function correlator_horizontal( ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, - O₁::AbstractTensorMap{T,S,1,2}, - O₂::AbstractTensorMap{T,S,2,1}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}; + O₁::AbstractTensorMap{T,S}, + O₂::AbstractTensorMap{T,S}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, ) where {T,S} @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." (r, c₁) = Tuple(inds[1]) @@ -15,56 +15,60 @@ function correlator_horizontal( (Nr, Nc) = size(ket) corr = T[] - @autoopt @tensor left_side[χS DEt Dstring DEb; χN] := - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * - env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * - ket[mod1(r, Nr), mod1(c₁, Nc)][dt; DNt DEt DSt DWt] * - conj(bra[r, c₁][db; DNb DEb DSb DWb]) * - O₁[db; dt Dstring] * - env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * - env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] - @autoopt @tensor left_side_norm[χS DEt DEb; χN] := - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * - env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * - ket[mod1(r, Nr), mod1(c₁, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c₁, Nc)][d; DNb DEb DSb DWb]) * - env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * - env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] + left_start = start_left( + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)], + env.edges[4, mod1(r, Nr), _prev(c₁, Nc)], + env.corners[4, _next(r, Nr), _prev(c₁, Nc)], + ) + left_side = transfer_left( + left_start, + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], + (ket[mod1(r, Nr), mod1(c₁, Nc)], O₁, bra[mod1(r, Nr), mod1(c₁, Nc)]), + env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + ) + left_side_norm = transfer_left( + left_start, + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], + (ket[mod1(r, Nr), mod1(c₁, Nc)], bra[mod1(r, Nr), mod1(c₁, Nc)]), + env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + ) + for c in (c₁ + 1):c₂ - final = @autoopt @tensor left_side[χ6 DWt Dstring DWb; χ1] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * - env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * - ket[mod1(r, Nr), mod1(c, Nc)][dt; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][db; DNb DEb DSb DWb]) * - O₂[Dstring db; dt] * - env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * - env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] - final_norm = @autoopt @tensor left_side_norm[χ6 DWt DWb; χ1] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * - env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * - env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] + left_side_final = + left_side * HorizontalTransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], O₂, bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) + final = end_right( + left_side_final, + env.corners[2, _prev(r, Nr), _next(c, Nc)], + env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[3, _next(r, Nr), _next(c, Nc)], + ) + + left_side_norm_final = + left_side_norm * HorizontalTransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) + final_norm = end_right( + left_side_norm_final, + env.corners[2, _prev(r, Nr), _next(c, Nc)], + env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[3, _next(r, Nr), _next(c, Nc)], + ) push!(corr, final / final_norm) if c ≠ c₂ - @autoopt @tensor left_side[χS DEt Dstring DEb; χN] = - left_side[χ1 DWt Dstring DWb; χ4] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] - @autoopt @tensor left_side_norm[χS DEt DEb; χN] = - left_side_norm[χ1 DWt DWb; χ4] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] + (left_side, left_side_norm) = [ + l * HorizontalTransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) for l in (left_side, left_side_norm) + ] end end return corr @@ -74,81 +78,51 @@ function correlator_horizontal( ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, - O₁::AbstractTensorMap{T,S,1,1}, - O₂::AbstractTensorMap{T,S,1,1}, + O::AbstractTensorMap{T,S,2,2}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, +) where {T,S} + U, Σ, V = tsvd(O, ((1, 3), (2, 4))) + O₁ = permute(U * sqrt(Σ), ((1,), (2, 3))) + O₂ = permute(sqrt(Σ) * V, ((1, 2), (3,))) + return correlator_horizontal(ket, bra, env, O₁, O₂, inds) +end + +function correlator_vertical( + ket::InfinitePEPS, + bra::InfinitePEPS, + env::CTMRGEnv, + O::AbstractTensorMap{T,S,2,2}, inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, ) where {T,S} - corr = T[] - (Nr, Nc) = size(ket) - @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." (r, c₁) = Tuple(inds[1]) (r₂, c₂) = Tuple(inds[2]) - @assert r == r₂ "Not a horizontal correlation function." - @assert c₁ < c₂ "The first column index must be less than the second." - @autoopt @tensor left_side[χS DEt DEb; χN] := - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * - env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * - ket[mod1(r, Nr), mod1(c₁, Nc)][dt; DNt DEt DSt DWt] * - conj(bra[r, c₁][db; DNb DEb DSb DWb]) * - O₁[db; dt] * - env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * - env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] - @autoopt @tensor left_side_norm[χS DEt DEb; χN] := - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)][χ3; χ4] * - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)][χ4 DNt DNb; χN] * - env.edges[4, mod1(r, Nr), _prev(c₁, Nc)][χ2 DWt DWb; χ3] * - ket[mod1(r, Nr), mod1(c₁, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c₁, Nc)][d; DNb DEb DSb DWb]) * - env.corners[4, _next(r, Nr), _prev(c₁, Nc)][χ1; χ2] * - env.edges[3, _next(r, Nr), mod1(c₁, Nc)][χS DSt DSb; χ1] - for c in (c₁ + 1):c₂ - final = @autoopt @tensor left_side[χ6 DWt DWb; χ1] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * - env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * - ket[mod1(r, Nr), mod1(c, Nc)][dt; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][db; DNb DEb DSb DWb]) * - O₂[db; dt] * - env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * - env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] - final_norm = @autoopt @tensor left_side_norm[χ6 DWt DWb; χ1] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ1 DNt DNb; χ2] * - env.corners[2, _prev(r, Nr), _next(c, Nc)][χ2; χ3] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[2, mod1(r, Nr), _next(c, Nc)][χ3 DEt DEb; χ4] * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χ5 DSt DSb; χ6] * - env.corners[3, _next(r, Nr), _next(c, Nc)][χ4; χ5] - push!(corr, final / final_norm) - if c ≠ c₂ - @autoopt @tensor left_side[χS DEt DEb; χN] = - left_side[χ1 DWt DWb; χ4] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] - @autoopt @tensor left_side_norm[χS DEt DEb; χN] = - left_side_norm[χ1 DWt DWb; χ4] * - env.edges[1, _prev(r, Nr), mod1(c, Nc)][χ4 DNt DNb; χN] * - ket[mod1(r, Nr), mod1(c, Nc)][d; DNt DEt DSt DWt] * - conj(bra[mod1(r, Nr), mod1(c, Nc)][d; DNb DEb DSb DWb]) * - env.edges[3, _next(r, Nr), mod1(c, Nc)][χS DSt DSb; χ1] - end - end - return corr + return correlator_horizontal( + rotr90(ket), + rotr90(bra), + rotr90(env), + O, + (CartesianIndex(c₁, r), CartesianIndex(c₂, r₂)), + ) end -function correlator_horizontal( +function correlator_vertical( ket::InfinitePEPS, bra::InfinitePEPS, env::CTMRGEnv, - O::AbstractTensorMap{T,S,2,2}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, + O₁::AbstractTensorMap{T,S}, + O₂::AbstractTensorMap{T,S}, + inds::Tuple{CartesianIndex{2},CartesianIndex{2}}; ) where {T,S} - U, Σ, V = tsvd(O, ((1, 3), (2, 4))) - O₁ = permute(U * sqrt(Σ), ((1,), (2, 3))) - O₂ = permute(sqrt(Σ) * V, ((1, 2), (3,))) - return correlator_horizontal(ket, bra, env, O₁, O₂, inds) + (r, c₁) = Tuple(inds[1]) + (r₂, c₂) = Tuple(inds[2]) + + return correlator_horizontal( + rotr90(ket), + rotr90(bra), + rotr90(env), + O₁, + O₂, + (CartesianIndex(c₁, r), CartesianIndex(c₂, r₂)), + ) end diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl new file mode 100644 index 000000000..44b6b95b8 --- /dev/null +++ b/src/algorithms/transfermatrix.jl @@ -0,0 +1,101 @@ +abstract type AbstractTransferMatrix end; + +struct HorizontalTransferMatrix{T,S} <: AbstractTransferMatrix + above::TensorMap{T,S} + network + below::TensorMap{T,S} +end + +# TransferMatrix acting on a vector using * +Base.:*(vec, tm::AbstractTransferMatrix) = tm(vec); +# Base.:*(vec, tm::TransferMatrix) = flip(tm)(vec); + +# TransferMatrix acting as a function +function (d::HorizontalTransferMatrix)(vec) + return transfer_left(vec, d.above, d.network, d.below) +end; + +function transfer_left( + vec::TensorMap{T,S,4,1}, above, network::Tuple{PEPSTensor,PEPSTensor}, below +) where {T,S} + @autoopt @tensor vec[χS DEt Dstring DEb; χN] := + vec[χ1 DWt Dstring DWb; χ4] * + above[χ4 DNt DNb; χN] * + network[1][d; DNt DEt DSt DWt] * + conj(network[2][d; DNb DEb DSb DWb]) * + below[χS DSt DSb; χ1] + return vec +end + +function transfer_left( + vec::TensorMap{T,S,3,1}, above, network::Tuple{PEPSTensor,PEPSTensor}, below +) where {T,S} + @autoopt @tensor vec[χS DEt DEb; χN] := + vec[χ1 DWt DWb; χ4] * + above[χ4 DNt DNb; χN] * + network[1][d; DNt DEt DSt DWt] * + conj(network[2][d; DNb DEb DSb DWb]) * + below[χS DSt DSb; χ1] + return vec +end + +function transfer_left( + vec::TensorMap{T,S,3,1}, + above, + network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,1},PEPSTensor}, + below, +) where {T,S} + @autoopt @tensor vec[χS DEt DEb; χN] := + vec[χ1 DWt DWb; χ4] * + above[χ4 DNt DNb; χN] * + network[1][dt; DNt DEt DSt DWt] * + network[2][db; dt] * + conj(network[3][db; DNb DEb DSb DWb]) * + below[χS DSt DSb; χ1] + return vec +end + +function transfer_left( + vec::TensorMap{T,S,3,1}, + above, + network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,2},PEPSTensor}, + below, +) where {T,S} + @autoopt @tensor vec[χS DEt Dstring DEb; χN] := + vec[χ1 DWt DWb; χ4] * + above[χ4 DNt DNb; χN] * + network[1][dt; DNt DEt DSt DWt] * + network[2][db; dt Dstring] * + conj(network[3][db; DNb DEb DSb DWb]) * + below[χS DSt DSb; χ1] + return vec +end + +function transfer_left( + vec::TensorMap{T,S,4,1}, + above, + network::Tuple{PEPSTensor,AbstractTensorMap{T,S,2,1},PEPSTensor}, + below, +) where {T,S} + @autoopt @tensor vec[χS DEt DEb; χN] := + vec[χ1 DWt Dstring DWb; χ4] * + above[χ4 DNt DNb; χN] * + network[1][dt; DNt DEt DSt DWt] * + network[2][Dstring db; dt] * + conj(network[3][db; DNb DEb DSb DWb]) * + below[χS DSt DSb; χ1] + return vec +end + +function start_left(above, middle, below) + @autoopt @tensor vec[χ1 DWt DWb; χ4] := + above[χ3; χ4] * middle[χ2 DWt DWb; χ3] * below[χ1; χ2] + return vec +end + +function end_right(vec, above, middle, below) + return @autoopt @tensor vec[χ5 DEt DEb; χ2] * + above[χ2; χ3] * + middle[χ3 DEt DEb; χ4] * + below[χ4; χ5] +end diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index 2255643f1..b85b13fe7 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -52,10 +52,21 @@ corr_2 = correlator_horizontal( peps, peps, env, σz ⊗ σz, (CartesianIndex(1, 1), CartesianIndex(1, 21)) ) .- magnz^2 +corrv = + correlator_vertical( + peps, peps, env, σz, σz, (CartesianIndex(1, 1), CartesianIndex(21, 1)) + ) .- magnz^2 +corrv_2 = + correlator_vertical( + peps, peps, env, σz ⊗ σz, (CartesianIndex(1, 1), CartesianIndex(21, 1)) + ) .- magnz^2 @test corr[end] ≈ 0.0 atol = 1e-5 @test 1 / log(corr[18] / corr[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects -@test maximum(abs.(corr - corr_2)) < 1e-14 +@test corrv[end] ≈ 0.0 atol = 1e-5 +@test 1 / log(corrv[18] / corrv[19]) ≈ ξ_v[1] atol = 3e-2 # test correlation length far away from short-range effects +@test maximum(abs.(corrv - corrv_2)) < 1e-14 +@test maximum(abs.(corrv - corr)) < 1e-4 # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5) From 39b4b22b08b246feba239e19d7e799081f11c430 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 5 Jun 2025 18:39:51 +0200 Subject: [PATCH 03/10] change arguments --- src/algorithms/correlators.jl | 146 ++++++++++++++++++---------------- test/examples/tf_ising.jl | 61 +++++++++++--- 2 files changed, 130 insertions(+), 77 deletions(-) diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 6e5e272ac..19c979f72 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -1,16 +1,16 @@ function correlator_horizontal( - ket::InfinitePEPS, bra::InfinitePEPS, + O::Tuple{AbstractTensorMap{T,S},AbstractTensorMap{T,S}}, + i::CartesianIndex{2}, + j::CartesianIndices{2}, + ket::InfinitePEPS, env::CTMRGEnv, - O₁::AbstractTensorMap{T,S}, - O₂::AbstractTensorMap{T,S}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, ) where {T,S} @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." - (r, c₁) = Tuple(inds[1]) - (r₂, c₂) = Tuple(inds[2]) - @assert r == r₂ "Not a horizontal correlation function." - @assert c₁ < c₂ "The first column index must be less than the second." + (r, c₁) = Tuple(i) + cs = sort([ind[2] for ind in j]; dims=2) + @assert all([r == ind[1] for ind in j]) "Not a horizontal correlation function." + @assert all(c₁ .< cs) "The first column index must be less than the second." (Nr, Nc) = size(ket) corr = T[] @@ -23,7 +23,7 @@ function correlator_horizontal( left_side = transfer_left( left_start, env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], - (ket[mod1(r, Nr), mod1(c₁, Nc)], O₁, bra[mod1(r, Nr), mod1(c₁, Nc)]), + (ket[mod1(r, Nr), mod1(c₁, Nc)], O[1], bra[mod1(r, Nr), mod1(c₁, Nc)]), env.edges[3, _next(r, Nr), mod1(c₁, Nc)], ) left_side_norm = transfer_left( @@ -33,35 +33,36 @@ function correlator_horizontal( env.edges[3, _next(r, Nr), mod1(c₁, Nc)], ) - for c in (c₁ + 1):c₂ - left_side_final = - left_side * HorizontalTransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], O₂, bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], + for c in (c₁ + 1):cs[end] + if c ∈ cs + left_side_final = + left_side * HorizontalTransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], O[2], bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) + final = end_right( + left_side_final, + env.corners[2, _prev(r, Nr), _next(c, Nc)], + env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[3, _next(r, Nr), _next(c, Nc)], ) - final = end_right( - left_side_final, - env.corners[2, _prev(r, Nr), _next(c, Nc)], - env.edges[2, mod1(r, Nr), _next(c, Nc)], - env.corners[3, _next(r, Nr), _next(c, Nc)], - ) - left_side_norm_final = - left_side_norm * HorizontalTransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], + left_side_norm_final = + left_side_norm * HorizontalTransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) + final_norm = end_right( + left_side_norm_final, + env.corners[2, _prev(r, Nr), _next(c, Nc)], + env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[3, _next(r, Nr), _next(c, Nc)], ) - final_norm = end_right( - left_side_norm_final, - env.corners[2, _prev(r, Nr), _next(c, Nc)], - env.edges[2, mod1(r, Nr), _next(c, Nc)], - env.corners[3, _next(r, Nr), _next(c, Nc)], - ) - - push!(corr, final / final_norm) - if c ≠ c₂ + push!(corr, final / final_norm) + end + if c ≠ cs[end] (left_side, left_side_norm) = [ l * HorizontalTransferMatrix( env.edges[1, _prev(r, Nr), mod1(c, Nc)], @@ -75,54 +76,65 @@ function correlator_horizontal( end function correlator_horizontal( - ket::InfinitePEPS, bra::InfinitePEPS, - env::CTMRGEnv, O::AbstractTensorMap{T,S,2,2}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, + i::CartesianIndex{2}, + j::CartesianIndices{2}, + ket::InfinitePEPS, + env::CTMRGEnv, ) where {T,S} U, Σ, V = tsvd(O, ((1, 3), (2, 4))) O₁ = permute(U * sqrt(Σ), ((1,), (2, 3))) O₂ = permute(sqrt(Σ) * V, ((1, 2), (3,))) - return correlator_horizontal(ket, bra, env, O₁, O₂, inds) + return correlator_horizontal(bra, (O₁, O₂), i, j, ket, env) end function correlator_vertical( - ket::InfinitePEPS, bra::InfinitePEPS, + O, + i::CartesianIndex{2}, + j::CartesianIndices{2}, + ket::InfinitePEPS, env::CTMRGEnv, - O::AbstractTensorMap{T,S,2,2}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}, -) where {T,S} - (r, c₁) = Tuple(inds[1]) - (r₂, c₂) = Tuple(inds[2]) +) + i_rot = CartesianIndex(i[2], i[1]) + j_rot = CartesianIndex(j[1][2], j[1][1]):CartesianIndex(j[end][2], j[end][1]) - return correlator_horizontal( - rotr90(ket), - rotr90(bra), - rotr90(env), - O, - (CartesianIndex(c₁, r), CartesianIndex(c₂, r₂)), - ) + return correlator_horizontal(rotr90(bra), O, i_rot, j_rot, rotr90(ket), rotr90(env)) end -function correlator_vertical( +function correlator_horizontal( + bra::InfinitePEPS, + O, + i::CartesianIndex{2}, + j::CartesianIndex{2}, ket::InfinitePEPS, + env::CTMRGEnv, +) + return first(correlator_horizontal(bra, O, i, j:j, ket, env)) +end + +function correlator_vertical( bra::InfinitePEPS, + O, + i::CartesianIndex{2}, + j::CartesianIndex{2}, + ket::InfinitePEPS, env::CTMRGEnv, - O₁::AbstractTensorMap{T,S}, - O₂::AbstractTensorMap{T,S}, - inds::Tuple{CartesianIndex{2},CartesianIndex{2}}; -) where {T,S} - (r, c₁) = Tuple(inds[1]) - (r₂, c₂) = Tuple(inds[2]) +) + return first(correlator_vertical(bra, O, i, j:j, ket, env)) +end - return correlator_horizontal( - rotr90(ket), - rotr90(bra), - rotr90(env), - O₁, - O₂, - (CartesianIndex(c₁, r), CartesianIndex(c₂, r₂)), - ) +function correlator_horizontal( + bra::InfinitePEPS, O, i::CartesianIndex{2}, dist::Int, ket::InfinitePEPS, env::CTMRGEnv +) + j = (i + CartesianIndex(0, 1)):(i + CartesianIndex(0, dist)) + return correlator_horizontal(bra, O, i, j, ket, env) +end + +function correlator_vertical( + bra::InfinitePEPS, O, i::CartesianIndex{2}, dist::Int, ket::InfinitePEPS, env::CTMRGEnv +) + j = (i + CartesianIndex(1, 0)):(i + CartesianIndex(dist, 0)) + return correlator_vertical(bra, O, i, j, ket, env) end diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index b85b13fe7..981721845 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -43,30 +43,71 @@ magnz = expectation_value(peps, Mz, env) @test imag(magnx) ≈ 0 atol = 1e-6 @test abs(magnx) ≈ mˣ atol = 5e-2 -# compute horizontal connected correlation function -corr = +# compute connected correlation functions +corrh = correlator_horizontal( - peps, peps, env, σz, σz, (CartesianIndex(1, 1), CartesianIndex(1, 21)) + peps, + (σz, σz), + CartesianIndex(1, 1), + CartesianIndex(1, 2):CartesianIndex(1, 21), + peps, + env, ) .- magnz^2 -corr_2 = +corrh_2 = correlator_horizontal( - peps, peps, env, σz ⊗ σz, (CartesianIndex(1, 1), CartesianIndex(1, 21)) + peps, + σz ⊗ σz, + CartesianIndex(1, 1), + CartesianIndex(1, 2):CartesianIndex(1, 21), + peps, + env, ) .- magnz^2 corrv = correlator_vertical( - peps, peps, env, σz, σz, (CartesianIndex(1, 1), CartesianIndex(21, 1)) + peps, + (σz, σz), + CartesianIndex(1, 1), + CartesianIndex(2, 1):CartesianIndex(21, 1), + peps, + env, ) .- magnz^2 corrv_2 = correlator_vertical( - peps, peps, env, σz ⊗ σz, (CartesianIndex(1, 1), CartesianIndex(21, 1)) + peps, + σz ⊗ σz, + CartesianIndex(1, 1), + CartesianIndex(2, 1):CartesianIndex(21, 1), + peps, + env, ) .- magnz^2 -@test corr[end] ≈ 0.0 atol = 1e-5 -@test 1 / log(corr[18] / corr[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects + +@test corrh[end] ≈ 0.0 atol = 1e-5 +@test 1 / log(corrh[18] / corrh[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects @test corrv[end] ≈ 0.0 atol = 1e-5 @test 1 / log(corrv[18] / corrv[19]) ≈ ξ_v[1] atol = 3e-2 # test correlation length far away from short-range effects +@test maximum(abs.(corrh - corrh_2)) < 1e-14 @test maximum(abs.(corrv - corrv_2)) < 1e-14 -@test maximum(abs.(corrv - corr)) < 1e-4 +@test maximum(abs.(corrv - corrh)) < 1e-4 + +# Change from specific values and distances to a range +corrh_int = + correlator_horizontal( + peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(1, 21), peps, env + ) .- magnz^2 +corrv_int = + correlator_vertical( + peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(21, 1), peps, env + ) .- magnz^2 +corrh_dist = + correlator_horizontal(peps, (σz, σz), CartesianIndex(1, 1), 20, peps, env) .- magnz^2 +corrv_dist = + correlator_vertical(peps, (σz, σz), CartesianIndex(1, 1), 20, peps, env) .- magnz^2 + +@test corrh_int ≈ corrh[20] +@test maximum(abs.(corrh - corrh_dist)) < 1e-14 +@test corrv_int ≈ corrv[20] +@test maximum(abs.(corrv - corrv_dist)) < 1e-14 # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5) From e343f4d93e99922c58adaf6d552c5a092172d03e Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Thu, 5 Jun 2025 18:40:20 +0200 Subject: [PATCH 04/10] fix formatting --- test/examples/tf_ising.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index 981721845..6288a6b4d 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -81,7 +81,6 @@ corrv_2 = env, ) .- magnz^2 - @test corrh[end] ≈ 0.0 atol = 1e-5 @test 1 / log(corrh[18] / corrh[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects @test corrv[end] ≈ 0.0 atol = 1e-5 From ce7860d4e6a4938e6444b0bc1065840e7d2cbf0d Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 6 Jun 2025 11:20:16 +0200 Subject: [PATCH 05/10] use MPSKit.TransferMatrix change to TransferMatrix as defined in MPSKit change the argument order of transfer_left to be compatible with MPSKit add function MPSKit.correlator that checks whether the correlation function is horizontal or vertical only export correlator --- src/PEPSKit.jl | 2 +- src/algorithms/correlators.jl | 100 +++++++++++++++---------------- src/algorithms/transfermatrix.jl | 53 +++++++--------- test/examples/tf_ising.jl | 40 ++++--------- 4 files changed, 82 insertions(+), 113 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index a110884fa..66795182b 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -80,7 +80,7 @@ export CTMRGEnv, SequentialCTMRG, SimultaneousCTMRG export FixedSpaceTruncation, HalfInfiniteProjector, FullInfiniteProjector export LocalOperator, physicalspace export expectation_value, cost_function, product_peps, correlation_length, network_value -export correlator_horizontal, correlator_vertical +export correlator export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 19c979f72..c5d1457b6 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -16,61 +16,53 @@ function correlator_horizontal( corr = T[] left_start = start_left( - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)], env.edges[4, mod1(r, Nr), _prev(c₁, Nc)], + env.corners[1, _prev(r, Nr), _prev(c₁, Nc)], env.corners[4, _next(r, Nr), _prev(c₁, Nc)], ) - left_side = transfer_left( - left_start, - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], - (ket[mod1(r, Nr), mod1(c₁, Nc)], O[1], bra[mod1(r, Nr), mod1(c₁, Nc)]), - env.edges[3, _next(r, Nr), mod1(c₁, Nc)], - ) - left_side_norm = transfer_left( - left_start, - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], - (ket[mod1(r, Nr), mod1(c₁, Nc)], bra[mod1(r, Nr), mod1(c₁, Nc)]), - env.edges[3, _next(r, Nr), mod1(c₁, Nc)], - ) - + left_side = + left_start * MPSKit.TransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], + (ket[mod1(r, Nr), mod1(c₁, Nc)], O[1], bra[mod1(r, Nr), mod1(c₁, Nc)]), + env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + ) + left_side_norm = + left_start * MPSKit.TransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], + (ket[mod1(r, Nr), mod1(c₁, Nc)], bra[mod1(r, Nr), mod1(c₁, Nc)]), + env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + ) for c in (c₁ + 1):cs[end] if c ∈ cs left_side_final = - left_side * HorizontalTransferMatrix( + left_side * MPSKit.TransferMatrix( env.edges[1, _prev(r, Nr), mod1(c, Nc)], (ket[mod1(r, Nr), mod1(c, Nc)], O[2], bra[mod1(r, Nr), mod1(c, Nc)]), env.edges[3, _next(r, Nr), mod1(c, Nc)], ) final = end_right( left_side_final, - env.corners[2, _prev(r, Nr), _next(c, Nc)], env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[2, _prev(r, Nr), _next(c, Nc)], env.corners[3, _next(r, Nr), _next(c, Nc)], ) - - left_side_norm_final = - left_side_norm * HorizontalTransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], - ) + end + (left_side, left_side_norm) = [ + l * MPSKit.TransferMatrix( + env.edges[1, _prev(r, Nr), mod1(c, Nc)], + (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), + env.edges[3, _next(r, Nr), mod1(c, Nc)], + ) for l in (left_side, left_side_norm) + ] + if c ∈ cs final_norm = end_right( - left_side_norm_final, - env.corners[2, _prev(r, Nr), _next(c, Nc)], + left_side_norm, env.edges[2, mod1(r, Nr), _next(c, Nc)], + env.corners[2, _prev(r, Nr), _next(c, Nc)], env.corners[3, _next(r, Nr), _next(c, Nc)], ) push!(corr, final / final_norm) end - if c ≠ cs[end] - (left_side, left_side_norm) = [ - l * HorizontalTransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], - ) for l in (left_side, left_side_norm) - ] - end end return corr end @@ -103,18 +95,24 @@ function correlator_vertical( return correlator_horizontal(rotr90(bra), O, i_rot, j_rot, rotr90(ket), rotr90(env)) end -function correlator_horizontal( +function MPSKit.correlator( bra::InfinitePEPS, O, i::CartesianIndex{2}, - j::CartesianIndex{2}, + j::CartesianIndices{2}, ket::InfinitePEPS, env::CTMRGEnv, ) - return first(correlator_horizontal(bra, O, i, j:j, ket, env)) + if i[1] == j[1][1] + return correlator_horizontal(bra, O, i, j, ket, env) + elseif i[2] == j[1][2] + return correlator_vertical(bra, O, i, j, ket, env) + else + error("The indices must be either horizontal or vertical.") + end end -function correlator_vertical( +function MPSKit.correlator( bra::InfinitePEPS, O, i::CartesianIndex{2}, @@ -122,19 +120,21 @@ function correlator_vertical( ket::InfinitePEPS, env::CTMRGEnv, ) - return first(correlator_vertical(bra, O, i, j:j, ket, env)) -end - -function correlator_horizontal( - bra::InfinitePEPS, O, i::CartesianIndex{2}, dist::Int, ket::InfinitePEPS, env::CTMRGEnv -) - j = (i + CartesianIndex(0, 1)):(i + CartesianIndex(0, dist)) - return correlator_horizontal(bra, O, i, j, ket, env) + if i[1] == j[1] + return first(correlator_horizontal(bra, O, i, j:j, ket, env)) + elseif i[2] == j[2] + return first(correlator_vertical(bra, O, i, j:j, ket, env)) + else + error("Only horizontal and vertical correlators are implemented.") + end end -function correlator_vertical( - bra::InfinitePEPS, O, i::CartesianIndex{2}, dist::Int, ket::InfinitePEPS, env::CTMRGEnv +function MPSKit.correlator( + state::InfinitePEPS, + O, + i::CartesianIndex{2}, + j::Union{CartesianIndex{2},CartesianIndices{2}}, + env::CTMRGEnv, ) - j = (i + CartesianIndex(1, 0)):(i + CartesianIndex(dist, 0)) - return correlator_vertical(bra, O, i, j, ket, env) + return MPSKit.correlator(state, O, i, j, state, env) end diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl index 44b6b95b8..b10d07558 100644 --- a/src/algorithms/transfermatrix.jl +++ b/src/algorithms/transfermatrix.jl @@ -1,22 +1,8 @@ -abstract type AbstractTransferMatrix end; - -struct HorizontalTransferMatrix{T,S} <: AbstractTransferMatrix - above::TensorMap{T,S} - network - below::TensorMap{T,S} -end - -# TransferMatrix acting on a vector using * -Base.:*(vec, tm::AbstractTransferMatrix) = tm(vec); -# Base.:*(vec, tm::TransferMatrix) = flip(tm)(vec); - -# TransferMatrix acting as a function -function (d::HorizontalTransferMatrix)(vec) - return transfer_left(vec, d.above, d.network, d.below) -end; - -function transfer_left( - vec::TensorMap{T,S,4,1}, above, network::Tuple{PEPSTensor,PEPSTensor}, below +function MPSKit.transfer_left( + vec::TensorMap{T,S,4,1}, + network::Tuple{PEPSTensor,PEPSTensor}, + above::TensorMap{T,S,3,1}, + below::TensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt Dstring DEb; χN] := vec[χ1 DWt Dstring DWb; χ4] * @@ -27,8 +13,11 @@ function transfer_left( return vec end -function transfer_left( - vec::TensorMap{T,S,3,1}, above, network::Tuple{PEPSTensor,PEPSTensor}, below +function MPSKit.transfer_left( + vec::TensorMap{T,S,3,1}, + network::Tuple{PEPSTensor,PEPSTensor}, + above::TensorMap{T,S,3,1}, + below::TensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt DEb; χN] := vec[χ1 DWt DWb; χ4] * @@ -39,11 +28,11 @@ function transfer_left( return vec end -function transfer_left( +function MPSKit.transfer_left( vec::TensorMap{T,S,3,1}, - above, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,1},PEPSTensor}, - below, + above::TensorMap{T,S,3,1}, + below::TensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt DEb; χN] := vec[χ1 DWt DWb; χ4] * @@ -55,11 +44,11 @@ function transfer_left( return vec end -function transfer_left( +function MPSKit.transfer_left( vec::TensorMap{T,S,3,1}, - above, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,2},PEPSTensor}, - below, + above::TensorMap{T,S,3,1}, + below::TensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt Dstring DEb; χN] := vec[χ1 DWt DWb; χ4] * @@ -71,11 +60,11 @@ function transfer_left( return vec end -function transfer_left( +function MPSKit.transfer_left( vec::TensorMap{T,S,4,1}, - above, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,2,1},PEPSTensor}, - below, + above::TensorMap{T,S,3,1}, + below::TensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt DEb; χN] := vec[χ1 DWt Dstring DWb; χ4] * @@ -87,13 +76,13 @@ function transfer_left( return vec end -function start_left(above, middle, below) +function start_left(middle, above, below) @autoopt @tensor vec[χ1 DWt DWb; χ4] := above[χ3; χ4] * middle[χ2 DWt DWb; χ3] * below[χ1; χ2] return vec end -function end_right(vec, above, middle, below) +function end_right(vec, middle, above, below) return @autoopt @tensor vec[χ5 DEt DEb; χ2] * above[χ2; χ3] * middle[χ3 DEt DEb; χ4] * diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index 6288a6b4d..284e7ce25 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -45,40 +45,28 @@ magnz = expectation_value(peps, Mz, env) # compute connected correlation functions corrh = - correlator_horizontal( + correlator( peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(1, 2):CartesianIndex(1, 21), - peps, env, ) .- magnz^2 corrh_2 = - correlator_horizontal( - peps, - σz ⊗ σz, - CartesianIndex(1, 1), - CartesianIndex(1, 2):CartesianIndex(1, 21), - peps, - env, + correlator( + peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(1, 2):CartesianIndex(1, 21), env ) .- magnz^2 corrv = - correlator_vertical( + correlator( peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(2, 1):CartesianIndex(21, 1), - peps, env, ) .- magnz^2 corrv_2 = - correlator_vertical( - peps, - σz ⊗ σz, - CartesianIndex(1, 1), - CartesianIndex(2, 1):CartesianIndex(21, 1), - peps, - env, + correlator( + peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(2, 1):CartesianIndex(21, 1), env ) .- magnz^2 @test corrh[end] ≈ 0.0 atol = 1e-5 @@ -91,22 +79,14 @@ corrv_2 = # Change from specific values and distances to a range corrh_int = - correlator_horizontal( - peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(1, 21), peps, env - ) .- magnz^2 + correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(1, 21), env) - magnz^2 corrv_int = - correlator_vertical( - peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(21, 1), peps, env - ) .- magnz^2 -corrh_dist = - correlator_horizontal(peps, (σz, σz), CartesianIndex(1, 1), 20, peps, env) .- magnz^2 -corrv_dist = - correlator_vertical(peps, (σz, σz), CartesianIndex(1, 1), 20, peps, env) .- magnz^2 + correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(21, 1), env) - magnz^2 @test corrh_int ≈ corrh[20] -@test maximum(abs.(corrh - corrh_dist)) < 1e-14 @test corrv_int ≈ corrv[20] -@test maximum(abs.(corrv - corrv_dist)) < 1e-14 + +@test_broken correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(2, 2), env) # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5) From 493c7dd70ef9856a38f782e90df4bb4eb4911ad9 Mon Sep 17 00:00:00 2001 From: sanderdemeyer Date: Fri, 6 Jun 2025 15:09:53 +0200 Subject: [PATCH 06/10] use conjugated version of transfer_left --- src/algorithms/correlators.jl | 8 +++--- src/algorithms/transfermatrix.jl | 47 +++++++++++--------------------- 2 files changed, 20 insertions(+), 35 deletions(-) diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index c5d1457b6..b10c60a17 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -24,13 +24,13 @@ function correlator_horizontal( left_start * MPSKit.TransferMatrix( env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], (ket[mod1(r, Nr), mod1(c₁, Nc)], O[1], bra[mod1(r, Nr), mod1(c₁, Nc)]), - env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + _dag(env.edges[3, _next(r, Nr), mod1(c₁, Nc)]), ) left_side_norm = left_start * MPSKit.TransferMatrix( env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], (ket[mod1(r, Nr), mod1(c₁, Nc)], bra[mod1(r, Nr), mod1(c₁, Nc)]), - env.edges[3, _next(r, Nr), mod1(c₁, Nc)], + _dag(env.edges[3, _next(r, Nr), mod1(c₁, Nc)]), ) for c in (c₁ + 1):cs[end] if c ∈ cs @@ -38,7 +38,7 @@ function correlator_horizontal( left_side * MPSKit.TransferMatrix( env.edges[1, _prev(r, Nr), mod1(c, Nc)], (ket[mod1(r, Nr), mod1(c, Nc)], O[2], bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], + _dag(env.edges[3, _next(r, Nr), mod1(c, Nc)]), ) final = end_right( left_side_final, @@ -51,7 +51,7 @@ function correlator_horizontal( l * MPSKit.TransferMatrix( env.edges[1, _prev(r, Nr), mod1(c, Nc)], (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), - env.edges[3, _next(r, Nr), mod1(c, Nc)], + _dag(env.edges[3, _next(r, Nr), mod1(c, Nc)]), ) for l in (left_side, left_side_norm) ] if c ∈ cs diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl index b10d07558..3b31e6c22 100644 --- a/src/algorithms/transfermatrix.jl +++ b/src/algorithms/transfermatrix.jl @@ -1,38 +1,23 @@ function MPSKit.transfer_left( - vec::TensorMap{T,S,4,1}, + vec::AbstractTensorMap{T,S,4,1}, network::Tuple{PEPSTensor,PEPSTensor}, - above::TensorMap{T,S,3,1}, - below::TensorMap{T,S,3,1}, + above::AbstractTensorMap{T,S,3,1}, + below::AbstractTensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt Dstring DEb; χN] := vec[χ1 DWt Dstring DWb; χ4] * above[χ4 DNt DNb; χN] * network[1][d; DNt DEt DSt DWt] * conj(network[2][d; DNb DEb DSb DWb]) * - below[χS DSt DSb; χ1] + conj(below[χ1 DSt DSb; χS]) return vec end function MPSKit.transfer_left( - vec::TensorMap{T,S,3,1}, - network::Tuple{PEPSTensor,PEPSTensor}, - above::TensorMap{T,S,3,1}, - below::TensorMap{T,S,3,1}, -) where {T,S} - @autoopt @tensor vec[χS DEt DEb; χN] := - vec[χ1 DWt DWb; χ4] * - above[χ4 DNt DNb; χN] * - network[1][d; DNt DEt DSt DWt] * - conj(network[2][d; DNb DEb DSb DWb]) * - below[χS DSt DSb; χ1] - return vec -end - -function MPSKit.transfer_left( - vec::TensorMap{T,S,3,1}, + vec::AbstractTensorMap{T,S,3,1}, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,1},PEPSTensor}, - above::TensorMap{T,S,3,1}, - below::TensorMap{T,S,3,1}, + above::AbstractTensorMap{T,S,3,1}, + below::AbstractTensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt DEb; χN] := vec[χ1 DWt DWb; χ4] * @@ -40,15 +25,15 @@ function MPSKit.transfer_left( network[1][dt; DNt DEt DSt DWt] * network[2][db; dt] * conj(network[3][db; DNb DEb DSb DWb]) * - below[χS DSt DSb; χ1] + conj(below[χ1 DSt DSb; χS]) return vec end function MPSKit.transfer_left( - vec::TensorMap{T,S,3,1}, + vec::AbstractTensorMap{T,S,3,1}, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,2},PEPSTensor}, - above::TensorMap{T,S,3,1}, - below::TensorMap{T,S,3,1}, + above::AbstractTensorMap{T,S,3,1}, + below::AbstractTensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt Dstring DEb; χN] := vec[χ1 DWt DWb; χ4] * @@ -56,15 +41,15 @@ function MPSKit.transfer_left( network[1][dt; DNt DEt DSt DWt] * network[2][db; dt Dstring] * conj(network[3][db; DNb DEb DSb DWb]) * - below[χS DSt DSb; χ1] + conj(below[χ1 DSt DSb; χS]) return vec end function MPSKit.transfer_left( - vec::TensorMap{T,S,4,1}, + vec::AbstractTensorMap{T,S,4,1}, network::Tuple{PEPSTensor,AbstractTensorMap{T,S,2,1},PEPSTensor}, - above::TensorMap{T,S,3,1}, - below::TensorMap{T,S,3,1}, + above::AbstractTensorMap{T,S,3,1}, + below::AbstractTensorMap{T,S,3,1}, ) where {T,S} @autoopt @tensor vec[χS DEt DEb; χN] := vec[χ1 DWt Dstring DWb; χ4] * @@ -72,7 +57,7 @@ function MPSKit.transfer_left( network[1][dt; DNt DEt DSt DWt] * network[2][Dstring db; dt] * conj(network[3][db; DNb DEb DSb DWb]) * - below[χS DSt DSb; χ1] + conj(below[χ1 DSt DSb; χS]) return vec end From 1f550bbf56726c08008bf1e0ae10fd916636d70b Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 12 Jun 2025 10:23:03 -0400 Subject: [PATCH 07/10] Add excited transfer --- src/algorithms/contractions/vumps_contractions.jl | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 9510dcbe2..019cf2a52 100644 --- a/src/algorithms/contractions/vumps_contractions.jl +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -40,6 +40,21 @@ function _transfer_left( A[χ_NW D_N_above D_N_below; χ_NE] end +# transfer with excited GL +function MPSKit.transfer_left( + GL::GenericMPSTensor{S,4}, + O::PEPSSandwich, + A::GenericMPSTensor{S,3}, + Ā::GenericMPSTensor{S,3}, +) where {S} + @autoopt @tensor GL′[χ_SE D_E_above d_string D_E_below; χ_NE] := + GL[χ_SW D_W_above d_string D_W_below; χ_NW] * + conj(Ā[χ_SW D_S_above D_S_below; χ_SE]) * + ket(O)[d; D_N_above D_E_above D_S_above D_W_above] * + conj(bra(O)[d; D_N_below D_E_below D_S_below D_W_below]) * + A[χ_NW D_N_above D_N_below; χ_NE] +end + function _transfer_right( GR::GenericMPSTensor{S,3}, O::PEPSSandwich, From 15ee9133ce8a859ee7371c6b9c08a819a2443196 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 12 Jun 2025 10:24:06 -0400 Subject: [PATCH 08/10] Rewrite correlators --- src/PEPSKit.jl | 7 +- .../contractions/vumps_contractions.jl | 2 - src/algorithms/correlators.jl | 237 +++++++++++------- src/algorithms/transfermatrix.jl | 75 ------ 4 files changed, 156 insertions(+), 165 deletions(-) delete mode 100644 src/algorithms/transfermatrix.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 162a49fe0..391a7e151 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -4,10 +4,14 @@ using LinearAlgebra, Statistics, Base.Threads, Base.Iterators, Printf using Compat using Accessors: @set, @reset using VectorInterface -using TensorKit, KrylovKit, MPSKit, OptimKit, TensorOperations +using TensorKit, KrylovKit, OptimKit, TensorOperations using ChainRulesCore, Zygote using LoggingExtras + +using MPSKit +using MPSKit: MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix import MPSKit: leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace + using MPSKitModels using FiniteDifferences using OhMyThreads: tmap @@ -65,7 +69,6 @@ include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") include("algorithms/toolbox.jl") -include("algorithms/transfermatrix.jl") include("algorithms/correlators.jl") include("utility/symmetrization.jl") diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 019cf2a52..c886da6ef 100644 --- a/src/algorithms/contractions/vumps_contractions.jl +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -1,5 +1,3 @@ -using MPSKit: GenericMPSTensor, MPSBondTensor - # # Environment transfer functions # diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index b10c60a17..75f5a0f4d 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -1,114 +1,185 @@ function correlator_horizontal( bra::InfinitePEPS, - O::Tuple{AbstractTensorMap{T,S},AbstractTensorMap{T,S}}, + operator, i::CartesianIndex{2}, - j::CartesianIndices{2}, + js::AbstractVector{CartesianIndex{2}}, ket::InfinitePEPS, env::CTMRGEnv, -) where {T,S} - @assert size(ket) == size(bra) "The ket and bra must have the same unit cell." - (r, c₁) = Tuple(i) - cs = sort([ind[2] for ind in j]; dims=2) - @assert all([r == ind[1] for ind in j]) "Not a horizontal correlation function." - @assert all(c₁ .< cs) "The first column index must be less than the second." - - (Nr, Nc) = size(ket) - corr = T[] - - left_start = start_left( - env.edges[4, mod1(r, Nr), _prev(c₁, Nc)], - env.corners[1, _prev(r, Nr), _prev(c₁, Nc)], - env.corners[4, _next(r, Nr), _prev(c₁, Nc)], +) + size(ket) == size(bra) || + throw(DimensionMismatch("The ket and bra must have the same unit cell.")) + all(==(i[1]) ∘ first ∘ Tuple, js) || + throw(ArgumentError("Not a horizontal correlation function")) + issorted(vcat(i, js); by=last ∘ Tuple) || + throw(ArgumentError("Not an increasing sequence of coordinates")) + + O = FiniteMPO(operator) + length(O) == 2 || throw(ArgumentError("Operator must act on two sites")) + + # preallocate with correct scalartype + G = similar( + js, + TensorOperations.promote_contract( + scalartype(bra), scalartype(ket), scalartype(env), scalartype.(O)... + ), ) - left_side = - left_start * MPSKit.TransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], - (ket[mod1(r, Nr), mod1(c₁, Nc)], O[1], bra[mod1(r, Nr), mod1(c₁, Nc)]), - _dag(env.edges[3, _next(r, Nr), mod1(c₁, Nc)]), - ) - left_side_norm = - left_start * MPSKit.TransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c₁, Nc)], - (ket[mod1(r, Nr), mod1(c₁, Nc)], bra[mod1(r, Nr), mod1(c₁, Nc)]), - _dag(env.edges[3, _next(r, Nr), mod1(c₁, Nc)]), - ) - for c in (c₁ + 1):cs[end] - if c ∈ cs - left_side_final = - left_side * MPSKit.TransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], O[2], bra[mod1(r, Nr), mod1(c, Nc)]), - _dag(env.edges[3, _next(r, Nr), mod1(c, Nc)]), - ) - final = end_right( - left_side_final, - env.edges[2, mod1(r, Nr), _next(c, Nc)], - env.corners[2, _prev(r, Nr), _next(c, Nc)], - env.corners[3, _next(r, Nr), _next(c, Nc)], - ) - end - (left_side, left_side_norm) = [ - l * MPSKit.TransferMatrix( - env.edges[1, _prev(r, Nr), mod1(c, Nc)], - (ket[mod1(r, Nr), mod1(c, Nc)], bra[mod1(r, Nr), mod1(c, Nc)]), - _dag(env.edges[3, _next(r, Nr), mod1(c, Nc)]), - ) for l in (left_side, left_side_norm) - ] - if c ∈ cs - final_norm = end_right( - left_side_norm, - env.edges[2, mod1(r, Nr), _next(c, Nc)], - env.corners[2, _prev(r, Nr), _next(c, Nc)], - env.corners[3, _next(r, Nr), _next(c, Nc)], + + # left start for operator and norm contractions + Vn, Vo = start_correlator(i, bra, O[1], ket, env) + i += CartesianIndex(0, 1) + + for (k, j) in enumerate(js) + # transfer until left of site j + while j > i + Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + sandwich = ( + ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)] ) - push!(corr, final / final_norm) + T = TransferMatrix(Atop, sandwich, _dag(Abot)) + Vo = Vo * T + Vn = Vn * T + i += CartesianIndex(0, 1) end + + # compute overlap with operator + numerator = end_correlator_numerator(j, Vo, bra, O[2], ket, env) + + # transfer right of site j + Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + sandwich = ( + ket[mod1(i[1], end), mod1(i[2], end)], bra[mod1(i[1], end), mod1(i[2], end)] + ) + T = TransferMatrix(Atop, sandwich, _dag(Abot)) + Vo = Vo * T + Vn = Vn * T + i += CartesianIndex(0, 1) + + # compute overlap without operator + denominator = end_correlator_denominator(j, Vn, env) + + G[k] = numerator / denominator end - return corr + + return G end -function correlator_horizontal( - bra::InfinitePEPS, - O::AbstractTensorMap{T,S,2,2}, +function start_correlator( i::CartesianIndex{2}, - j::CartesianIndices{2}, - ket::InfinitePEPS, + below::InfinitePEPS, + O::MPOTensor, + above::InfinitePEPS, env::CTMRGEnv, -) where {T,S} - U, Σ, V = tsvd(O, ((1, 3), (2, 4))) - O₁ = permute(U * sqrt(Σ), ((1,), (2, 3))) - O₂ = permute(sqrt(Σ) * V, ((1, 2), (3,))) - return correlator_horizontal(bra, (O₁, O₂), i, j, ket, env) +) + r, c = Tuple(i) + E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] + E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] + E_west = env.edges[WEST, mod1(r, end), _prev(c, end)] + C_northwest = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] + C_southwest = env.corners[SOUTHWEST, _next(r, end), _prev(c, end)] + sandwich = (below[mod1(r, end), mod1(c, end)], above[mod1(r, end), mod1(c, end)]) + + # TODO: part of these contractions is duplicated between the two output tensors, + # so could be optimized further + @autoopt @tensor Vn[χSE Detop Debot; χNE] := + E_south[χSE Dstop Dsbot; χSW2] * + C_southwest[χSW2; χSW] * + E_west[χSW Dwtop Dwbot; χNW] * + C_northwest[χNW; χN] * + conj(bra(sandwich)[d; Dnbot Debot Dsbot Dwbot]) * + ket(sandwich)[d; Dntop Detop Dstop Dwtop] * + E_north[χN Dntop Dnbot; χNE] + + @autoopt @tensor Vo[χSE Detop dstring Debot; χNE] := + E_south[χSE Dstop Dsbot; χSW2] * + C_southwest[χSW2; χSW] * + E_west[χSW Dwtop Dwbot; χNW] * + C_northwest[χNW; χN] * + conj(bra(sandwich)[d1; Dnbot Debot Dsbot Dwbot]) * + removeunit(O, 1)[d1; d2 dstring] * + ket(sandwich)[d2; Dntop Detop Dstop Dwtop] * + E_north[χN Dntop Dnbot; χNE] + + return Vn, Vo +end + +function end_correlator_numerator( + j::CartesianIndex{2}, + V, + above::InfinitePEPS, + O::MPOTensor, + below::InfinitePEPS, + env::CTMRGEnv, +) + r, c = Tuple(j) + E_north = env.edges[NORTH, _prev(r, end), mod1(c, end)] + E_east = env.edges[EAST, mod1(r, end), _next(c, end)] + E_south = env.edges[SOUTH, _next(r, end), mod1(c, end)] + C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] + C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] + sandwich = (above[mod1(r, end), mod1(c, end)], below[mod1(r, end), mod1(c, end)]) + + return @autoopt @tensor contractcheck = true V[χSW DWt dstring DWb; χNW] * + E_south[χSSE DSt DSb; χSW] * + E_east[χNEE DEt DEb; χSEE] * + E_north[χNW DNt DNb; χNNE] * + C_northeast[χNNE; χNEE] * + C_southeast[χSEE; χSSE] * + conj(bra(sandwich)[db; DNb DEb DSb DWb]) * + ket(sandwich)[dt; DNt DEt DSt DWt] * + removeunit(O, 4)[dstring db; dt] +end + +function end_correlator_denominator(j::CartesianIndex{2}, V, env::CTMRGEnv) + r, c = Tuple(j) + C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] + E_east = env.edges[EAST, mod1(r, end), _next(c, end)] + C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] + + return @autoopt @tensor V[χS DEt DEb; χN] * + C_northeast[χN; χNE] * + E_east[χNE DEt DEb; χSE] * + C_southeast[χSE; χS] end function correlator_vertical( bra::InfinitePEPS, O, i::CartesianIndex{2}, - j::CartesianIndices{2}, + js::AbstractVector{CartesianIndex{2}}, ket::InfinitePEPS, env::CTMRGEnv, ) - i_rot = CartesianIndex(i[2], i[1]) - j_rot = CartesianIndex(j[1][2], j[1][1]):CartesianIndex(j[end][2], j[end][1]) + rotated_bra = rotl90(bra) + rotated_ket = bra === ket ? rotated_bra : rotl90(ket) - return correlator_horizontal(rotr90(bra), O, i_rot, j_rot, rotr90(ket), rotr90(env)) + rotated_i = rotl90(i) + rotated_j = map(rotl90, js) + + return correlator_horizontal( + rotated_bra, O, rotated_i, rotated_j, rotated_ket, rotl90(env) + ) end +const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}},CartesianIndices{N}} + function MPSKit.correlator( bra::InfinitePEPS, O, i::CartesianIndex{2}, - j::CartesianIndices{2}, + js::CoordCollection{2}, ket::InfinitePEPS, env::CTMRGEnv, ) - if i[1] == j[1][1] - return correlator_horizontal(bra, O, i, j, ket, env) - elseif i[2] == j[1][2] - return correlator_vertical(bra, O, i, j, ket, env) + js = vec(js) # map CartesianIndices to actual Vector instead of Matrix + + if all(==(i[1]) ∘ first ∘ Tuple, js) + return correlator_horizontal(bra, O, i, js, ket, env) + elseif all(==(i[2]) ∘ last ∘ Tuple, js) + return correlator_vertical(bra, O, i, js, ket, env) else - error("The indices must be either horizontal or vertical.") + error("Only horizontal or vertical correlators are implemented") end end @@ -120,20 +191,14 @@ function MPSKit.correlator( ket::InfinitePEPS, env::CTMRGEnv, ) - if i[1] == j[1] - return first(correlator_horizontal(bra, O, i, j:j, ket, env)) - elseif i[2] == j[2] - return first(correlator_vertical(bra, O, i, j:j, ket, env)) - else - error("Only horizontal and vertical correlators are implemented.") - end + return only(correlator(bra, O, i, j:j, ket, env)) end function MPSKit.correlator( state::InfinitePEPS, O, i::CartesianIndex{2}, - j::Union{CartesianIndex{2},CartesianIndices{2}}, + j, env::CTMRGEnv, ) return MPSKit.correlator(state, O, i, j, state, env) diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl deleted file mode 100644 index 3b31e6c22..000000000 --- a/src/algorithms/transfermatrix.jl +++ /dev/null @@ -1,75 +0,0 @@ -function MPSKit.transfer_left( - vec::AbstractTensorMap{T,S,4,1}, - network::Tuple{PEPSTensor,PEPSTensor}, - above::AbstractTensorMap{T,S,3,1}, - below::AbstractTensorMap{T,S,3,1}, -) where {T,S} - @autoopt @tensor vec[χS DEt Dstring DEb; χN] := - vec[χ1 DWt Dstring DWb; χ4] * - above[χ4 DNt DNb; χN] * - network[1][d; DNt DEt DSt DWt] * - conj(network[2][d; DNb DEb DSb DWb]) * - conj(below[χ1 DSt DSb; χS]) - return vec -end - -function MPSKit.transfer_left( - vec::AbstractTensorMap{T,S,3,1}, - network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,1},PEPSTensor}, - above::AbstractTensorMap{T,S,3,1}, - below::AbstractTensorMap{T,S,3,1}, -) where {T,S} - @autoopt @tensor vec[χS DEt DEb; χN] := - vec[χ1 DWt DWb; χ4] * - above[χ4 DNt DNb; χN] * - network[1][dt; DNt DEt DSt DWt] * - network[2][db; dt] * - conj(network[3][db; DNb DEb DSb DWb]) * - conj(below[χ1 DSt DSb; χS]) - return vec -end - -function MPSKit.transfer_left( - vec::AbstractTensorMap{T,S,3,1}, - network::Tuple{PEPSTensor,AbstractTensorMap{T,S,1,2},PEPSTensor}, - above::AbstractTensorMap{T,S,3,1}, - below::AbstractTensorMap{T,S,3,1}, -) where {T,S} - @autoopt @tensor vec[χS DEt Dstring DEb; χN] := - vec[χ1 DWt DWb; χ4] * - above[χ4 DNt DNb; χN] * - network[1][dt; DNt DEt DSt DWt] * - network[2][db; dt Dstring] * - conj(network[3][db; DNb DEb DSb DWb]) * - conj(below[χ1 DSt DSb; χS]) - return vec -end - -function MPSKit.transfer_left( - vec::AbstractTensorMap{T,S,4,1}, - network::Tuple{PEPSTensor,AbstractTensorMap{T,S,2,1},PEPSTensor}, - above::AbstractTensorMap{T,S,3,1}, - below::AbstractTensorMap{T,S,3,1}, -) where {T,S} - @autoopt @tensor vec[χS DEt DEb; χN] := - vec[χ1 DWt Dstring DWb; χ4] * - above[χ4 DNt DNb; χN] * - network[1][dt; DNt DEt DSt DWt] * - network[2][Dstring db; dt] * - conj(network[3][db; DNb DEb DSb DWb]) * - conj(below[χ1 DSt DSb; χS]) - return vec -end - -function start_left(middle, above, below) - @autoopt @tensor vec[χ1 DWt DWb; χ4] := - above[χ3; χ4] * middle[χ2 DWt DWb; χ3] * below[χ1; χ2] - return vec -end - -function end_right(vec, middle, above, below) - return @autoopt @tensor vec[χ5 DEt DEb; χ2] * - above[χ2; χ3] * - middle[χ3 DEt DEb; χ4] * - below[χ4; χ5] -end From 2d69d80de60ef3cc61734588e1f7133ad3281ce9 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 12 Jun 2025 10:24:12 -0400 Subject: [PATCH 09/10] Update tests --- test/examples/tf_ising.jl | 44 ++++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index 284e7ce25..c0229cb7d 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -45,48 +45,54 @@ magnz = expectation_value(peps, Mz, env) # compute connected correlation functions corrh = - correlator( - peps, - (σz, σz), - CartesianIndex(1, 1), - CartesianIndex(1, 2):CartesianIndex(1, 21), - env, - ) .- magnz^2 -corrh_2 = correlator( peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(1, 2):CartesianIndex(1, 21), env ) .- magnz^2 corrv = correlator( peps, - (σz, σz), + σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(2, 1):CartesianIndex(21, 1), env, ) .- magnz^2 -corrv_2 = - correlator( - peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(2, 1):CartesianIndex(21, 1), env - ) .- magnz^2 @test corrh[end] ≈ 0.0 atol = 1e-5 @test 1 / log(corrh[18] / corrh[19]) ≈ ξ_h[1] atol = 2e-2 # test correlation length far away from short-range effects @test corrv[end] ≈ 0.0 atol = 1e-5 @test 1 / log(corrv[18] / corrv[19]) ≈ ξ_v[1] atol = 3e-2 # test correlation length far away from short-range effects -@test maximum(abs.(corrh - corrh_2)) < 1e-14 -@test maximum(abs.(corrv - corrv_2)) < 1e-14 -@test maximum(abs.(corrv - corrh)) < 1e-4 +@test corrv ≈ corrh rtol = 1e-2 + +# compute weird geometries +corrh_2 = + correlator( + peps, + σz ⊗ σz, + CartesianIndex(3, 2), + CartesianIndex(3, 3):CartesianIndex(1, 2):CartesianIndex(3, 7), + env, + ) .- magnz^2 +@test corrh_2 ≈ corrh[1:2:5] +corrv_2 = + correlator( + peps, + σz ⊗ σz, + CartesianIndex(2, 3), + CartesianIndex(3, 3):CartesianIndex(2, 1):CartesianIndex(7, 3), + env, + ) .- magnz^2 +@test corrv_2 ≈ corrv[1:2:5] # Change from specific values and distances to a range corrh_int = - correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(1, 21), env) - magnz^2 + correlator(peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(1, 21), env) - magnz^2 corrv_int = - correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(21, 1), env) - magnz^2 + correlator(peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(21, 1), env) - magnz^2 @test corrh_int ≈ corrh[20] @test corrv_int ≈ corrv[20] -@test_broken correlator(peps, (σz, σz), CartesianIndex(1, 1), CartesianIndex(2, 2), env) +@test_broken correlator(peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(2, 2), env) # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5) From 84d659d3bfb21de9100ee9d676d66706d135e08d Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 12 Jun 2025 10:26:23 -0400 Subject: [PATCH 10/10] Formatter --- src/algorithms/correlators.jl | 8 +------- test/examples/tf_ising.jl | 6 +----- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 75f5a0f4d..9b5ce5a0e 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -194,12 +194,6 @@ function MPSKit.correlator( return only(correlator(bra, O, i, j:j, ket, env)) end -function MPSKit.correlator( - state::InfinitePEPS, - O, - i::CartesianIndex{2}, - j, - env::CTMRGEnv, -) +function MPSKit.correlator(state::InfinitePEPS, O, i::CartesianIndex{2}, j, env::CTMRGEnv) return MPSKit.correlator(state, O, i, j, state, env) end diff --git a/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index c0229cb7d..9ef71b62c 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -50,11 +50,7 @@ corrh = ) .- magnz^2 corrv = correlator( - peps, - σz ⊗ σz, - CartesianIndex(1, 1), - CartesianIndex(2, 1):CartesianIndex(21, 1), - env, + peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(2, 1):CartesianIndex(21, 1), env ) .- magnz^2 @test corrh[end] ≈ 0.0 atol = 1e-5