Add missing fermionic parity matrices (#174)

* Fix supertrace to trace in norm

* Add `twistdual`

* Fix fermionic boundarymps contractions

* Add peps sandwich to mpo converter

* Replace trace with planar contraction

* Fix `@tensor` call

* Generalize twist handling in VUMPS contractions a bit more

* Avoid type piracy

* Remove sandwich alias, just hardcode union

* Add test

---------

Co-authored-by: leburgel <lander.burgelman@ugent.be>

Author: lkdvos

src/algorithms/contractions/vumps_contractions.jl

@@ -4,9 +4,29 @@ using MPSKit: GenericMPSTensor, MPSBondTensor
 # Environment transfer functions
 #
 
+function MPSKit.transfer_left(
+    GL::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    A::GenericMPSTensor{S,N},
+    Ā::GenericMPSTensor{S,N},
+) where {S,N}
+    Ā = twistdual(Ā, 2:N)
+    return _transfer_left(GL, O, A, Ā)
+end
+
+function MPSKit.transfer_right(
+    GR::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    A::GenericMPSTensor{S,N},
+    Ā::GenericMPSTensor{S,N},
+) where {S,N}
+    Ā = twistdual(Ā, 2:N)
+    return _transfer_right(GR, O, A, Ā)
+end
+
 ## PEPS
 
-function MPSKit.transfer_left(
+function _transfer_left(
     GL::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     A::GenericMPSTensor{S,3},
@@ -20,7 +40,7 @@ function MPSKit.transfer_left(
         A[χ_NW D_N_above D_N_below; χ_NE]
 end
 
-function MPSKit.transfer_right(
+function _transfer_right(
     GR::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     A::GenericMPSTensor{S,3},
@@ -36,7 +56,7 @@ end
 
 ## PEPO
 
-@generated function MPSKit.transfer_left(
+@generated function _transfer_left(
     GL::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     A::GenericMPSTensor{S,N},
@@ -65,7 +85,7 @@ end
     return macroexpand(@__MODULE__, :(return @autoopt @tensor $GL´_e := $rhs))
 end
 
-@generated function MPSKit.transfer_right(
+@generated function _transfer_right(
     GR::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     A::GenericMPSTensor{S,N},
@@ -118,7 +138,14 @@ end
 # Derivative contractions
 #
 
-@generated function MPSKit.∂C(
+function MPSKit.∂C(
+    C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N}
+) where {S,N}
+    GL = twistdual(GL, 1)
+    GR = twistdual(GR, numind(GR))
+    return _∂C(C, GL, GR)
+end
+@generated function _∂C(
     C::MPSBondTensor{S}, GL::GenericMPSTensor{S,N}, GR::GenericMPSTensor{S,N}
 ) where {S,N}
     C´_e = tensorexpr(:C´, -1, -2)
@@ -128,9 +155,20 @@ end
     return macroexpand(@__MODULE__, :(return @tensor $C´_e := $GL_e * $C_e * $GR_e))
 end
 
+function MPSKit.∂AC(
+    AC::GenericMPSTensor{S,N},
+    O::Union{PEPSSandwich,PEPOSandwich},
+    GL::GenericMPSTensor{S,N},
+    GR::GenericMPSTensor{S,N},
+) where {S,N}
+    GL = twistdual(GL, 1)
+    GR = twistdual(GR, numind(GR))
+    return _∂AC(AC, O, GL, GR)
+end
+
 ## PEPS
 
-function MPSKit.∂AC(
+function _∂AC(
     AC::GenericMPSTensor{S,3},
     O::PEPSSandwich,
     GL::GenericMPSTensor{S,3},
@@ -162,7 +200,7 @@ end
 
 ## PEPO
 
-@generated function MPSKit.∂AC(
+@generated function _∂AC(
     AC::GenericMPSTensor{S,N},
     O::PEPOSandwich{H},
     GL::GenericMPSTensor{S,N},

src/algorithms/toolbox.jl

@@ -142,12 +142,11 @@ environment `env`.
 """
 function _contract_corners(ind::Tuple{Int,Int}, env::CTMRGEnv)
     r, c = ind
-    return tr(
-        env.corners[NORTHWEST, _prev(r, end), _prev(c, end)] *
-        env.corners[NORTHEAST, _prev(r, end), c] *
-        env.corners[SOUTHEAST, r, c] *
-        env.corners[SOUTHWEST, r, _prev(c, end)],
-    )
+    C_NW = env.corners[NORTHWEST, _prev(r, end), _prev(c, end)]
+    C_NE = env.corners[NORTHEAST, _prev(r, end), c]
+    C_SE = env.corners[SOUTHEAST, r, c]
+    C_SW = env.corners[SOUTHWEST, r, _prev(c, end)]
+    return @tensor C_NW[1; 2] * C_NE[2; 3] * C_SE[3; 4] * C_SW[4; 1]
 end
 
 """

src/networks/local_sandwich.jl

@@ -55,6 +55,39 @@ function virtualspace(O::PEPSSandwich, dir)
     return virtualspace(ket(O), dir) ⊗ virtualspace(bra(O), dir)'
 end
 
+# not overloading MPOTensor because that defines AbstractTensorMap{<:Any,S,2,2}(::PEPSTensor, ::PEPSTensor)
+# ie type piracy
+mpotensor(top::PEPSTensor) = mpotensor((top, top))
+function mpotensor((top, bot)::PEPSSandwich)
+    @assert virtualspace(top, NORTH) == dual(virtualspace(top, SOUTH)) &&
+        virtualspace(bot, NORTH) == dual(virtualspace(bot, SOUTH)) &&
+        virtualspace(top, EAST) == dual(virtualspace(top, WEST)) &&
+        virtualspace(bot, EAST) == dual(virtualspace(bot, WEST)) &&
+        isdual(virtualspace(top, NORTH)) &&
+        isdual(virtualspace(bot, NORTH)) &&
+        isdual(virtualspace(top, EAST)) &&
+        isdual(virtualspace(bot, EAST)) "Method not yet implemented for given virtual spaces"
+
+    F_west = isomorphism(
+        storagetype(top),
+        fuse(virtualspace(top, WEST), virtualspace(bot, WEST)'),
+        virtualspace(top, WEST) ⊗ virtualspace(bot, WEST)',
+    )
+    F_south = isomorphism(
+        storagetype(top),
+        fuse(virtualspace(top, SOUTH), virtualspace(bot, SOUTH)'),
+        virtualspace(top, SOUTH) ⊗ virtualspace(bot, SOUTH)',
+    )
+    @tensor O[west south; north east] :=
+        top[phys; top_north top_east top_south top_west] *
+        conj(bot[phys; bot_north bot_east bot_south bot_west]) *
+        twist(F_west, 3)[west; top_west bot_west] *
+        twist(F_south, 3)[south; top_south bot_south] *
+        conj(F_west[east; top_east bot_east]) *
+        conj(F_south[north; top_north bot_north])
+    return O
+end
+
 ## PEPO
 
 const PEPOSandwich{N,T<:PEPSTensor,P<:PEPOTensor} = Tuple{T,T,Vararg{P,N}}

src/utility/util.jl

@@ -102,6 +102,22 @@ function absorb_s(u::AbstractTensorMap, s::DiagonalTensorMap, vh::AbstractTensor
     return u * sqrt_s, sqrt_s * vh
 end
 
+"""
+    twistdual(t::AbstractTensorMap, i)
+    twistdual!(t::AbstractTensorMap, i)
+
+Twist the i-th leg of a tensor `t` if it represents a dual space.
+"""
+function twistdual!(t::AbstractTensorMap, i::Int)
+    isdual(space(t, i)) || return t
+    return twist!(t, i)
+end
+function twistdual!(t::AbstractTensorMap, is)
+    is′ = filter(i -> isdual(space(t, i)), is)
+    return twist!(t, is′)
+end
+twistdual(t::AbstractTensorMap, is) = twistdual!(copy(t), is)
+
 # Check whether diagonals contain degenerate values up to absolute or relative tolerance
 function is_degenerate_spectrum(
     S; atol::Real=0, rtol::Real=atol > 0 ? 0 : sqrt(eps(scalartype(S)))

test/boundarymps/vumps.jl

@@ -8,6 +8,7 @@ using LinearAlgebra
 Random.seed!(29384293742893)
 
 const vumps_alg = VUMPS(; alg_eigsolve=MPSKit.Defaults.alg_eigsolve(; ishermitian=false))
+
 @testset "(1, 1) PEPS" begin
     psi = InfinitePEPS(ComplexSpace(2), ComplexSpace(2))
 
@@ -37,6 +38,41 @@ end
     @test N ≈ N´ rtol = 1e-2
 end
 
+@testset "Fermionic PEPS" begin
+    D = Vect[fℤ₂](0 => 1, 1 => 1)
+    d = Vect[fℤ₂](0 => 1, 1 => 1)
+    χ = Vect[fℤ₂](0 => 10, 1 => 10)
+
+    psi = InfinitePEPS(D, d; unitcell=(1, 1))
+    n = InfiniteSquareNetwork(psi)
+    T = PEPSKit.InfiniteTransferPEPS(psi, 1, 1)
+
+    # compare boundary MPS contraction to CTMRG contraction
+    mps = PEPSKit.initializeMPS(T, [χ])
+    mps, env, ϵ = leading_boundary(mps, T, vumps_alg)
+    N_vumps = abs(prod(expectation_value(mps, T)))
+
+    ctm, = leading_boundary(CTMRGEnv(psi, χ), psi)
+    N_ctm = abs(norm(psi, ctm))
+
+    @test N_vumps ≈ N_ctm rtol = 1e-2
+
+    # and again after blocking the local sandwiches
+    n´ = InfiniteSquareNetwork(map(PEPSKit.mpotensor, PEPSKit.unitcell(n)))
+    T´ = InfiniteMPO(map(PEPSKit.mpotensor, T.O))
+
+    mps´ = InfiniteMPS(randn, ComplexF64, [physicalspace(T´, 1)], [χ])
+    mps´, env´, ϵ = leading_boundary(mps´, T´, vumps_alg)
+    N_vumps´ = abs(prod(expectation_value(mps´, T´)))
+
+    ctm´, = leading_boundary(CTMRGEnv(n´, χ), n´)
+    N_ctm´ = abs(network_value(n´, ctm´))
+
+    @show N_vumps´
+    @test N_vumps´ ≈ N_vumps rtol = 1e-2
+    @test N_vumps´ ≈ N_ctm´ rtol = 1e-2
+end
+
 @testset "PEPO runthrough" begin
     function ising_pepo(beta; unitcell=(1, 1, 1))
         t = ComplexF64[exp(beta) exp(-beta); exp(-beta) exp(beta)]