diff --git a/src/PEPSKit.jl b/src/PEPSKit.jl index 02253a719..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,6 +69,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") @@ -80,6 +85,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 export leading_boundary export PEPSOptimize, GeomSum, ManualIter, LinSolver, EigSolver export fixedpoint diff --git a/src/algorithms/contractions/vumps_contractions.jl b/src/algorithms/contractions/vumps_contractions.jl index 9510dcbe2..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 # @@ -40,6 +38,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, diff --git a/src/algorithms/correlators.jl b/src/algorithms/correlators.jl new file mode 100644 index 000000000..9b5ce5a0e --- /dev/null +++ b/src/algorithms/correlators.jl @@ -0,0 +1,199 @@ +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 + 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 G +end + +function start_correlator( + i::CartesianIndex{2}, + below::InfinitePEPS, + O::MPOTensor, + above::InfinitePEPS, + env::CTMRGEnv, +) + 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}, + js::AbstractVector{CartesianIndex{2}}, + ket::InfinitePEPS, + env::CTMRGEnv, +) + rotated_bra = rotl90(bra) + rotated_ket = bra === ket ? rotated_bra : rotl90(ket) + + 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}, + 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/test/examples/tf_ising.jl b/test/examples/tf_ising.jl index cc6758db5..9ef71b62c 100644 --- a/test/examples/tf_ising.jl +++ b/test/examples/tf_ising.jl @@ -33,12 +33,62 @@ 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 connected correlation functions +corrh = + correlator( + peps, σz ⊗ σz, CartesianIndex(1, 1), CartesianIndex(1, 2):CartesianIndex(1, 21), env + ) .- magnz^2 +corrv = + 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 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 +corrv_int = + 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) # find fixedpoint in polarized phase and compute correlations lengths H_polar = transverse_field_ising(InfiniteSquare(); g=4.5)