From c8419cf72dfb6c53c4601cf589217cff44929cdb Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 5 Aug 2025 08:00:46 -0400 Subject: [PATCH 1/9] Correct hubbard spaces with SU2 particle symmetry --- src/hubbardoperators.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index ad20f1e..c22d1ad 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -56,15 +56,17 @@ function hubbard_space(::Type{U1Irrep}, ::Type{SU2Irrep}) ) end function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial}) - return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) + return Vect[FermionParity ⊠ SU2Irrep]((0, 1 // 2) => 1, (1, 0) => 2) end function hubbard_space(::Type{SU2Irrep}, ::Type{U1Irrep}) return Vect[FermionParity ⊠ SU2Irrep ⊠ U1Irrep]( - (0, 0, 0) => 1, (1, 1 // 2, 1) => 1, (0, 0, 2) => 1 + (0, 1 // 2, 0) => 1, (1, 0, -1 // 2) => 1, (1, 0, 1 // 2) => 1 ) end function hubbard_space(::Type{SU2Irrep}, ::Type{SU2Irrep}) - return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep]((1, 1 // 2, 1 // 2) => 1) + return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep]( + (0, 1 // 2, 0) => 1, (1, 0, 1 // 2) => 1 + ) end # Single-site operators From 9ee4477bc6b18bd0c0048b6e9ac5e3c739b62367 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 5 Aug 2025 08:16:29 -0400 Subject: [PATCH 2/9] Add half_ud_num --- src/hubbardoperators.jl | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index c22d1ad..a6aa9aa 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -3,7 +3,7 @@ module HubbardOperators using TensorKit export hubbard_space -export e_num, u_num, d_num, ud_num +export e_num, u_num, d_num, ud_num, half_ud_num export S_x, S_y, S_z, S_plus, S_min export u_plus_u_min, d_plus_d_min export u_min_u_plus, d_min_d_plus @@ -241,6 +241,28 @@ function ud_num(elt::Type{<:Number}, ::Type{U1Irrep}, ::Type{SU2Irrep}) end const nꜛꜜ = ud_num +@doc """ + half_ud_num([elt::Type{<:Number}], [particle_symmetry::Type{<:Sector}], [spin_symmetry::Type{<:Sector}]) + +Return the one-body operator that counts the number of singly occupied sites. +This is equivalent to `(nꜛ - 1/2)(nꜜ - 1/2)`, and respects the particle-hole symmetry. +""" +half_ud_num(P::Type{<:Sector}, S::Type{<:Sector}) = half_ud_num(ComplexF64, P, S) +function half_ud_num( + elt::Type{<:Number}, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector} + ) + I = id(hubbard_space(particle_symmetry, spin_symmetry)) + return (u_num(elt, particle_symmetry, spin_symmetry) - I / 2) * + (d_num(elt, particle_symmetry, spin_symmetry) - I / 2) +end +function half_ud_num(elt::Type{<:Number}, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + t = single_site_operator(elt, SU2Irrep, SU2Irrep) + block(t, sectortype(t)(0, 1 // 2, 0)) .= 1 // 4 + block(t, sectortype(t)(1, 0, 1 // 2)) .= -1 // 4 + return t +end + @doc """ S_plus(elt::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) S⁺(elt::Type{<:Number}, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) From ca3218b249e920bed5ddbe9e53d3011bebd33d66 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 5 Aug 2025 08:20:52 -0400 Subject: [PATCH 3/9] Add e_hopping for SU2Irrep x SU2Irrep --- src/hubbardoperators.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index a6aa9aa..33130bb 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -610,6 +610,22 @@ function e_hopping( return e_plus_e_min(elt, particle_symmetry, spin_symmetry) - e_min_e_plus(elt, particle_symmetry, spin_symmetry) end +function e_hopping(elt::Type{<:Number}, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + elt <: Complex || throw(DomainError(elt, "SU₂ × SU₂ symmetry requires complex entries")) + t = two_site_operator(elt, SU2Irrep, SU2Irrep) + I = sectortype(t) + even = I(0, 1 // 2, 0) + odd = I(1, 0, 1 // 2) + f1 = only(fusiontrees((odd, odd), one(I))) + f2 = only(fusiontrees((even, even), one(I))) + t[f1, f2] .= 2im + t[f2, f1] .= -2im + f3 = only(fusiontrees((even, odd), I((1, 1 // 2, 1 // 2)))) + f4 = only(fusiontrees((odd, even), I((1, 1 // 2, 1 // 2)))) + t[f3, f4] .= im + t[f4, f3] .= -im + return t +end const e_hop = e_hopping @doc """ From 0f022fbef975bf5d3fb5df245fa5be20134676c3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 5 Aug 2025 09:24:07 -0400 Subject: [PATCH 4/9] Add some tests --- test/hubbardoperators.jl | 37 ++++++++++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 8d214a5..a902d24 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -16,7 +16,7 @@ implemented_symmetries = [ spin_symmetry in [Trivial, U1Irrep, SU2Irrep] if (particle_symmetry, spin_symmetry) in implemented_symmetries - space = hubbard_space(particle_symmetry, spin_symmetry) + space = @inferred hubbard_space(particle_symmetry, spin_symmetry) O = e_plus_e_min(ComplexF64, particle_symmetry, spin_symmetry) O_triv = e_plus_e_min(ComplexF64, Trivial, Trivial) @@ -40,6 +40,14 @@ implemented_symmetries = [ @test_broken S_exchange(ComplexF64, particle_symmetry, spin_symmetry) end end + + # special case for SU2Irrep x SU2Irrep + O = e_hopping(ComplexF64, SU2Irrep, SU2Irrep) + O_triv = e_hopping(ComplexF64, Trivial, Trivial) + test_operator(O, O_triv) + + O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep) + O_triv = half_ud_num(ComplexF64, Trivial, Trivial) end @testset "basic properties" begin @@ -184,6 +192,19 @@ function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) end return H end +function hubbard_hamiltonian(::Type{SU2Irrep}, ::Type{SU2Irrep}; t, U, mu=U / 2, L) + @assert mu ≈ U / 2 + hopping = -t * e_hopping(SU2Irrep, SU2Irrep) + interaction = U * half_ud_num(SU2Irrep, SU2Irrep) + I = id(hubbard_space(SU2Irrep, SU2Irrep)) + H = sum(1:(L - 1)) do i + return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping)) + end + + sum(1:L) do i + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) + end + return H +end @testset "spectrum" begin L = 4 @@ -208,6 +229,20 @@ end sort!(vals_symm) @test vals_triv ≈ vals_symm end + + mu = U / 2 + H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu, L) + vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v) + return repeat(real.(v), dim(c)) .+ U * L / 4 + end + sort!(vals_triv) + + H_symm = hubbard_hamiltonian(SU2Irrep, SU2Irrep; t, U, mu, L) + vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v) + return repeat(real.(v), dim(c)) + end + sort!(vals_symm) + @test vals_triv ≈ vals_symm end @testset "Exact diagonalisation" begin From 3b01d6d1e6d284668c5092ac898399b4124e6641 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 5 Aug 2025 09:30:52 -0400 Subject: [PATCH 5/9] break translation invariance --- test/hubbardoperators.jl | 54 ++++++++++++++++++++++------------------ 1 file changed, 30 insertions(+), 24 deletions(-) diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index a902d24..938d2e1 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -176,43 +176,46 @@ end end end -function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) - hopping = -t * e_hopping(particle_symmetry, spin_symmetry) - interaction = U * ud_num(particle_symmetry, spin_symmetry) - chemical_potential = -mu * e_num(particle_symmetry, spin_symmetry) +function hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu) + L = length(t) + 1 + @assert length(t) + 1 == length(U) == length(mu) + hopping = e_hopping(particle_symmetry, spin_symmetry) + interaction = ud_num(particle_symmetry, spin_symmetry) + chemical_potential = e_num(particle_symmetry, spin_symmetry) I = id(hubbard_space(particle_symmetry, spin_symmetry)) H = sum(1:(L - 1)) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping)) + return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping * -t[i])) end + sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction * U[i])) end + sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential * -mu[i])) end return H end -function hubbard_hamiltonian(::Type{SU2Irrep}, ::Type{SU2Irrep}; t, U, mu=U / 2, L) +function hubbard_hamiltonian(::Type{SU2Irrep}, ::Type{SU2Irrep}; t, U, mu = U ./ 2) + L = length(t) + 1 + @assert length(t) + 1 == length(U) == length(mu) @assert mu ≈ U / 2 - hopping = -t * e_hopping(SU2Irrep, SU2Irrep) - interaction = U * half_ud_num(SU2Irrep, SU2Irrep) + hopping = e_hopping(SU2Irrep, SU2Irrep) + interaction = half_ud_num(SU2Irrep, SU2Irrep) I = id(hubbard_space(SU2Irrep, SU2Irrep)) H = sum(1:(L - 1)) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping)) - end + - sum(1:L) do i - return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) + return reduce(⊗, insert!(collect(Any, fill(I, L - 2)), i, hopping * -t[i])) + end + sum(1:L) do i + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction * U[i])) end return H end @testset "spectrum" begin L = 4 - t = randn() - U = randn() - mu = randn() + t = randn(L - 1) + U = randn(L) + mu = randn(L) - H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu, L) + H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu) vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v) return repeat(real.(v), dim(c)) end @@ -222,7 +225,7 @@ end if (particle_symmetry, spin_symmetry) == (Trivial, Trivial) continue end - H_symm = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) + H_symm = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu) vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v) return repeat(real.(v), dim(c)) end @@ -230,14 +233,14 @@ end @test vals_triv ≈ vals_symm end - mu = U / 2 - H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu, L) + mu = U ./ 2 + H_triv = hubbard_hamiltonian(Trivial, Trivial; t, U, mu) vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v) - return repeat(real.(v), dim(c)) .+ U * L / 4 + return repeat(real.(v), dim(c)) .+ sum(U) / 4 end sort!(vals_triv) - H_symm = hubbard_hamiltonian(SU2Irrep, SU2Irrep; t, U, mu, L) + H_symm = hubbard_hamiltonian(SU2Irrep, SU2Irrep; t, U, mu) vals_symm = mapreduce(vcat, eigvals(H_symm)) do (c, v) return repeat(real.(v), dim(c)) end @@ -257,7 +260,10 @@ end mu = 0.0 E⁻ = U / 2 - sqrt((U / 2)^2 + 4 * t^2) E⁺ = U / 2 + sqrt((U / 2)^2 + 4 * t^2) - H_triv = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) + H_triv = hubbard_hamiltonian( + particle_symmetry, spin_symmetry; + t = fill(t, L - 1), U = fill(U, L), mu = fill(mu, L) + ) # Values based on https://arxiv.org/pdf/0807.4878. Introduction to Hubbard Model and Exact Diagonalization true_eigenvals = sort( From 530086b49b32aae1f453a23c797de61575e52cad Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Wed, 6 Aug 2025 09:14:00 +0800 Subject: [PATCH 6/9] Add missing `test_operator(O, O_triv)` --- test/hubbardoperators.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 938d2e1..47cfa2e 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -48,6 +48,7 @@ implemented_symmetries = [ O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep) O_triv = half_ud_num(ComplexF64, Trivial, Trivial) + test_operator(O, O_triv) end @testset "basic properties" begin From 9b7699f62c4be7a9fb55c8875193ad2f31a37cd8 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 6 Aug 2025 12:15:01 -0400 Subject: [PATCH 7/9] Update docstring --- src/hubbardoperators.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index 33130bb..3954b9d 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -244,8 +244,7 @@ const nꜛꜜ = ud_num @doc """ half_ud_num([elt::Type{<:Number}], [particle_symmetry::Type{<:Sector}], [spin_symmetry::Type{<:Sector}]) -Return the one-body operator that counts the number of singly occupied sites. -This is equivalent to `(nꜛ - 1/2)(nꜜ - 1/2)`, and respects the particle-hole symmetry. +Return the the one-body operator that is equivalent to `(nꜛ - 1/2)(nꜜ - 1/2)`, which respects the particle-hole symmetry. """ half_ud_num(P::Type{<:Sector}, S::Type{<:Sector}) = half_ud_num(ComplexF64, P, S) function half_ud_num( From c94db566ac60ab1627cf632997a76df942fbc53c Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Wed, 6 Aug 2025 12:17:45 -0400 Subject: [PATCH 8/9] add to implemented symmetries --- test/hubbardoperators.jl | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 47cfa2e..5812c5e 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -9,6 +9,7 @@ using StableRNGs implemented_symmetries = [ (Trivial, Trivial), (Trivial, U1Irrep), (Trivial, SU2Irrep), (U1Irrep, Trivial), (U1Irrep, U1Irrep), (U1Irrep, SU2Irrep), + (SU2Irrep, SU2Irrep), ] @testset "Compare symmetric with trivial tensors" begin @@ -18,6 +19,18 @@ implemented_symmetries = [ if (particle_symmetry, spin_symmetry) in implemented_symmetries space = @inferred hubbard_space(particle_symmetry, spin_symmetry) + if particle_symmetry == spin_symmetry == SU2Irrep + O = e_hopping(ComplexF64, SU2Irrep, SU2Irrep) + O_triv = e_hopping(ComplexF64, Trivial, Trivial) + test_operator(O, O_triv) + + O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep) + O_triv = half_ud_num(ComplexF64, Trivial, Trivial) + test_operator(O, O_triv) + else + continue + end + O = e_plus_e_min(ComplexF64, particle_symmetry, spin_symmetry) O_triv = e_plus_e_min(ComplexF64, Trivial, Trivial) test_operator(O, O_triv) @@ -40,15 +53,6 @@ implemented_symmetries = [ @test_broken S_exchange(ComplexF64, particle_symmetry, spin_symmetry) end end - - # special case for SU2Irrep x SU2Irrep - O = e_hopping(ComplexF64, SU2Irrep, SU2Irrep) - O_triv = e_hopping(ComplexF64, Trivial, Trivial) - test_operator(O, O_triv) - - O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep) - O_triv = half_ud_num(ComplexF64, Trivial, Trivial) - test_operator(O, O_triv) end @testset "basic properties" begin @@ -59,6 +63,7 @@ end @test dim(space) == 4 if (particle_symmetry, spin_symmetry) in implemented_symmetries + particle_symmetry == spin_symmetry == SU2Irrep && continue # test hopping operator epem = e_plus_e_min(particle_symmetry, spin_symmetry) emep = e_min_e_plus(particle_symmetry, spin_symmetry) @@ -254,6 +259,7 @@ end spin_symmetry in [Trivial, U1Irrep, SU2Irrep] if (particle_symmetry, spin_symmetry) in implemented_symmetries + particle_symmetry == spin_symmetry == SU2Irrep && continue rng = StableRNG(123) L = 2 From 0c87c9c78eb0b0bd07e6afaf3653d703944788ac Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Tue, 12 Aug 2025 10:25:49 +0200 Subject: [PATCH 9/9] fix some oopsies --- test/hubbardoperators.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 5812c5e..fec0d55 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -27,7 +27,6 @@ implemented_symmetries = [ O = half_ud_num(ComplexF64, SU2Irrep, SU2Irrep) O_triv = half_ud_num(ComplexF64, Trivial, Trivial) test_operator(O, O_triv) - else continue end @@ -228,7 +227,8 @@ end sort!(vals_triv) for (particle_symmetry, spin_symmetry) in implemented_symmetries - if (particle_symmetry, spin_symmetry) == (Trivial, Trivial) + if particle_symmetry == spin_symmetry == Trivial || + particle_symmetry == spin_symmetry == SU2Irrep continue end H_symm = hubbard_hamiltonian(particle_symmetry, spin_symmetry; t, U, mu)