Fix expectation values for PEPS/PEPO with dual physical spaces

src/operators/infinitepepo.jl

@@ -141,9 +141,18 @@ end
 ## Spaces
 
 TensorKit.spacetype(::Type{P}) where {P <: InfinitePEPO} = spacetype(eltype(P))
-virtualspace(T::InfinitePEPO, r::Int, c::Int, h::Int, dir) = virtualspace(T[r, c, h], dir)
-domain_physicalspace(T::InfinitePEPO, r::Int, c::Int) = domain_physicalspace(T[r, c, 1])
-codomain_physicalspace(T::InfinitePEPO, r::Int, c::Int) = codomain_physicalspace(T[r, c, end])
+function virtualspace(T::InfinitePEPO, r::Int, c::Int, h::Int, dir)
+    Nr, Nc, Nh = size(T)
+    return virtualspace(T[mod1(r, Nr), mod1(c, Nc), mod1(h, Nh)], dir)
+end
+function domain_physicalspace(T::InfinitePEPO, r::Int, c::Int)
+    Nr, Nc, = size(T)
+    return domain_physicalspace(T[mod1(r, Nr), mod1(c, Nc), 1])
+end
+function codomain_physicalspace(T::InfinitePEPO, r::Int, c::Int)
+    Nr, Nc, = size(T)
+    return codomain_physicalspace(T[mod1(r, Nr), mod1(c, Nc), end])
+end
 physicalspace(T::InfinitePEPO) = [physicalspace(T, row, col) for row in axes(T, 1), col in axes(T, 2)]
 function physicalspace(T::InfinitePEPO, r::Int, c::Int)
     codomain_physicalspace(T, r, c) == domain_physicalspace(T, r, c) || throw(

src/states/infinitepartitionfunction.jl

@@ -139,7 +139,10 @@ end
 ## Spaces
 
 TensorKit.spacetype(::Type{T}) where {T <: InfinitePartitionFunction} = spacetype(eltype(T))
-virtualspace(n::InfinitePartitionFunction, r::Int, c::Int, dir) = virtualspace(n[r, c], dir)
+function virtualspace(n::InfinitePartitionFunction, r::Int, c::Int, dir)
+    Nr, Nc = size(n)
+    return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir)
+end
 
 ## InfiniteSquareNetwork interface
 

src/states/infinitepeps.jl

@@ -140,9 +140,15 @@ end
 
 TensorKit.spacetype(::Type{T}) where {T <: InfinitePEPS} = spacetype(eltype(T))
 virtualspace(n::InfinitePEPS, dir) = virtualspace.(unitcell(n), dir)
-virtualspace(n::InfinitePEPS, r::Int, c::Int, dir) = virtualspace(n[r, c], dir)
+function virtualspace(n::InfinitePEPS, r::Int, c::Int, dir)
+    Nr, Nc = size(n)
+    return virtualspace(n[mod1(r, Nr), mod1(c, Nc)], dir)
+end
 physicalspace(n::InfinitePEPS) = physicalspace.(unitcell(n))
-physicalspace(n::InfinitePEPS, r::Int, c::Int) = physicalspace(n[r, c])
+function physicalspace(n::InfinitePEPS, r::Int, c::Int)
+    Nr, Nc = size(n)
+    return physicalspace(n[mod1(r, Nr), mod1(c, Nc)])
+end
 
 ## InfiniteSquareNetwork interface
 

test/runtests.jl

@@ -72,6 +72,9 @@ end
         @time @safetestset "Density matrix from double-layer PEPO" begin
             include("toolbox/densitymatrices.jl")
         end
+        @time @safetestset "Expectation value of identity operator" begin
+            include("toolbox/expval_id.jl")
+        end
     end
     if GROUP == "ALL" || GROUP == "UTILITY"
         @time @safetestset "LocalOperator" begin

test/toolbox/expval_id.jl

