Fix and tests fermionic PBC

Author: lkdvos

src/algorithms/toolbox.jl

@@ -274,9 +274,11 @@ function periodic_boundary_conditions(mpo::InfiniteMPO{O},
         F_left = i == 1 ? cup : F_right
         F_right = i == L ? cup :
                   isomorphism(ST, V_right ← V_wrap' * right_virtualspace(mpo, i))
+        # account for fermionic boundary conditions
+        o = i == L ? twist(mpo[i], 4) : mpo[i]
         @plansor output[i][-1 -2; -3 -4] = F_left[-1; 1 2] *
                                            τ[-3 1; 4 3] *
-                                           mpo[i][2 -2; 3 5] *
+                                           o[2 -2; 3 5] *
                                            conj(F_right[-4; 4 5])
     end
 
@@ -421,7 +423,9 @@ function periodic_boundary_conditions(H::InfiniteMPOHamiltonian, L=length(H))
             F_left = fusers[end][j, k, chi]
             k′ = indmap[j, k, chi]
             k′ == size(output[end], 1) && continue
-            @plansor o[-1 -2; -3 -4] := F_left[-1; 1 2 6] * h[1 3; -3 4] * τ[3 2; 4 5] *
+            # apply twist for fermionic boundary conditions
+            h′ = twist(h, 4)
+            @plansor o[-1 -2; -3 -4] := F_left[-1; 1 2 6] * h′[1 3; -3 4] * τ[3 2; 4 5] *
                                         τ[5 6; -4 -2]
             output[end][k′, 1, 1, 1] = o
         end

test/algorithms.jl

@@ -759,13 +759,24 @@ end
 
 @testset "periodic boundary conditions" begin
     Hs = [transverse_field_ising(), heisenberg_XXX(), classical_ising(), sixvertex()]
-    for N in 2:6
+    for N in 2:5
         for H in Hs
             TH = convert(TensorMap, periodic_boundary_conditions(H, N))
             @test TH ≈
                   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)