From 79c88ab30e54300147dd16244a8488cf303a30e5 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 11 Oct 2025 09:54:33 +0800 Subject: [PATCH 01/13] Refactor iPEPS correlator --- src/PEPSKit.jl | 3 +- .../contractions/correlator/pepo_1layer.jl | 147 ++++++++++++++++++ .../correlator/peps.jl} | 15 +- src/operators/infinitepepo.jl | 9 +- test/utility/correlator.jl | 62 +++++--- 5 files changed, 207 insertions(+), 29 deletions(-) create mode 100644 src/algorithms/contractions/correlator/pepo_1layer.jl rename src/algorithms/{correlators.jl => contractions/correlator/peps.jl} (94%) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 9e1c24b74..b279315d6 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -57,6 +57,8 @@ include("algorithms/contractions/bondenv/benv_tools.jl") include("algorithms/contractions/bondenv/gaugefix.jl") include("algorithms/contractions/bondenv/als_solve.jl") include("algorithms/contractions/bondenv/benv_ctm.jl") +include("algorithms/contractions/correlator/peps.jl") +include("algorithms/contractions/correlator/pepo_1layer.jl") include("algorithms/ctmrg/sparse_environments.jl") include("algorithms/ctmrg/ctmrg.jl") @@ -74,7 +76,6 @@ include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") include("algorithms/toolbox.jl") -include("algorithms/correlators.jl") include("utility/symmetrization.jl") diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl new file mode 100644 index 000000000..8cdff8519 --- /dev/null +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -0,0 +1,147 @@ +function correlator_horizontal( + ρ::InfinitePEPO, operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + env::CTMRGEnv + ) + (size(ρ, 3) == 1) || + throw(ArgumentError("The input PEPO ρ must have only one layer.")) + 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(ρ), scalartype(env), scalartype.(O)... + ), + ) + # left start for operator and norm contractions + Vn, Vo = start_correlator(i, ρ, O[1], 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)] + t = ρ[mod1(i[1], end), mod1(i[2], end)] + @tensor Amid[w s; n e] := t[d d; n e s w] + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = TransferMatrix(Atop, Amid, _dag(Abot)) + twistdual!(T.below, 2:numout(T.below)) + Vn = Vn * T + Vo = Vo * T + i += CartesianIndex(0, 1) + end + # compute overlap with operator + numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) + # transfer right of site j + Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + t = ρ[mod1(i[1], end), mod1(i[2], end)] + @tensor Amid[w s; n e] := t[d d; n e s w] + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = TransferMatrix(Atop, Amid, _dag(Abot)) + twistdual!(T.below, 2:numout(T.below)) + Vn = Vn * T + i += CartesianIndex(0, 1) + # compute overlap without operator + denominator = end_correlator_denominator(j, Vn, env) + G[k] = numerator / denominator + end + return G +end + +function start_correlator( + i::CartesianIndex{2}, ρ::InfinitePEPO, + O::MPOTensor, env::CTMRGEnv + ) + (size(ρ, 3) == 1) || + throw(ArgumentError("The input PEPO ρ must have only one layer.")) + 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)] + t = ρ[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 De; χNE] := + E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * + E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * + t[d d; Dn De Ds Dw] * E_north[χN Dn; χNE] + @autoopt @tensor Vo[χSE dstring De; χNE] := + E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * + E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * + removeunit(O, 1)[d1; d2 dstring] * + t[d2 d1; Dn De Ds Dw] * E_north[χN Dn; χNE] + return Vn, Vo +end + +function end_correlator_numerator( + j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1}, + ρ::InfinitePEPO, + O::MPOTensor, env::CTMRGEnv + ) where {T, S} + (size(ρ, 3) == 1) || + throw(ArgumentError("The input PEPO ρ must have only one layer.")) + 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)] + t = ρ[mod1(r, end), mod1(c, end)] + return @autoopt @tensor V[χSW dstring DW; χNW] * + E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * + C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * + t[dt db; DN DE DS DW] * removeunit(O, 4)[dstring db; dt] +end + +function end_correlator_denominator( + j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 2, 1}, env::CTMRGEnv + ) where {T, S} + 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 DE; χN] * C_northeast[χN; χNE] * + E_east[χNE DE; χSE] * C_southeast[χSE; χS] +end + +function correlator_vertical( + ρ::InfinitePEPO, operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + env::CTMRGEnv, + ) + rotated_ρ = rotl90(ρ) + rotated_i = siterotl90(i, size(bra)) + rotated_js = map(j -> siterotl90(j, size(bra)), js) + return correlator_horizontal( + rotated_ρ, operator, rotated_i, rotated_js, rotl90(env) + ) +end + +function MPSKit.correlator( + ρ::InfinitePEPO, O, + i::CartesianIndex{2}, js::CoordCollection{2}, + env::CTMRGEnv, + ) + js = vec(js) # map CartesianIndices to Vector instead of Matrix + if all(==(i[1]) ∘ first ∘ Tuple, js) + return correlator_horizontal(ρ, O, i, js, env) + elseif all(==(i[2]) ∘ last ∘ Tuple, js) + return correlator_vertical(ρ, O, i, js, env) + else + error("Only horizontal or vertical correlators are implemented") + end +end + +function MPSKit.correlator( + ρ::InfinitePEPO, O, + i::CartesianIndex{2}, j::CartesianIndex{2}, + env::CTMRGEnv, + ) + return only(correlator(ρ, O, i, j:j, env)) +end diff --git a/src/algorithms/correlators.jl b/src/algorithms/contractions/correlator/peps.jl similarity index 94% rename from src/algorithms/correlators.jl rename to src/algorithms/contractions/correlator/peps.jl index f97d03b84..f4c7bc93a 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -107,12 +107,12 @@ end function end_correlator_numerator( j::CartesianIndex{2}, - V, + V::AbstractTensorMap{T, S, 4, 1}, above::InfinitePEPS, O::MPOTensor, below::InfinitePEPS, env::CTMRGEnv, - ) + ) where {T, S} 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)] @@ -132,7 +132,10 @@ function end_correlator_numerator( removeunit(O, 4)[dstring db; dt] end -function end_correlator_denominator(j::CartesianIndex{2}, V, env::CTMRGEnv) +function end_correlator_denominator( + j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1}, + env::CTMRGEnv + ) where {T, S} 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)] @@ -146,7 +149,7 @@ end function correlator_vertical( bra::InfinitePEPS, - O, + operator, i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, ket::InfinitePEPS, env::CTMRGEnv, @@ -155,10 +158,10 @@ function correlator_vertical( rotated_ket = bra === ket ? rotated_bra : rotl90(ket) rotated_i = siterotl90(i, size(bra)) - rotated_j = map(j -> siterotl90(j, size(bra)), js) + rotated_js = map(j -> siterotl90(j, size(bra)), js) return correlator_horizontal( - rotated_bra, O, rotated_i, rotated_j, rotated_ket, rotl90(env) + rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env) ) end diff --git a/src/operators/infinitepepo.jl b/src/operators/infinitepepo.jl index 92e28ae71..80043f5a4 100644 --- a/src/operators/infinitepepo.jl +++ b/src/operators/infinitepepo.jl @@ -65,14 +65,17 @@ end function InfinitePEPO( Pspaces::A, Nspaces::A, Espaces::A = Nspaces ) where {A <: AbstractArray{<:ElementarySpace, 2}} + return InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) +end +function InfinitePEPO( + f, T, Pspaces::A, Nspaces::A, Espaces::A = Nspaces + ) where {A <: AbstractArray{<:ElementarySpace, 2}} size(Pspaces) == size(Nspaces) == size(Espaces) || throw(ArgumentError("Input spaces should have equal sizes.")) - Pspaces = reshape(Pspaces, (size(Pspaces)..., 1)) Nspaces = reshape(Pspaces, (size(Nspaces)..., 1)) Espaces = reshape(Pspaces, (size(Espaces)..., 1)) - - return InfinitePEPO(Pspaces, Nspaces, Espaces) + return InfinitePEPO(f, T, Pspaces, Nspaces, Espaces) end """ diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl index 3ba9786ee..7658901b4 100644 --- a/test/utility/correlator.jl +++ b/test/utility/correlator.jl @@ -2,27 +2,51 @@ using Test using Random using TensorKit using PEPSKit -import MPSKitModels: TJOperators as tJ -Pspace = tJ.tj_space(Trivial, Trivial) -Vspace = Vect[FermionParity](0 => 2, 1 => 2) -Espace = Vect[FermionParity](0 => 3, 1 => 3) -Random.seed!(100) -peps = InfinitePEPS(rand, ComplexF64, Pspace, Vspace; unitcell = (2, 2)); -env = CTMRGEnv(rand, ComplexF64, peps, Espace); -lattice = collect(space(t, 1) for t in peps.A) +Nr, Nc = 2, 2 +Vphy = Vect[FermionParity](0 => 1, 1 => 1) +Pspaces = fill(Vphy, (Nr, Nc)) +op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) + +V = Vect[FermionParity](0 => 1, 1 => 2) +Venv = Vect[FermionParity](0 => 2, 1 => 2) +Nspaces = [V' V; V V'] +Espaces = [V V'; V' V] site0 = CartesianIndex(1, 1) -maxsep = 8 -site1s = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep) +maxsep = 6 +site1xs = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep) +site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep) -op = tJ.S_exchange(ComplexF64, Trivial, Trivial); +@testset "Correlator in InfinitePEPS" begin + Random.seed!(100) + peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) + env = CTMRGEnv(randn, ComplexF64, peps, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(peps, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(peps, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 + end +end -vals1 = correlator(peps, op, site0, site1s, env) -vals2 = collect( - begin - O = LocalOperator(lattice, (site0, site1) => op) - val = expectation_value(peps, O, env) - end for site1 in site1s -) -@test vals1 ≈ vals2 +@testset "Correlator in InfinitePEPO (1-layer)" begin + Random.seed!(100) + pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) + pf = InfinitePartitionFunction(pepo) + env = CTMRGEnv(randn, ComplexF64, pf, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(pepo, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(pepo, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 + end +end From 0ce76376fe55c3cf83ea6c2ec92fbdecd67b1aca Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 11 Oct 2025 10:49:52 +0800 Subject: [PATCH 02/13] Add mixed state (iPEPO) correlator --- .../contractions/correlator/pepo_1layer.jl | 36 +++++++++++++------ .../contractions/correlator/peps.jl | 4 ++- test/utility/correlator.jl | 28 +++++++++------ 3 files changed, 47 insertions(+), 21 deletions(-) diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index 8cdff8519..adf803507 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -1,3 +1,17 @@ +function _transfer_left( + vec::GenericMPSTensor{S, 2}, d::MPSKit.SingleTransferMatrix + ) where {S} + return @tensor y[-1 -2; -3] := vec[1 2; 4] * + d.above[4 5; -3] * d.middle[2 3; 5 -2] * conj(d.below[1 3; -1]) +end + +function _transfer_left( + vec::GenericMPSTensor{S, 3}, d::MPSKit.SingleTransferMatrix + ) where {S} + return @tensor y[-1 dstring -2; -3] := vec[1 dstring 2; 4] * + d.above[4 5; -3] * d.middle[2 3; 5 -2] * conj(d.below[1 3; -1]) +end + function correlator_horizontal( ρ::InfinitePEPO, operator, i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, @@ -29,9 +43,8 @@ function correlator_horizontal( @tensor Amid[w s; n e] := t[d d; n e s w] Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) - twistdual!(T.below, 2:numout(T.below)) - Vn = Vn * T - Vo = Vo * T + Vn = _transfer_left(Vn, T) + Vo = _transfer_left(Vo, T) i += CartesianIndex(0, 1) end # compute overlap with operator @@ -42,8 +55,10 @@ function correlator_horizontal( @tensor Amid[w s; n e] := t[d d; n e s w] Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) - twistdual!(T.below, 2:numout(T.below)) - Vn = Vn * T + Vn = _transfer_left(Vn, T) + if k < length(js) + Vo = _transfer_left(Vo, T) + end i += CartesianIndex(0, 1) # compute overlap without operator denominator = end_correlator_denominator(j, Vn, env) @@ -74,8 +89,8 @@ function start_correlator( @autoopt @tensor Vo[χSE dstring De; χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * - removeunit(O, 1)[d1; d2 dstring] * - t[d2 d1; Dn De Ds Dw] * E_north[χN Dn; χNE] + removeunit(O, 1)[d2; d1 dstring] * + t[d1 d2; Dn De Ds Dw] * E_north[χN Dn; χNE] return Vn, Vo end @@ -96,7 +111,7 @@ function end_correlator_numerator( return @autoopt @tensor V[χSW dstring DW; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * - t[dt db; DN DE DS DW] * removeunit(O, 4)[dstring db; dt] + t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1] end function end_correlator_denominator( @@ -116,8 +131,9 @@ function correlator_vertical( env::CTMRGEnv, ) rotated_ρ = rotl90(ρ) - rotated_i = siterotl90(i, size(bra)) - rotated_js = map(j -> siterotl90(j, size(bra)), js) + unitcell = size(ρ)[1:2] + rotated_i = siterotl90(i, unitcell) + rotated_js = map(j -> siterotl90(j, unitcell), js) return correlator_horizontal( rotated_ρ, operator, rotated_i, rotated_js, rotl90(env) ) diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index f4c7bc93a..3122f45ac 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -52,7 +52,9 @@ function correlator_horizontal( 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 + if k < length(js) + Vo = Vo * T + end twistdual!(T.below, 2:numout(T.below)) Vn = Vn * T i += CartesianIndex(0, 1) diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl index 7658901b4..d45a0a0de 100644 --- a/test/utility/correlator.jl +++ b/test/utility/correlator.jl @@ -3,23 +3,29 @@ using Random using TensorKit using PEPSKit -Nr, Nc = 2, 2 -Vphy = Vect[FermionParity](0 => 1, 1 => 1) -Pspaces = fill(Vphy, (Nr, Nc)) -op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) +const syms = (Z2Irrep, FermionParity) -V = Vect[FermionParity](0 => 1, 1 => 2) -Venv = Vect[FermionParity](0 => 2, 1 => 2) -Nspaces = [V' V; V V'] -Espaces = [V V'; V' V] +function get_spaces(sym::Type{<:Sector}) + @assert sym in syms + Nr, Nc = 2, 2 + Vphy = Vect[sym](0 => 1, 1 => 1) + Pspaces = fill(Vphy, (Nr, Nc)) + V = Vect[sym](0 => 1, 1 => 2) + Venv = Vect[sym](0 => 2, 1 => 2) + Nspaces = [V' V; V V'] + Espaces = [V V'; V' V] + return Vphy, Venv, Pspaces, Nspaces, Espaces +end site0 = CartesianIndex(1, 1) maxsep = 6 site1xs = collect(site0 + CartesianIndex(0, i) for i in 2:2:maxsep) site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep) -@testset "Correlator in InfinitePEPS" begin +@testset "Correlator in InfinitePEPS ($(sym))" for sym in syms Random.seed!(100) + Vphy, Venv, Pspaces, Nspaces, Espaces = get_spaces(sym) + op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) env = CTMRGEnv(randn, ComplexF64, peps, Venv) for site1s in (site1xs, site1ys) @@ -34,8 +40,10 @@ site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep) end end -@testset "Correlator in InfinitePEPO (1-layer)" begin +@testset "Correlator in 1-layer InfinitePEPO ($(sym))" for sym in syms Random.seed!(100) + Vphy, Venv, Pspaces, Nspaces, Espaces = get_spaces(sym) + op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) pf = InfinitePartitionFunction(pepo) env = CTMRGEnv(randn, ComplexF64, pf, Venv) From 783f6a096bb14793a74f0c433288927736d774f4 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 11 Oct 2025 11:05:52 +0800 Subject: [PATCH 03/13] Fix fermionic twists --- src/algorithms/contractions/correlator/pepo_1layer.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index adf803507..ae7313b3c 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -39,8 +39,7 @@ function correlator_horizontal( # transfer until left of site j while j > i Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - t = ρ[mod1(i[1], end), mod1(i[2], end)] - @tensor Amid[w s; n e] := t[d d; n e s w] + Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) Vn = _transfer_left(Vn, T) @@ -51,8 +50,7 @@ function correlator_horizontal( numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) # transfer right of site j Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - t = ρ[mod1(i[1], end), mod1(i[2], end)] - @tensor Amid[w s; n e] := t[d d; n e s w] + Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) Vn = _transfer_left(Vn, T) @@ -79,7 +77,7 @@ function start_correlator( 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)] - t = ρ[mod1(r, end), mod1(c, end)] + t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further @autoopt @tensor Vn[χSE De; χNE] := @@ -107,7 +105,7 @@ function end_correlator_numerator( 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)] - t = ρ[mod1(r, end), mod1(c, end)] + t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) return @autoopt @tensor V[χSW dstring DW; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * From 2429c38bb56e2f1a8b92b0d0c23b40525cad8a35 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Oct 2025 15:34:21 +0800 Subject: [PATCH 04/13] Separate correlator contractions and logic --- src/PEPSKit.jl | 1 + .../contractions/correlator/pepo_1layer.jl | 23 ------- .../contractions/correlator/peps.jl | 34 ---------- src/algorithms/correlators.jl | 62 +++++++++++++++++++ 4 files changed, 63 insertions(+), 57 deletions(-) create mode 100644 src/algorithms/correlators.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index b279315d6..458bca995 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -76,6 +76,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") diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index ae7313b3c..0fb9e2737 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -136,26 +136,3 @@ function correlator_vertical( rotated_ρ, operator, rotated_i, rotated_js, rotl90(env) ) end - -function MPSKit.correlator( - ρ::InfinitePEPO, O, - i::CartesianIndex{2}, js::CoordCollection{2}, - env::CTMRGEnv, - ) - js = vec(js) # map CartesianIndices to Vector instead of Matrix - if all(==(i[1]) ∘ first ∘ Tuple, js) - return correlator_horizontal(ρ, O, i, js, env) - elseif all(==(i[2]) ∘ last ∘ Tuple, js) - return correlator_vertical(ρ, O, i, js, env) - else - error("Only horizontal or vertical correlators are implemented") - end -end - -function MPSKit.correlator( - ρ::InfinitePEPO, O, - i::CartesianIndex{2}, j::CartesianIndex{2}, - env::CTMRGEnv, - ) - return only(correlator(ρ, O, i, j:j, env)) -end diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index 3122f45ac..d27710c52 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -166,37 +166,3 @@ function correlator_vertical( rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env) ) end - -const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianIndices{N}} - -function MPSKit.correlator( - bra::InfinitePEPS, - O, - i::CartesianIndex{2}, js::CoordCollection{2}, - ket::InfinitePEPS, - env::CTMRGEnv, - ) - 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("Only horizontal or vertical correlators are implemented") - end -end - -function MPSKit.correlator( - bra::InfinitePEPS, - O, - i::CartesianIndex{2}, j::CartesianIndex{2}, - ket::InfinitePEPS, - env::CTMRGEnv, - ) - return only(correlator(bra, O, i, j:j, ket, env)) -end - -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/src/algorithms/correlators.jl b/src/algorithms/correlators.jl new file mode 100644 index 000000000..2f7322ed6 --- /dev/null +++ b/src/algorithms/correlators.jl @@ -0,0 +1,62 @@ +const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianIndices{N}} + +# Correlators in InfinitePEPS + +function MPSKit.correlator( + bra::InfinitePEPS, + O, + i::CartesianIndex{2}, js::CoordCollection{2}, + ket::InfinitePEPS, + env::CTMRGEnv, + ) + 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("Only horizontal or vertical correlators are implemented") + end +end + +function MPSKit.correlator( + bra::InfinitePEPS, + O, + i::CartesianIndex{2}, j::CartesianIndex{2}, + ket::InfinitePEPS, + env::CTMRGEnv, + ) + return only(correlator(bra, O, i, j:j, ket, env)) +end + +function MPSKit.correlator(state::InfinitePEPS, O, i::CartesianIndex{2}, j, env::CTMRGEnv) + return MPSKit.correlator(state, O, i, j, state, env) +end + +# Correlators in InfinitePEPO (tr(ρO)) + +function MPSKit.correlator( + ρ::InfinitePEPO, O, + i::CartesianIndex{2}, js::CoordCollection{2}, + env::CTMRGEnv, + ) + js = vec(js) # map CartesianIndices to Vector instead of Matrix + if all(==(i[1]) ∘ first ∘ Tuple, js) + return correlator_horizontal(ρ, O, i, js, env) + elseif all(==(i[2]) ∘ last ∘ Tuple, js) + return correlator_vertical(ρ, O, i, js, env) + else + error("Only horizontal or vertical correlators are implemented") + end +end + +function MPSKit.correlator( + ρ::InfinitePEPO, O, + i::CartesianIndex{2}, j::CartesianIndex{2}, + env::CTMRGEnv, + ) + return only(correlator(ρ, O, i, j:j, env)) +end + +# TODO: Correlators in InfinitePEPO (⟨ρ|O|ρ⟩) From 6a19152b198428db0666d7212582829b9fb2d8fb Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Oct 2025 16:42:45 +0800 Subject: [PATCH 05/13] [WIP] Use `V * T` syntax [skip ci] --- src/PEPSKit.jl | 2 +- .../contractions/correlator/pepo_1layer.jl | 37 ++++++++++--------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 458bca995..c9f3d93f8 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -14,7 +14,7 @@ using ChainRulesCore, Zygote using LoggingExtras using MPSKit -using MPSKit: MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix +using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace import MPSKit: infinite_temperature_density_matrix diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index 0fb9e2737..b89776869 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -1,15 +1,16 @@ -function _transfer_left( - vec::GenericMPSTensor{S, 2}, d::MPSKit.SingleTransferMatrix - ) where {S} - return @tensor y[-1 -2; -3] := vec[1 2; 4] * - d.above[4 5; -3] * d.middle[2 3; 5 -2] * conj(d.below[1 3; -1]) +function MPSKit.transfer_left( + vec::AbstractTensorMap{T, S, 1, 2}, + O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} + ) where {T, S} + return @tensor y[-1; -2 -3] := vec[1 2; 4] * + A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end -function _transfer_left( - vec::GenericMPSTensor{S, 3}, d::MPSKit.SingleTransferMatrix +function MPSKit.transfer_left( + v::GenericMPSTensor{S, 3}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} ) where {S} - return @tensor y[-1 dstring -2; -3] := vec[1 dstring 2; 4] * - d.above[4 5; -3] * d.middle[2 3; 5 -2] * conj(d.below[1 3; -1]) + return @tensor t[d_string -1 -2; -3] := v[d_string 1 2; 4] * + A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end function correlator_horizontal( @@ -42,8 +43,8 @@ function correlator_horizontal( Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) - Vn = _transfer_left(Vn, T) - Vo = _transfer_left(Vo, T) + Vo = Vo * T + Vn = Vn * T i += CartesianIndex(0, 1) end # compute overlap with operator @@ -53,10 +54,10 @@ function correlator_horizontal( Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] T = TransferMatrix(Atop, Amid, _dag(Abot)) - Vn = _transfer_left(Vn, T) if k < length(js) - Vo = _transfer_left(Vo, T) + Vo = Vo * T end + Vn = Vn * T i += CartesianIndex(0, 1) # compute overlap without operator denominator = end_correlator_denominator(j, Vn, env) @@ -80,11 +81,11 @@ function start_correlator( t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further - @autoopt @tensor Vn[χSE De; χNE] := + @autoopt @tensor Vn[χSE; De χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * t[d d; Dn De Ds Dw] * E_north[χN Dn; χNE] - @autoopt @tensor Vo[χSE dstring De; χNE] := + @autoopt @tensor Vo[dstring χSE De; χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * removeunit(O, 1)[d2; d1 dstring] * @@ -106,20 +107,20 @@ function end_correlator_numerator( C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) - return @autoopt @tensor V[χSW dstring DW; χNW] * + return @autoopt @tensor V[dstring χSW DW; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1] end function end_correlator_denominator( - j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 2, 1}, env::CTMRGEnv + j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 1, 2}, env::CTMRGEnv ) where {T, S} 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 DE; χN] * C_northeast[χN; χNE] * + return @autoopt @tensor V[χS; DE χN] * C_northeast[χN; χNE] * E_east[χNE DE; χSE] * C_southeast[χSE; χS] end From f0cf27fa45fd2efe91c1cc246ca894d4006908df Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Oct 2025 17:07:22 +0800 Subject: [PATCH 06/13] Move `correlator_horizontal/vertical` [skip ci] --- .../contractions/correlator/pepo_1layer.jl | 71 +-------- .../contractions/correlator/peps.jl | 88 ----------- src/algorithms/correlators.jl | 144 ++++++++++++++++++ 3 files changed, 146 insertions(+), 157 deletions(-) diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index b89776869..e1db97ab1 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -1,8 +1,8 @@ function MPSKit.transfer_left( - vec::AbstractTensorMap{T, S, 1, 2}, + vec::AbstractTensorMap{T, S, 1, 2}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} ) where {T, S} - return @tensor y[-1; -2 -3] := vec[1 2; 4] * + return @tensor y[-1; -2 -3] := vec[1; 2 4] * A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end @@ -13,59 +13,6 @@ function MPSKit.transfer_left( A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end -function correlator_horizontal( - ρ::InfinitePEPO, operator, - i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, - env::CTMRGEnv - ) - (size(ρ, 3) == 1) || - throw(ArgumentError("The input PEPO ρ must have only one layer.")) - 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(ρ), scalartype(env), scalartype.(O)... - ), - ) - # left start for operator and norm contractions - Vn, Vo = start_correlator(i, ρ, O[1], 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)] - Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] - T = TransferMatrix(Atop, Amid, _dag(Abot)) - Vo = Vo * T - Vn = Vn * T - i += CartesianIndex(0, 1) - end - # compute overlap with operator - numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) - # transfer right of site j - Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] - T = TransferMatrix(Atop, Amid, _dag(Abot)) - if k < length(js) - Vo = Vo * T - end - Vn = Vn * T - i += CartesianIndex(0, 1) - # compute overlap without operator - denominator = end_correlator_denominator(j, Vn, env) - G[k] = numerator / denominator - end - return G -end - function start_correlator( i::CartesianIndex{2}, ρ::InfinitePEPO, O::MPOTensor, env::CTMRGEnv @@ -123,17 +70,3 @@ function end_correlator_denominator( return @autoopt @tensor V[χS; DE χN] * C_northeast[χN; χNE] * E_east[χNE DE; χSE] * C_southeast[χSE; χS] end - -function correlator_vertical( - ρ::InfinitePEPO, operator, - i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, - env::CTMRGEnv, - ) - rotated_ρ = rotl90(ρ) - unitcell = size(ρ)[1:2] - rotated_i = siterotl90(i, unitcell) - rotated_js = map(j -> siterotl90(j, unitcell), js) - return correlator_horizontal( - rotated_ρ, operator, rotated_i, rotated_js, rotl90(env) - ) -end diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index d27710c52..697ad2a72 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -1,73 +1,3 @@ -function correlator_horizontal( - bra::InfinitePEPS, - operator, - i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, - ket::InfinitePEPS, - env::CTMRGEnv, - ) - 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 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)], - ) - T = TransferMatrix(Atop, sandwich, _dag(Abot)) - Vo = Vo * T - twistdual!(T.below, 2:numout(T.below)) - 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)) - if k < length(js) - Vo = Vo * T - end - twistdual!(T.below, 2:numout(T.below)) - Vn = Vn * T - i += CartesianIndex(0, 1) - - # compute overlap without operator - denominator = end_correlator_denominator(j, Vn, env) - - G[k] = numerator / denominator - end - - return G -end - function start_correlator( i::CartesianIndex{2}, below::InfinitePEPS, @@ -148,21 +78,3 @@ function end_correlator_denominator( E_east[χNE DEt DEb; χSE] * C_southeast[χSE; χS] end - -function correlator_vertical( - bra::InfinitePEPS, - operator, - i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, - ket::InfinitePEPS, - env::CTMRGEnv, - ) - rotated_bra = rotl90(bra) - rotated_ket = bra === ket ? rotated_bra : rotl90(ket) - - rotated_i = siterotl90(i, size(bra)) - rotated_js = map(j -> siterotl90(j, size(bra)), js) - - return correlator_horizontal( - rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env) - ) -end diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 2f7322ed6..951fdecce 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -2,6 +2,83 @@ const CoordCollection{N} = Union{AbstractVector{CartesianIndex{N}}, CartesianInd # Correlators in InfinitePEPS +function correlator_horizontal( + bra::InfinitePEPS, + operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + ket::InfinitePEPS, + env::CTMRGEnv, + ) + 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 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)], + ) + T = TransferMatrix(Atop, sandwich, _dag(Abot)) + Vo = Vo * T + twistdual!(T.below, 2:numout(T.below)) + 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)) + if k < length(js) + Vo = Vo * T + end + twistdual!(T.below, 2:numout(T.below)) + Vn = Vn * T + i += CartesianIndex(0, 1) + # compute overlap without operator + denominator = end_correlator_denominator(j, Vn, env) + G[k] = numerator / denominator + end + return G +end + +function correlator_vertical( + bra::InfinitePEPS, + operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + ket::InfinitePEPS, + env::CTMRGEnv, + ) + rotated_bra = rotl90(bra) + rotated_ket = bra === ket ? rotated_bra : rotl90(ket) + rotated_i = siterotl90(i, size(bra)) + rotated_js = map(j -> siterotl90(j, size(bra)), js) + return correlator_horizontal( + rotated_bra, operator, rotated_i, rotated_js, rotated_ket, rotl90(env) + ) +end + function MPSKit.correlator( bra::InfinitePEPS, O, @@ -36,6 +113,73 @@ end # Correlators in InfinitePEPO (tr(ρO)) +function correlator_horizontal( + ρ::InfinitePEPO, operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + env::CTMRGEnv + ) + (size(ρ, 3) == 1) || + throw(ArgumentError("The input PEPO ρ must have only one layer.")) + 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(ρ), scalartype(env), scalartype.(O)... + ), + ) + # left start for operator and norm contractions + Vn, Vo = start_correlator(i, ρ, O[1], 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)] + Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = TransferMatrix(Atop, Amid, _dag(Abot)) + Vo = Vo * T + Vn = Vn * T + i += CartesianIndex(0, 1) + end + # compute overlap with operator + numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) + # transfer right of site j + Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = TransferMatrix(Atop, Amid, _dag(Abot)) + if k < length(js) + Vo = Vo * T + end + Vn = Vn * T + i += CartesianIndex(0, 1) + # compute overlap without operator + denominator = end_correlator_denominator(j, Vn, env) + G[k] = numerator / denominator + end + return G +end + +function correlator_vertical( + ρ::InfinitePEPO, operator, + i::CartesianIndex{2}, js::AbstractVector{CartesianIndex{2}}, + env::CTMRGEnv, + ) + rotated_ρ = rotl90(ρ) + unitcell = size(ρ)[1:2] + rotated_i = siterotl90(i, unitcell) + rotated_js = map(j -> siterotl90(j, unitcell), js) + return correlator_horizontal( + rotated_ρ, operator, rotated_i, rotated_js, rotl90(env) + ) +end + function MPSKit.correlator( ρ::InfinitePEPO, O, i::CartesianIndex{2}, js::CoordCollection{2}, From 1acded8c5ad0bbeb5e546a94272e4fceb58982c3 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 24 Oct 2025 18:56:50 +0800 Subject: [PATCH 07/13] Reorganize `transfer_left`s used by correlators --- .../contractions/correlator/pepo_1layer.jl | 26 +++++++++---------- .../contractions/correlator/peps.jl | 13 ++++++++++ .../contractions/vumps_contractions.jl | 13 ---------- 3 files changed, 25 insertions(+), 27 deletions(-) diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index e1db97ab1..bedcbd634 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -1,15 +1,14 @@ function MPSKit.transfer_left( - vec::AbstractTensorMap{T, S, 1, 2}, - O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} + vec::AbstractTensor{T, S, 3}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} ) where {T, S} - return @tensor y[-1; -2 -3] := vec[1; 2 4] * + return @tensor y[-1 -2 -3] := vec[1 2 4] * A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end function MPSKit.transfer_left( - v::GenericMPSTensor{S, 3}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} - ) where {S} - return @tensor t[d_string -1 -2; -3] := v[d_string 1 2; 4] * + v::AbstractTensor{T, S, 4}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} + ) where {T, S} + return @tensor t[d_string -1 -2 -3] := v[d_string 1 2 4] * A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) end @@ -28,11 +27,11 @@ function start_correlator( t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further - @autoopt @tensor Vn[χSE; De χNE] := + @autoopt @tensor Vn[χSE De χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * t[d d; Dn De Ds Dw] * E_north[χN Dn; χNE] - @autoopt @tensor Vo[dstring χSE De; χNE] := + @autoopt @tensor Vo[dstring χSE De χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * removeunit(O, 1)[d2; d1 dstring] * @@ -41,9 +40,8 @@ function start_correlator( end function end_correlator_numerator( - j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 3, 1}, - ρ::InfinitePEPO, - O::MPOTensor, env::CTMRGEnv + j::CartesianIndex{2}, V::AbstractTensor{T, S, 4}, + ρ::InfinitePEPO, O::MPOTensor, env::CTMRGEnv ) where {T, S} (size(ρ, 3) == 1) || throw(ArgumentError("The input PEPO ρ must have only one layer.")) @@ -54,19 +52,19 @@ function end_correlator_numerator( C_northeast = env.corners[NORTHEAST, _prev(r, end), _next(c, end)] C_southeast = env.corners[SOUTHEAST, _next(r, end), _next(c, end)] t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) - return @autoopt @tensor V[dstring χSW DW; χNW] * + return @autoopt @tensor V[dstring χSW DW χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1] end function end_correlator_denominator( - j::CartesianIndex{2}, V::AbstractTensorMap{T, S, 1, 2}, env::CTMRGEnv + j::CartesianIndex{2}, V::AbstractTensor{T, S, 3}, env::CTMRGEnv ) where {T, S} 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; DE χN] * C_northeast[χN; χNE] * + return @autoopt @tensor V[χS DE χN] * C_northeast[χN; χNE] * E_east[χNE DE; χSE] * C_southeast[χSE; χS] end diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index 697ad2a72..26102934d 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -1,3 +1,16 @@ +# transfer with excited GL +function MPSKit.transfer_left( + GL::GenericMPSTensor{S, 4}, O::PEPSSandwich, + A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, + ) where {S} + return @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 start_correlator( i::CartesianIndex{2}, below::InfinitePEPS, diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 97f231905..91d4bb06e 100644 --- a/src/algorithms/contractions/vumps_contractions.jl +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -32,19 +32,6 @@ 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} - return @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, A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, From d773a22867826ed70ac4d6f46a9d787fb0f9b6e3 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 29 Oct 2025 13:26:43 +0100 Subject: [PATCH 08/13] Attempt at edge transfer matrices --- src/PEPSKit.jl | 4 +- .../contractions/correlator/pepo_1layer.jl | 32 +--- .../contractions/correlator/peps.jl | 13 -- src/algorithms/contractions/transfer.jl | 178 ++++++++++++++++++ .../contractions/vumps_contractions.jl | 12 +- src/algorithms/correlators.jl | 30 ++- src/algorithms/ctmrg/gaugefix.jl | 4 +- src/algorithms/toolbox.jl | 48 +++-- src/algorithms/transfermatrix.jl | 38 ++++ 9 files changed, 285 insertions(+), 74 deletions(-) create mode 100644 src/algorithms/contractions/transfer.jl create mode 100644 src/algorithms/transfermatrix.jl diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index c9f3d93f8..989ef7bda 100644 --- a/src/PEPSKit.jl +++ b/src/PEPSKit.jl @@ -14,7 +14,7 @@ using ChainRulesCore, Zygote using LoggingExtras using MPSKit -using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, TransferMatrix +using MPSKit: MPSTensor, MPOTensor, GenericMPSTensor, MPSBondTensor, ProductTransferMatrix import MPSKit: tensorexpr, leading_boundary, loginit!, logiter!, logfinish!, logcancel!, physicalspace import MPSKit: infinite_temperature_density_matrix @@ -51,6 +51,7 @@ include("environments/vumps_environments.jl") include("environments/suweight.jl") include("algorithms/contractions/ctmrg_contractions.jl") +include("algorithms/contractions/transfer.jl") include("algorithms/contractions/localoperator.jl") include("algorithms/contractions/vumps_contractions.jl") include("algorithms/contractions/bondenv/benv_tools.jl") @@ -75,6 +76,7 @@ include("algorithms/time_evolution/evoltools.jl") include("algorithms/time_evolution/simpleupdate.jl") include("algorithms/time_evolution/simpleupdate3site.jl") +include("algorithms/transfermatrix.jl") include("algorithms/toolbox.jl") include("algorithms/correlators.jl") diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index bedcbd634..80b53a226 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -1,20 +1,6 @@ -function MPSKit.transfer_left( - vec::AbstractTensor{T, S, 3}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} - ) where {T, S} - return @tensor y[-1 -2 -3] := vec[1 2 4] * - A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) -end - -function MPSKit.transfer_left( - v::AbstractTensor{T, S, 4}, O::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S} - ) where {T, S} - return @tensor t[d_string -1 -2 -3] := v[d_string 1 2 4] * - A[4 5; -3] * O[2 3; 5 -2] * conj(Ab[1 3; -1]) -end - function start_correlator( i::CartesianIndex{2}, ρ::InfinitePEPO, - O::MPOTensor, env::CTMRGEnv + O::PFTensor, env::CTMRGEnv ) (size(ρ, 3) == 1) || throw(ArgumentError("The input PEPO ρ must have only one layer.")) @@ -27,11 +13,11 @@ function start_correlator( t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: part of these contractions is duplicated between the two output tensors, # so could be optimized further - @autoopt @tensor Vn[χSE De χNE] := + @autoopt @tensor Vn[χSE De; χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * t[d d; Dn De Ds Dw] * E_north[χN Dn; χNE] - @autoopt @tensor Vo[dstring χSE De χNE] := + @autoopt @tensor Vo[χSE De dstring; χNE] := E_south[χSE Ds; χSW2] * C_southwest[χSW2; χSW] * E_west[χSW Dw; χNW] * C_northwest[χNW; χN] * removeunit(O, 1)[d2; d1 dstring] * @@ -40,8 +26,8 @@ function start_correlator( end function end_correlator_numerator( - j::CartesianIndex{2}, V::AbstractTensor{T, S, 4}, - ρ::InfinitePEPO, O::MPOTensor, env::CTMRGEnv + j::CartesianIndex{2}, V::CTMRGEdgeTensor{T, S, 3}, + ρ::InfinitePEPO, O::PFTensor, env::CTMRGEnv ) where {T, S} (size(ρ, 3) == 1) || throw(ArgumentError("The input PEPO ρ must have only one layer.")) @@ -51,20 +37,20 @@ function end_correlator_numerator( 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)] - t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) - return @autoopt @tensor V[dstring χSW DW χNW] * + t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: is this still needed? + return @autoopt @tensor V[χSW DW dstring; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * t[d1 d2; DN DE DS DW] * removeunit(O, 4)[dstring d2; d1] end function end_correlator_denominator( - j::CartesianIndex{2}, V::AbstractTensor{T, S, 3}, env::CTMRGEnv + j::CartesianIndex{2}, V::CTMRGEdgeTensor{T, S, 2}, env::CTMRGEnv ) where {T, S} 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 DE χN] * C_northeast[χN; χNE] * + return @autoopt @tensor V[χS DE; χN] * C_northeast[χN; χNE] * E_east[χNE DE; χSE] * C_southeast[χSE; χS] end diff --git a/src/algorithms/contractions/correlator/peps.jl b/src/algorithms/contractions/correlator/peps.jl index 26102934d..697ad2a72 100644 --- a/src/algorithms/contractions/correlator/peps.jl +++ b/src/algorithms/contractions/correlator/peps.jl @@ -1,16 +1,3 @@ -# transfer with excited GL -function MPSKit.transfer_left( - GL::GenericMPSTensor{S, 4}, O::PEPSSandwich, - A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, - ) where {S} - return @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 start_correlator( i::CartesianIndex{2}, below::InfinitePEPS, diff --git a/src/algorithms/contractions/transfer.jl b/src/algorithms/contractions/transfer.jl new file mode 100644 index 000000000..3743396ed --- /dev/null +++ b/src/algorithms/contractions/transfer.jl @@ -0,0 +1,178 @@ +# +# Transfer function for (CTMRG) edges +# + +edge_transfer_left(v, ::Nothing, A, B) = edge_transfer_left(v, A, B) +edge_transfer_right(v, ::Nothing, A, B) = edge_transfer_right(v, A, B) + +""" + edge_transfer_left(v, Et, Eb) + +Apply an edge transfer matrix to the left. + +``` + ┌─Et─ +-v │ + └─qƎ─ +``` +""" +@generated function edge_transfer_left( + v::AbstractTensorMap{<:Any, S, 1, N₁}, + Etop::CTMRGEdgeTensor{<:Any, S, N₂}, + Ebot::CTMRGEdgeTensor{<:Any, S, N₂} + ) where {S, N₁, N₂} + t_out = tensorexpr(:v, -1, -(2:(N₁ + 1))) + t_top = tensorexpr(:Etop, 2:(N₂ + 1), -(N₁ + 1)) + t_bot = tensorexpr(:Ebot, (-1, (3:(N₂ + 1))...), 1) + t_in = tensorexpr(:v, 1, (-(2:N₁)..., 2)) + return macroexpand( + @__MODULE__, :(return @tensor $t_out := $t_in * $t_top * $t_bot) + ) +end + + +""" + edge_transfer_right(v, Et, Eb) + +Apply an edge transfer matrix to the right. + +``` +─Et─┐ + │ v- +─qƎ─┘ +``` +""" +@generated function edge_transfer_right( + v::AbstractTensorMap{<:Any, S, 1, N₁}, + Etop::CTMRGEdgeTensor{<:Any, S, N₂}, + Ebot::CTMRGEdgeTensor{<:Any, S, N₂} + ) where {S, N₁, N₂} + t_out = tensorexpr(:v, -1, -(2:(N₁ + 1))) + t_top = tensorexpr(:Etop, (-1, (3:(N₂ + 1))...), 1) + t_bot = tensorexpr(:Ebot, (2, (3:(N₂ + 1))...), -(N₁ + 1)) + t_in = tensorexpr(:v, 1, (-(2:N₁)..., 2)) + return macroexpand( + @__MODULE__, :(return @tensor $t_out := $t_top * $t_bot * $t_in) + ) +end + +""" + edge_transfer_left(v, O, Et, Eb) + +Apply an edge transfer matrix to the left. + +``` + ┌──Et─ + │ │ + v──O── + │ │ + └──qƎ─ +``` +""" +function edge_transfer_left( + v::CTMRGEdgeTensor{<:Any, S, 3}, + O::PEPSSandwich, + Etop::CTMRGEdgeTensor{<:Any, S, 3}, + Ebot::CTMRGEdgeTensor{<:Any, S, 3}, + ) where {S} + @autoopt @tensor v´[χ_SE D_E_above D_E_below; χ_NE] := + v[χ_SW D_W_above D_W_below; χ_NW] * + Etop[χ_NW D_N_above D_N_below; χ_NE] * + Ebot[χ_SE D_S_above D_S_below; χ_SW] * + 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]) + + return v´ +end +function edge_transfer_left( + v::CTMRGEdgeTensor{<:Any, S, 2}, + O::PFTensor, + Etop::CTMRGEdgeTensor{<:Any, S, 2}, + Ebot::CTMRGEdgeTensor{<:Any, S, 2}, + ) where {S} + @autoopt @tensor v´[χ_SE D_E; χ_NE] := + v[χ_SW D_W; χ_NW] * + Etop[χ_NW D_N; χ_NE] * + Ebot[χ_SE D_S; χ_SW] * + O[D_W D_S; D_N D_E] + + return v´ +end + +""" + transfer_right(v, Et, Eb) + +Apply an edge transfer matrix to the right. + +``` +──Et─┐ + │ │ +──O──v + │ │ +──qƎ─┘ +``` +""" +function edge_transfer_right( + v::CTMRGEdgeTensor{<:Any, S, 3}, + O::PEPSSandwich, + Etop::CTMRGEdgeTensor{<:Any, S, 3}, + Ebot::CTMRGEdgeTensor{<:Any, S, 3}, + ) where {S} + @autoopt @tensor v′[χ_NW D_W_above D_W_below; χ_SW] := + v[χ_NE D_E_above D_E_below; χ_SE] * + Etop[χ_NW D_N_above D_N_below; χ_NE] * + Ebot[χ_SE D_S_below D_S_above; χ_SW] * + 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]) + + return v′ +end +function edge_transfer_right( + v::CTMRGEdgeTensor{<:Any, S, 2}, + O::PFTensor, + Etop::CTMRGEdgeTensor{<:Any, S, 2}, + Ebot::CTMRGEdgeTensor{<:Any, S, 2}, + ) where {S} + return @autoopt @tensor v′[χ_NW D_W; χ_SW] := + v[χ_NE D_E; χ_SE] * + Etop[χ_NW D_N; χ_NE] * + Ebot[χ_SE D_S; χ_SW] * + O[D_W D_S; D_N D_E] + + return v′ +end + +""" + edge_transfer_left(v, O, Et, Eb) + +Apply an edge transfer matrix to the left on an excited vector. + +``` + ┌──Et─ + │ │ +-v──O── + │ │ + └──qƎ─ +``` +""" +function edge_transfer_left( + v::CTMRGEdgeTensor{<:Any, S, 4}, O::PEPSSandwich, + Etop::CTMRGEdgeTensor{<:Any, S, 3}, Ebot::CTMRGEdgeTensor{<:Any, S, 3}, + ) where {S} + return @autoopt @tensor v′[χ_SE D_E_above d_string D_E_below; χ_NE] := + v[χ_SW D_W_above d_string D_W_below; χ_NW] * + Etop[χ_NW D_N_above D_N_below; χ_NE] * + Ebot[χ_SE D_S_above D_S_below; χ_SW] * + 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]) +end +function edge_transfer_left( + v::CTMRGEdgeTensor{<:Any, S, 3}, O::PFTensor, + Etop::CTMRGEdgeTensor{<:Any, S, 2}, Ebot::CTMRGEdgeTensor{<:Any, S, 2}, + ) where {S} + return @autoopt @tensor v′[χ_SE D_E d_string; χ_NE] := + v[χ_SW D_W d_string; χ_NW] * + Etop[χ_NW D_N; χ_NE] * + Ebot[χ_SE D_S; χ_SW] * + O[D_W D_S; D_N D_E] +end diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 91d4bb06e..58584f7e6 100644 --- a/src/algorithms/contractions/vumps_contractions.jl +++ b/src/algorithms/contractions/vumps_contractions.jl @@ -7,7 +7,7 @@ function MPSKit.transfer_left( A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N} Ā = twistdual(Ā, 2:N) - return _transfer_left(GL, O, A, Ā) + return mps_transfer_left(GL, O, A, Ā) end function MPSKit.transfer_right( @@ -15,12 +15,12 @@ function MPSKit.transfer_right( A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N} Ā = twistdual(Ā, 2:N) - return _transfer_right(GR, O, A, Ā) + return mps_transfer_right(GR, O, A, Ā) end ## PEPS -function _transfer_left( +function mps_transfer_left( GL::GenericMPSTensor{S, 3}, O::PEPSSandwich, A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, ) where {S} @@ -32,7 +32,7 @@ function _transfer_left( A[χ_NW D_N_above D_N_below; χ_NE] end -function _transfer_right( +function mps_transfer_right( GR::GenericMPSTensor{S, 3}, O::PEPSSandwich, A::GenericMPSTensor{S, 3}, Ā::GenericMPSTensor{S, 3}, ) where {S} @@ -46,7 +46,7 @@ end ## PEPO -@generated function _transfer_left( +@generated function mps_transfer_left( GL::GenericMPSTensor{S, N}, O::PEPOSandwich{H}, A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N, H} @@ -70,7 +70,7 @@ end return macroexpand(@__MODULE__, :(return @autoopt @tensor $GL´_e := $rhs)) end -@generated function _transfer_right( +@generated function mps_transfer_right( GR::GenericMPSTensor{S, N}, O::PEPOSandwich{H}, A::GenericMPSTensor{S, N}, Ā::GenericMPSTensor{S, N}, ) where {S, N, H} diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl index 951fdecce..674644aab 100644 --- a/src/algorithms/correlators.jl +++ b/src/algorithms/correlators.jl @@ -30,30 +30,28 @@ function correlator_horizontal( 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)] + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Ebot = 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)) + T = edge_transfermatrix(Etop, sandwich, Ebot) Vo = Vo * T - twistdual!(T.below, 2:numout(T.below)) 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)] + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Ebot = 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)) + T = edge_transfermatrix(Etop, sandwich, Ebot) if k < length(js) Vo = Vo * T end - twistdual!(T.below, 2:numout(T.below)) Vn = Vn * T i += CartesianIndex(0, 1) # compute overlap without operator @@ -139,10 +137,10 @@ function correlator_horizontal( 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)] - Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] - T = TransferMatrix(Atop, Amid, _dag(Abot)) + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = edge_transfermatrix(Etop, Omid, Ebot) Vo = Vo * T Vn = Vn * T i += CartesianIndex(0, 1) @@ -150,10 +148,10 @@ function correlator_horizontal( # compute overlap with operator numerator = end_correlator_numerator(j, Vo, ρ, O[2], env) # transfer right of site j - Atop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] - Amid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) - Abot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] - T = TransferMatrix(Atop, Amid, _dag(Abot)) + Etop = env.edges[NORTH, _prev(i[1], end), mod1(i[2], end)] + Omid = trace_physicalspaces(ρ[mod1(i[1], end), mod1(i[2], end)]) + Ebot = env.edges[SOUTH, _next(i[1], end), mod1(i[2], end)] + T = edge_transfermatrix(Etop, Omid, Ebot) if k < length(js) Vo = Vo * T end diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index c5aa9e504..8e03b6405 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -56,7 +56,7 @@ end # this is a bit of a hack to get the fixed point of the mixed transfer matrix # because MPSKit is not compatible with AD -@generated function _transfer_right( +@generated function mps_transfer_right( v::AbstractTensorMap{<:Any, S, 1, N₁}, A::GenericMPSTensor{S, N₂}, Abar::GenericMPSTensor{S, N₂}, ) where {S, N₁, N₂} @@ -71,7 +71,7 @@ end function transfermatrix_fixedpoint(tops, bottoms, ρinit) _, vecs, info = eigsolve(ρinit, 1, :LM, Arnoldi()) do ρ return foldr(zip(tops, bottoms); init = ρ) do (top, bottom), ρ - return _transfer_right(ρ, top, bottom) + return mps_transfer_right(ρ, top, bottom) end end info.converged > 0 || @warn "eigsolve did not converge" diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index cdd2adb6d..496494263 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -240,14 +240,36 @@ end end """ -Adjoint of an MPS tensor, but permutes the physical spaces back into the codomain. -Intuitively, this conjugates a tensor and then reinterprets its 'direction' as an MPS -tensor. + edge_transfer_spectrum(top::Vector{E}, bot::Vector{E}; tol=Defaults.tol, num_vals=20, + sector=one(sectortype(E))) where {E<:CTMRGEdgeTensor} + +Calculate the partial spectrum of the left edge transfer matrix corresponding to the given +`top` vector of edges and a `bot` vector of edge. The `sector` keyword argument can be used +to specify a non-trivial total charge for the transfer matrix eigenvectors. Specifically, an +auxiliary space `ℂ[typeof(sector)](sector => 1)'` will be added to the domain of each +eigenvector. The `tol` and `num_vals` keyword arguments are passed to `KrylovKit.eigolve`. """ -function _dag(A::MPSKit.GenericMPSTensor{S, N}) where {S, N} - return permute(A', ((1, (3:(N + 1))...), (2,))) +function edge_transfer_spectrum( + top::Vector{E}, bot::Vector{E}; tol = MPSKit.Defaults.tol, num_vals = 20, + sector = one(sectortype(E)) + ) where {E <: CTMRGEdgeTensor} + init = randn( + scalartype(E), + MPSKit._lastspace(first(bot))' ← ℂ[typeof(sector)](sector => 1)' ⊗ MPSKit._firstspace(first(top)), + ) + + transferspace = fuse(MPSKit._firstspace(first(top)) * MPSKit._lastspace(first(bot))) + num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values + eigenvals, eigenvecs, convhist = eigsolve( + flip(edge_transfermatrix(top, bot)), init, num_vals, :LM; tol = tol + ) + convhist.converged < num_vals && + @warn "correlation length failed to converge: normres = $(convhist.normres)" + + return eigenvals end + # TODO: decide on appropriate signature and returns for the more generic case """ correlation_length(state, env::CTMRGEnv; num_vals=2, kwargs...) @@ -271,15 +293,15 @@ function _correlation_length( # Horizontal λ_h = map(1:n_rows) do r - above = InfiniteMPS(env.edges[NORTH, r, :]) - below = InfiniteMPS(_dag.(env.edges[SOUTH, _next(r, n_rows), :])) - vals = MPSKit.transfer_spectrum(above; below, num_vals, sector, kwargs...) + top = env.edges[NORTH, r, :] + bot = env.edges[SOUTH, _next(r, n_rows), :] + vals = edge_transfer_spectrum(top, bot; num_vals, sector, kwargs...) # normalize using largest eigenvalue in trivial sector if isone(sector) N = first(vals) else - vals_triv = MPSKit.transfer_spectrum(above; below, num_vals = 1, kwargs...) + vals_triv = edge_transfer_spectrum(top, bot; num_vals = 1, kwargs...) N = first(vals_triv) end return vals ./ N # normalize largest eigenvalue @@ -287,15 +309,15 @@ function _correlation_length( # Vertical λ_v = map(1:n_cols) do c - above = InfiniteMPS(env.edges[EAST, :, c]) - below = InfiniteMPS(_dag.(env.edges[WEST, :, _next(c, n_cols)])) - vals = MPSKit.transfer_spectrum(above; below, num_vals, sector, kwargs...) + top = env.edges[EAST, :, c] + bot = env.edges[WEST, :, _next(c, n_cols)] + vals = edge_transfer_spectrum(top, bot; num_vals, sector, kwargs...) # normalize using largest eigenvalue in trivial sector if isone(sector) N = first(vals) else - vals_triv = MPSKit.transfer_spectrum(above; below, num_vals = 1, kwargs...) + vals_triv = edge_transfer_spectrum(top, bot; num_vals = 1, kwargs...) N = first(vals_triv) end return vals ./ N # normalize largest eigenvalue diff --git a/src/algorithms/transfermatrix.jl b/src/algorithms/transfermatrix.jl new file mode 100644 index 000000000..354d3fe73 --- /dev/null +++ b/src/algorithms/transfermatrix.jl @@ -0,0 +1,38 @@ +# +# Edge transfer matrices +# + +# single site transfer +struct EdgeTransferMatrix{A <: CTMRGEdgeTensor, B, C <: CTMRGEdgeTensor} <: MPSKit.AbstractTransferMatrix + top::A + mid::B + bot::C + isflipped::Bool +end + +Base.:*(tm1::T, tm2::T) where {T <: EdgeTransferMatrix} = ProductTransferMatrix([tm1, tm2]) + +# TODO: really not sure it TensorKit.flip is the suitable method for this... +function TensorKit.flip(tm::EdgeTransferMatrix) + return EdgeTransferMatrix(tm.top, tm.mid, tm.bot, !tm.isflipped) +end + +# action on a vector using * is dispatched in MPSKit to regular function application +function (d::EdgeTransferMatrix)(vec) + return if d.isflipped + edge_transfer_left(vec, d.mid, d.top, d.bot) + else + edge_transfer_right(vec, d.mid, d.top, d.bot) + end +end + +# constructors +edge_transfermatrix(a) = edge_transfermatrix(a, nothing, a) +edge_transfermatrix(a, b) = edge_transfermatrix(a, nothing, b) +function edge_transfermatrix(a::CTMRGEdgeTensor, b, c::CTMRGEdgeTensor, isflipped = false) + return EdgeTransferMatrix(a, b, c, isflipped) +end +function edge_transfermatrix(a::AbstractVector, b, c::AbstractVector, isflipped = false) + tot = ProductTransferMatrix(convert(Vector, edge_transfermatrix.(a, b, c))) + return isflipped ? flip(tot) : tot +end From a25b06fbc8fd6fcd1287dd6b139f02325effd839 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 30 Oct 2025 16:22:26 +0800 Subject: [PATCH 09/13] Test correlators with dual physical space --- test/utility/correlator.jl | 59 +++++++++++++++++++++----------------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl index d45a0a0de..c22d3e88a 100644 --- a/test/utility/correlator.jl +++ b/test/utility/correlator.jl @@ -9,12 +9,11 @@ function get_spaces(sym::Type{<:Sector}) @assert sym in syms Nr, Nc = 2, 2 Vphy = Vect[sym](0 => 1, 1 => 1) - Pspaces = fill(Vphy, (Nr, Nc)) V = Vect[sym](0 => 1, 1 => 2) Venv = Vect[sym](0 => 2, 1 => 2) Nspaces = [V' V; V V'] Espaces = [V V'; V' V] - return Vphy, Venv, Pspaces, Nspaces, Espaces + return Vphy, Venv, Nspaces, Espaces end site0 = CartesianIndex(1, 1) @@ -24,37 +23,43 @@ site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep) @testset "Correlator in InfinitePEPS ($(sym))" for sym in syms Random.seed!(100) - Vphy, Venv, Pspaces, Nspaces, Espaces = get_spaces(sym) - op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) - peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) - env = CTMRGEnv(randn, ComplexF64, peps, Venv) - for site1s in (site1xs, site1ys) - vals1 = correlator(peps, op, site0, site1s, env) - vals2 = map(site1s) do site1 - O = LocalOperator(Pspaces, (site0, site1) => op) - return expectation_value(peps, O, env) + Vphy, Venv, Nspaces, Espaces = get_spaces(sym) + for Vp in [Vphy', Vphy] + op = randn(ComplexF64, Vp ⊗ Vp → Vp ⊗ Vp) + Pspaces = fill(Vp, size(Nspaces)) + peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) + env = CTMRGEnv(randn, ComplexF64, peps, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(peps, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(peps, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 end - @info vals1 - @info vals2 - @test vals1 ≈ vals2 end end @testset "Correlator in 1-layer InfinitePEPO ($(sym))" for sym in syms Random.seed!(100) - Vphy, Venv, Pspaces, Nspaces, Espaces = get_spaces(sym) - op = randn(ComplexF64, Vphy ⊗ Vphy → Vphy ⊗ Vphy) - pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) - pf = InfinitePartitionFunction(pepo) - env = CTMRGEnv(randn, ComplexF64, pf, Venv) - for site1s in (site1xs, site1ys) - vals1 = correlator(pepo, op, site0, site1s, env) - vals2 = map(site1s) do site1 - O = LocalOperator(Pspaces, (site0, site1) => op) - return expectation_value(pepo, O, env) + Vphy, Venv, Nspaces, Espaces = get_spaces(sym) + for Vp in [Vphy', Vphy] + op = randn(ComplexF64, Vp ⊗ Vp → Vp ⊗ Vp) + Pspaces = fill(Vp, size(Nspaces)) + pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) + pf = InfinitePartitionFunction(pepo) + env = CTMRGEnv(randn, ComplexF64, pf, Venv) + for site1s in (site1xs, site1ys) + vals1 = correlator(pepo, op, site0, site1s, env) + vals2 = map(site1s) do site1 + O = LocalOperator(Pspaces, (site0, site1) => op) + return expectation_value(pepo, O, env) + end + @info vals1 + @info vals2 + @test vals1 ≈ vals2 end - @info vals1 - @info vals2 - @test vals1 ≈ vals2 end end From 5e0046cf1939d5b3b39f5eccc249015f48287833 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Tue, 4 Nov 2025 14:11:55 +0800 Subject: [PATCH 10/13] Remove dual physical space test for now --- test/utility/correlator.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/utility/correlator.jl b/test/utility/correlator.jl index c22d3e88a..766457c23 100644 --- a/test/utility/correlator.jl +++ b/test/utility/correlator.jl @@ -24,7 +24,8 @@ site1ys = collect(site0 + CartesianIndex(i, 0) for i in 2:2:maxsep) @testset "Correlator in InfinitePEPS ($(sym))" for sym in syms Random.seed!(100) Vphy, Venv, Nspaces, Espaces = get_spaces(sym) - for Vp in [Vphy', Vphy] + # TODO: test dual physical space + for Vp in [Vphy] op = randn(ComplexF64, Vp ⊗ Vp → Vp ⊗ Vp) Pspaces = fill(Vp, size(Nspaces)) peps = InfinitePEPS(randn, ComplexF64, Pspaces, Nspaces, Espaces) @@ -45,7 +46,8 @@ end @testset "Correlator in 1-layer InfinitePEPO ($(sym))" for sym in syms Random.seed!(100) Vphy, Venv, Nspaces, Espaces = get_spaces(sym) - for Vp in [Vphy', Vphy] + # TODO: test dual physical space + for Vp in [Vphy] op = randn(ComplexF64, Vp ⊗ Vp → Vp ⊗ Vp) Pspaces = fill(Vp, size(Nspaces)) pepo = InfinitePEPO(randn, ComplexF64, Pspaces, Nspaces, Espaces) From a4468e952d1be220c483abadf9be522328124489 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 4 Nov 2025 07:41:17 +0100 Subject: [PATCH 11/13] Don't unnecessarily depend on MPSKit.jl internals --- src/algorithms/ctmrg/gaugefix.jl | 2 +- src/algorithms/toolbox.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index ded2ba525..d8e06a482 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -38,7 +38,7 @@ function gauge_fix(envprev::CTMRGEnv{C, T}, envfinal::CTMRGEnv{C, T}) where {C, # Find right fixed points of mixed transfer matrices ρinit = randn( - scalartype(T), MPSKit._lastspace(Tsfinal[end])' ← MPSKit._lastspace(M[end])' + scalartype(T), space(Tsfinal[end], numind(Tsfinal[end]))' ← space(M[end], numind(M[end]))' ) ρprev = transfermatrix_fixedpoint(Tsprev, M, ρinit) ρfinal = transfermatrix_fixedpoint(Tsfinal, M, ρinit) diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 496494263..13d4652ae 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -255,10 +255,10 @@ function edge_transfer_spectrum( ) where {E <: CTMRGEdgeTensor} init = randn( scalartype(E), - MPSKit._lastspace(first(bot))' ← ℂ[typeof(sector)](sector => 1)' ⊗ MPSKit._firstspace(first(top)), + space(first(bot), numind(first(bot)))' ← ℂ[typeof(sector)](sector => 1)' ⊗ space(first(top), 1), ) - transferspace = fuse(MPSKit._firstspace(first(top)) * MPSKit._lastspace(first(bot))) + transferspace = fuse(space(first(top), 1) * space(first(bot), numind(first(bot)))') num_vals = min(dim(transferspace, sector), num_vals) # we can ask at most this many values eigenvals, eigenvecs, convhist = eigsolve( flip(edge_transfermatrix(top, bot)), init, num_vals, :LM; tol = tol From ef134904947f7b72a11ad5fb0b4b436d95b23d11 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Thu, 6 Nov 2025 22:59:48 +0800 Subject: [PATCH 12/13] Remove unneeded TODO --- src/algorithms/contractions/correlator/pepo_1layer.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/algorithms/contractions/correlator/pepo_1layer.jl b/src/algorithms/contractions/correlator/pepo_1layer.jl index 80b53a226..a8e7f0490 100644 --- a/src/algorithms/contractions/correlator/pepo_1layer.jl +++ b/src/algorithms/contractions/correlator/pepo_1layer.jl @@ -37,7 +37,7 @@ function end_correlator_numerator( 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)] - t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) # TODO: is this still needed? + t = twistdual(ρ[mod1(r, end), mod1(c, end)], 1:2) return @autoopt @tensor V[χSW DW dstring; χNW] * E_south[χSSE DS; χSW] * E_east[χNEE DE; χSEE] * E_north[χNW DN; χNNE] * C_northeast[χNNE; χNEE] * C_southeast[χSEE; χSSE] * From a94747f984339ae7bdaea4f34be0dbd0aa7708fb Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 6 Nov 2025 16:08:34 +0100 Subject: [PATCH 13/13] Add note --- src/algorithms/ctmrg/gaugefix.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/algorithms/ctmrg/gaugefix.jl b/src/algorithms/ctmrg/gaugefix.jl index d8e06a482..5cfa17106 100644 --- a/src/algorithms/ctmrg/gaugefix.jl +++ b/src/algorithms/ctmrg/gaugefix.jl @@ -56,6 +56,9 @@ end # this is a bit of a hack to get the fixed point of the mixed transfer matrix # because MPSKit is not compatible with AD +# NOTE: the action of the transfer operator here is NOT the same as that of +# MPSKit.transfer_right; to match the MPSKit implementation, it would require a twist on any +# index where spaces(Abar, 2:end) is dual (which we can just ignore as long as we use a random Abar) @generated function mps_transfer_right( v::AbstractTensorMap{<:Any, S, 1, N₁}, A::GenericMPSTensor{S, N₂}, Abar::GenericMPSTensor{S, N₂},