Add tests for fermionic periodic boundary conditions (#274)

* Fix and tests fermionic PBC

* Refactor ambiguous braiding checks

* Circumvent unsupported tensorkit feature

* formatting

* Fixes attempt II

* Fix tests

* remove twists

* Clean up utility

Author: lkdvos

src/transfermatrix/transfer.jl

@@ -68,25 +68,21 @@ end
 # the transfer operation of a density matrix with a utility leg in its codomain is ill defined - how should one braid the utility leg?
 # hence the checks - to make sure that this operation is uniquely defined
 function transfer_left(v::MPSTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S}
-    _can_unambiguously_braid(space(v, 2)) ||
-        throw(ArgumentError("transfer is not uniquely defined with utility space $(space(v,2))"))
+    check_unambiguous_braiding(space(v, 2))
     @plansor v[-1 -2; -3] := v[1 2; 4] * A[4 5; -3] * τ[2 3; 5 -2] * conj(Ab[1 3; -1])
 end
 function transfer_right(v::MPSTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S}
-    _can_unambiguously_braid(space(v, 2)) ||
-        throw(ArgumentError("transfer is not uniquely defined with utility space $(space(v,2))"))
+    check_unambiguous_braiding(space(v, 2))
     @plansor v[-1 -2; -3] := A[-1 2; 1] * τ[-2 4; 2 3] * conj(Ab[-3 4; 5]) * v[1 3; 5]
 end
 
 # the transfer operation with a utility leg in both the domain and codomain is also ill defined - only due to the codomain utility space
 function transfer_left(v::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S}
-    _can_unambiguously_braid(space(v, 2)) ||
-        throw(ArgumentError("transfer is not uniquely defined with utility space $(space(v,2))"))
+    check_unambiguous_braiding(space(v, 2))
     @plansor t[-1 -2; -3 -4] := v[1 2; -3 4] * A[4 5; -4] * τ[2 3; 5 -2] * conj(Ab[1 3; -1])
 end
 function transfer_right(v::MPOTensor{S}, A::MPSTensor{S}, Ab::MPSTensor{S}) where {S}
-    _can_unambiguously_braid(space(v, 2)) ||
-        throw(ArgumentError("transfer is not uniquely defined with utility space $(space(v,2))"))
+    check_unambiguous_braiding(space(v, 2))
     @plansor t[-1 -2; -3 -4] := A[-1 2; 1] * τ[-2 4; 2 3] * conj(Ab[-4 4; 5]) * v[1 3; -3 5]
 end
 

src/utility/utility.jl

@@ -89,18 +89,6 @@ function _embedders(spaces)
     return maps
 end
 
-function _can_unambiguously_braid(sp::VectorSpace)
-    s = sectortype(sp)
-
-    BraidingStyle(s) isa SymmetricBraiding && return true
-
-    # if it's not symmetric, then we are only really garantueed that this is possible when only one irrep occurs - the trivial one
-    for sect in sectors(sp)
-        sect == one(sect) || return false
-    end
-    return true
-end
-
 #=
 map every element in the tensormap to dfun(E)
 allows us to create random tensormaps for any storagetype
@@ -143,3 +131,21 @@ end
 function fuser(::Type{T}, V1::S, V2::S) where {T,S<:IndexSpace}
     return isomorphism(T, fuse(V1 ⊗ V2), V1 ⊗ V2)
 end
+
+"""
+    check_unambiguous_braiding(::Type{Bool}, V::VectorSpace)::Bool
+    check_unambiguous_braiding(V::VectorSpace)
+
+Verify that the braiding of a vector space is unambiguous. This is the case if the braiding
+is symmetric or if all sectors are trivial. The signature with `Type{Bool}` is used to check
+while the signature without is used to throw an error if the braiding is ambiguous.
+"""
+function check_unambiguous_braiding(::Type{Bool}, V::VectorSpace)
+    I = sectortype(V)
+    BraidingStyle(I) isa SymmetricBraiding && return true
+    return all(isone, sectors(V))
+end
+function check_unambiguous_braiding(V::VectorSpace)
+    return check_unambiguous_braiding(Bool, V) ||
+           throw(ArgumentError("cannot unambiguously braid $V"))
+end

test/algorithms.jl

@@ -766,6 +766,17 @@ end
                   permute(TH, ((vcat(N, 1:(N - 1))...,), (vcat(2N, (N + 1):(2N - 1))...,)))
         end
     end
+
+    # fermionic tests
+    for N in 3:5
+        h = real(c_plusmin() + c_minplus())
+        H = InfiniteMPOHamiltonian([space(h, 1)], (1, 2) => h)
+        H_periodic = periodic_boundary_conditions(H, N)
+        terms = [(i, i + 1) => h for i in 1:(N - 1)]
+        push!(terms, (1, N) => h)
+        H_periodic2 = FiniteMPOHamiltonian(physicalspace(H_periodic), terms)
+        @test H_periodic ≈ H_periodic2
+    end
 end
 
 @testset "TaylorCluster time evolution" begin

test/setup.jl

@@ -13,9 +13,11 @@ using Combinatorics: permutations
 
 # exports
 export S_xx, S_yy, S_zz, S_x, S_y, S_z
+export c_plusmin, c_minplus, c_number
 export force_planar
 export symm_mul_mpo
-export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model
+export transverse_field_ising, heisenberg_XXX, bilinear_biquadratic_model, XY_model,
+       kitaev_model
 export classical_ising, finite_classical_ising, sixvertex
 
 # using TensorOperations
@@ -224,6 +226,61 @@ function bilinear_biquadratic_model(::Type{SU2Irrep}; θ=atan(1 / 3), L=Inf)
     end
 end
 
+function c_plusmin()
+    P = Vect[FermionParity](0 => 1, 1 => 1)
+    t = zeros(ComplexF64, P^2 ← P^2)
+    I = sectortype(P)
+    t[(I(1), I(0), dual(I(0)), dual(I(1)))] .= 1
+    return t
+end
+
+function c_minplus()
+    P = Vect[FermionParity](0 => 1, 1 => 1)
+    t = zeros(ComplexF64, P^2 ← P^2)
+    I = sectortype(t)
+    t[(I(0), I(1), dual(I(1)), dual(I(0)))] .= 1
+    return t
+end
+
+function c_plusplus()
+    P = Vect[FermionParity](0 => 1, 1 => 1)
+    t = zeros(ComplexF64, P^2 ← P^2)
+    I = sectortype(t)
+    t[(I(1), I(1), dual(I(0)), dual(I(0)))] .= 1
+    return t
+end
+
+function c_minmin()
+    P = Vect[FermionParity](0 => 1, 1 => 1)
+    t = zeros(ComplexF64, P^2 ← P^2)
+    I = sectortype(t)
+    t[(I(0), I(0), dual(I(1)), dual(I(1)))] .= 1
+    return t
+end
+
+function c_number()
+    P = Vect[FermionParity](0 => 1, 1 => 1)
+    t = zeros(ComplexF64, P ← P)
+    block(t, fℤ₂(1)) .= 1
+    return t
+end
+
+function kitaev_model(; t=1.0, mu=1.0, Delta=1.0, L=Inf)
+    TB = scale!(c_plusmin() + c_minplus(), -t / 2)     # tight-binding term
+    SC = scale!(c_plusplus() + c_minmin(), Delta / 2)  # superconducting term
+    CP = scale!(c_number(), -mu)                       # chemical potential term
+
+    if L == Inf
+        lattice = PeriodicArray([space(TB, 1)])
+        return InfiniteMPOHamiltonian(lattice, (1, 2) => TB + SC, (1,) => CP)
+    else
+        lattice = fill(space(TB, 1), L)
+        terms = Iterators.flatten((((i, i + 1) => TB + SC for i in 1:(L - 1)),
+                                   (i,) => CP for i in 1:L))
+        return FiniteMPOHamiltonian(lattice, terms)
+    end
+end
+
 function ising_bond_tensor(β)
     t = [exp(β) exp(-β); exp(-β) exp(β)]
     r = eigen(t)