@@ -0,0 +1,98 @@
+using Test
+using TensorKit
+using PEPSKit
+using Random
+
+const syms = (FermionParity, U1Irrep)
+Random.seed!(0)
+
+function get_spaces(sym::Type{<:Sector})
+    if sym == U1Irrep
+        Vp = Vect[U1Irrep](-1 => 1, 0 => 1, 1 => 2)
+        Vv = Vect[U1Irrep](-1 => 1, 0 => 3, 1 => 2)
+        Ve = Vect[U1Irrep](-1 => 1, 0 => 3, 1 => 2)
+        return Vp, Vv, Ve
+    elseif sym == FermionParity
+        Vp = Vect[FermionParity](0 => 1, 1 => 2)
+        Vv = Vect[FermionParity](0 => 3, 1 => 2)
+        Ve = Vect[FermionParity](0 => 3, 1 => 2)
+        return Vp, Vv, Ve
+    else
+        error("Got a sector not intended to be tested.")
+    end
+end
+
+@testset "⟨ψ|1|ψ⟩ (InfinitePEPS, $(sym))" for sym in syms
+    Vp, Vv, Ve = get_spaces(sym)
+    Pspaces = [Vp Vp'; Vp Vp']
+    Vspaces = fill(Vv, (2, 2))
+    ψ = InfinitePEPS(Pspaces, Vspaces)
+    env = CTMRGEnv(ψ, Ve)
+    for site1 in Tuple.(CartesianIndices((2, 2)))
+        # 1-site
+        id1 = TensorKit.id(physicalspace(ψ, site1...))
+        O1 = LocalOperator(Pspaces, (site1,) => id1)
+        val1 = expectation_value(ψ, O1, env)
+        @info "$((site1,)): $(val1)"
+        @test val1 ≈ 1
+        # 2-site
+        for d in [(1, 0), (0, 1), (1, 1), (-1, 1)]
+            site2 = site1 .+ d
+            id2 = TensorKit.id(physicalspace(ψ, site2...))
+            O2 = LocalOperator(Pspaces, (site1, site2) => id1 ⊗ id2)
+            val2 = expectation_value(ψ, O2, env)
+            @info "$((site1, site2)): $(val2)"
+            @test val2 ≈ 1
+        end
+    end
+end
+
+@testset "tr(ρ1) (one-layer InfinitePEPO, $(sym))" for sym in syms
+    Vp, Vv, Ve = get_spaces(sym)
+    Pspaces = [Vp Vp'; Vp Vp']
+    Vspaces = fill(Vv, (2, 2))
+    ρ = InfinitePEPO(Pspaces, Vspaces)
+    env = CTMRGEnv(InfinitePartitionFunction(ρ), Ve)
+    for site1 in Tuple.(CartesianIndices((2, 2)))
+        # 1-site
+        id1 = TensorKit.id(physicalspace(ρ, site1...))
+        O1 = LocalOperator(Pspaces, (site1,) => id1)
+        val1 = expectation_value(ρ, O1, env)
+        @info "$((site1,)): $(val1)"
+        @test val1 ≈ 1
+        # 2-site
+        for d in [(1, 0), (0, 1), (1, 1), (-1, 1)]
+            site2 = site1 .+ d
+            id2 = TensorKit.id(physicalspace(ρ, site2...))
+            O2 = LocalOperator(Pspaces, (site1, site2) => id1 ⊗ id2)
+            val2 = expectation_value(ρ, O2, env)
+            @info "$((site1, site2)): $(val2)"
+            @test val2 ≈ 1
+        end
+    end
+end
+
+@testset "⟨ρ|1|ρ⟩ (InfinitePEPS with ancilla, $(sym))" for sym in syms
+    Vp, Vv, Ve = get_spaces(sym)
+    Pspaces = [Vp Vp'; Vp Vp']
+    Vspaces = fill(Vv, (2, 2))
+    ρ = InfinitePEPO(Pspaces, Vspaces)
+    env = CTMRGEnv(InfinitePEPS(ρ), Ve)
+    for site1 in Tuple.(CartesianIndices((2, 2)))
+        # 1-site
+        id1 = TensorKit.id(physicalspace(ρ, site1...))
+        O1 = LocalOperator(Pspaces, (site1,) => id1)
+        val1 = expectation_value(ρ, O1, ρ, env)
+        @info "$((site1,)): $(val1)"
+        @test val1 ≈ 1
+        # 2-site
+        for d in [(1, 0), (0, 1), (1, 1), (-1, 1)]
+            site2 = site1 .+ d
+            id2 = TensorKit.id(physicalspace(ρ, site2...))
+            O2 = LocalOperator(Pspaces, (site1, site2) => id1 ⊗ id2)
+            val2 = expectation_value(ρ, O2, ρ, env)
+            @info "$((site1, site2)): $(val2)"
+            @test val2 ≈ 1
+        end
+    end
+end