Add correlation functions (#210)

---------

Co-authored-by: Lukas Devos <ldevos98@gmail.com>

Author: sanderdemeyer

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

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,

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

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)