From 5d8ef2836426d8a321e01bf17d9367f48e196834 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 12:44:28 -0500 Subject: [PATCH 1/6] Refactor and correct hubbard operators --- src/MPSKitModels.jl | 2 + src/operators/fermionoperators.jl | 163 --------------- src/operators/hubbardoperators.jl | 336 ++++++++++++++++++++++++++++++ src/utility.jl | 13 ++ 4 files changed, 351 insertions(+), 163 deletions(-) create mode 100644 src/operators/hubbardoperators.jl diff --git a/src/MPSKitModels.jl b/src/MPSKitModels.jl index e20926f..7582214 100644 --- a/src/MPSKitModels.jl +++ b/src/MPSKitModels.jl @@ -63,6 +63,8 @@ include("operators/mpoham.jl") include("operators/spinoperators.jl") include("operators/fermionoperators.jl") +include("operators/hubbardoperators.jl") +using .HubbardOperators include("operators/bosonoperators.jl") include("models/hamiltonians.jl") diff --git a/src/operators/fermionoperators.jl b/src/operators/fermionoperators.jl index 67ee7c3..ece6f98 100644 --- a/src/operators/fermionoperators.jl +++ b/src/operators/fermionoperators.jl @@ -66,166 +66,3 @@ function c_number(elt::Type{<:Number}=ComplexF64) blocks(n)[fℤ₂(1)] .= one(elt) return n end - -#=========================================================================================== - spin 1/2 fermions -===========================================================================================# - -""" - e_plus([elt::Type{<:Number}=ComplexF64], particle_symmetry, spin_symmetry; side=:L) - e⁺([elt::Type{<:Number}=ComplexF64], particle_symmetry, spin_symmetry; side=:L) - -The creation operator for electron-like fermions. -""" -function e_plus end -function e_plus(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; kwargs...) - return e_plus(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) -end -function e_plus(elt::Type{<:Number}=ComplexF64, - ::Type{Trivial}=Trivial, - ::Type{Trivial}=Trivial; - side=:L) - pspace = Vect[fℤ₂](0 => 2, 1 => 2) - vspace = Vect[fℤ₂](1 => 2) - if side == :L - e⁺ = TensorMap(zeros, elt, pspace ← pspace ⊗ vspace) - blocks(e⁺)[fℤ₂(0)][2, 2:3] .= [one(elt), -one(elt)] - blocks(e⁺)[fℤ₂(1)][:, 1:2] .= [one(elt) zero(elt); zero(elt) one(elt)] - elseif side == :R - E = e_plus(elt, Trivial, Trivial; side=:L) - F = isomorphism(storagetype(E), vspace, flip(vspace)) - @planar e⁺[-1 -2; -3] := E[-2; 1 2] * τ[1 2; 3 -3] * F[3; -1] - else - throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) - end - return e⁺ -end -function e_plus(elt::Type{<:Number}, ::Type{U1Irrep}, ::Type{SU2Irrep}; side=:L) - pspace = Vect[fℤ₂ ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (0, 2, 0) => 1) - vspace = Vect[fℤ₂ ⊠ U1Irrep ⊠ SU2Irrep]((1, 1, 1 // 2) => 1) - if side == :L - e⁺ = TensorMap(zeros, elt, pspace ← pspace ⊗ vspace) - blocks(e⁺)[fℤ₂(0) ⊠ U1Irrep(2) ⊠ SU2Irrep(0)] .= sqrt(2) - blocks(e⁺)[fℤ₂(1) ⊠ U1Irrep(1) ⊠ SU2Irrep(1 // 2)] .= 1 - elseif side == :R - E = e_plus(elt, U1Irrep, SU2Irrep; side=:L) - F = isomorphism(storagetype(E), vspace, flip(vspace)) - @planar e⁺[-1 -2; -3] := E[-2; 1 2] * τ[1 2; 3 -3] * F[3; -1] - end - return e⁺ -end -const e⁺ = e_plus - -""" - e_min([elt::Type{<:Number}=ComplexF64], particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; side=:L) - e⁻([elt::Type{<:Number}=ComplexF64], particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; side=:L) - -The annihilation operator for electron-like fermions. -""" -function e_min end -function e_min(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; kwargs...) - return e_min(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) -end -function e_min(elt::Type{<:Number}=ComplexF64, - particle_symmetry::Type{<:Sector}=Trivial, - spin_symmetry::Type{<:Sector}=Trivial; - side=:L) - if side === :L - E = e_plus(elt, particle_symmetry, spin_symmetry; side=:L)' - F = isomorphism(storagetype(E), flip(space(E, 2)), space(E, 2)) - @planar e⁻[-1; -2 -3] := E[-1 1; -2] * F[-3; 1] - elseif side === :R - e⁻ = permute(e_plus(elt, particle_symmetry, spin_symmetry; side=:L)', - ((2, 1), (3,))) - else - throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) - end - return e⁻ -end -const e⁻ = e_min - -function e_plusmin(elt::Type{<:Number}=ComplexF64, - particle_symmetry::Type{<:Sector}=Trivial, - spin_symmetry::Type{<:Sector}=Trivial) - return contract_twosite(e⁺(elt, particle_symmetry, spin_symmetry; side=:L), - e⁻(elt, particle_symmetry, spin_symmetry; side=:R)) -end -function e_plusmin(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) - return e_plusmin(ComplexF64, particle_symmetry, spin_symmetry) -end -const e⁺e⁻ = e_plusmin - -function e_minplus(elt::Type{<:Number}=ComplexF64, - particle_symmetry::Type{<:Sector}=Trivial, - spin_symmetry::Type{<:Sector}=Trivial) - return contract_twosite(e⁻(elt, particle_symmetry, spin_symmetry; side=:L), - e⁺(elt, particle_symmetry, spin_symmetry; side=:R)) -end -function e_minplus(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) - return e_minplus(ComplexF64, particle_symmetry, spin_symmetry) -end -const e⁻e⁺ = e_minplus - -""" - e_number([elt::Type{<:Number}=ComplexF64], particle_symmetry=fℤ₂, spin_symmetry=ℤ₁) - -The number operator for electron-like fermions. -""" -function e_number end -function e_number(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) - return e_number(ComplexF64, particle_symmetry, spin_symmetry) -end -function e_number(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial, - ::Type{Trivial}=Trivial) - pspace = Vect[fℤ₂](0 => 2, 1 => 2) - n = TensorMap(zeros, elt, pspace ← pspace) - blocks(n)[fℤ₂(1)][1, 1] = 1 - blocks(n)[fℤ₂(1)][2, 2] = 1 - blocks(n)[fℤ₂(0)][2, 2] = 2 - return n -end -function e_number(elt::Type{<:Number}, ::Type{U1Irrep}, ::Type{SU2Irrep}) - pspace = Vect[fℤ₂ ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (0, 2, 0) => 1) - n = TensorMap(zeros, elt, pspace ← pspace) - for (c, b) in blocks(n) - b .= c.sectors[2].charge - end - return n -end - -function e_number_up(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial, - ::Type{Trivial}=Trivial) - pspace = Vect[fℤ₂](0 => 2, 1 => 2) - n = TensorMap(zeros, elt, pspace ← pspace) - blocks(n)[fℤ₂(1)][1, 1] = 1 - blocks(n)[fℤ₂(0)][2, 2] = 1 - return n -end - -function e_number_down(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial, - ::Type{Trivial}=Trivial) - pspace = Vect[fℤ₂](0 => 2, 1 => 2) - n = TensorMap(zeros, elt, pspace ← pspace) - blocks(n)[fℤ₂(1)][2, 2] = 1 - blocks(n)[fℤ₂(0)][2, 2] = 1 - return n -end - -function e_number_updown(elt::Type{<:Number}=ComplexF64, ::Type{Trivial}=Trivial, - ::Type{Trivial}=Trivial) - pspace = Vect[fℤ₂](0 => 2, 1 => 2) - n = TensorMap(zeros, elt, pspace ← pspace) - blocks(n)[fℤ₂(0)][2, 2] = 1 - return n -end -function e_number_updown(elt::Type{<:Number}, ::Type{U1Irrep}, ::Type{SU2Irrep}) - pspace = Vect[fℤ₂ ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, - (0, 2, 0) => 1) - n = TensorMap(zeros, elt, pspace ← pspace) - blocks(n)[fℤ₂(0) ⊠ U1Irrep(2) ⊠ SU2Irrep(0)] .= 1 - return n -end - -const nꜛnꜜ = e_number_updown diff --git a/src/operators/hubbardoperators.jl b/src/operators/hubbardoperators.jl new file mode 100644 index 0000000..f9360e2 --- /dev/null +++ b/src/operators/hubbardoperators.jl @@ -0,0 +1,336 @@ +# Operators that act on Hubbard-type models +# i.e. the local hilbert space consists of |∅⟩, |↑⟩, |↓⟩, |↑↓⟩ +module HubbardOperators + +using TensorKit + +export hubbard_space +export e_plusmin, e_plusmin_up, e_plusmin_down +export e_minplus, e_minplus_up, e_minplus_down +export e_number, e_number_up, e_number_down, e_number_updown + +export e⁺e⁻, e⁺e⁻ꜛ, e⁺e⁻ꜜ, e⁻e⁺, e⁻e⁺ꜛ, e⁻e⁺ꜜ +export nꜛ, nꜜ, nꜛnꜜ +# not exported because namespace: export n + +""" + hubbard_space(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the local hilbert space for a Hubbard-type model with the given particle and spin symmetries. +The possible symmetries are `Trivial`, `U1Irrep`, and `SU2Irrep`, for both particle number and spin. +""" +function hubbard_space(::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) + return Vect[FermionParity](0 => 2, 1 => 2) +end +function hubbard_space(::Type{Trivial}, ::Type{U1Irrep}) + return Vect[FermionParity ⊠ U1Irrep]((0, 0) => 2, (1, 1 // 2) => 1, (1, -1 // 2) => 1) +end +function hubbard_space(::Type{Trivial}, ::Type{SU2Irrep}) + return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) +end +function hubbard_space(::Type{U1Irrep}, ::Type{Trivial}=Trivial) + return Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2, (0, 2) => 1) +end +function hubbard_space(::Type{U1Irrep}, ::Type{U1Irrep}) + return Vect[FermionParity ⊠ U1Irrep ⊠ U1Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, + (1, 1, -1 // 2) => 1, (0, 2, 0) => 1) +end +function hubbard_space(::Type{U1Irrep}, ::Type{SU2Irrep}) + return Vect[FermionParity ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, + (0, 2, 0) => 1) +end +function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial}=Trivial) + return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) +end +function hubbard_space(::Type{SU2Irrep}, ::Type{U1Irrep}) + return Vect[FermionParity ⊠ SU2Irrep ⊠ U1Irrep]((0, 0, 0) => 1, (1, 1 // 2, 1) => 1) +end +function hubbard_space(::Type{SU2Irrep}, ::Type{SU2Irrep}) + return Vect[FermionParity ⊠ SU2Irrep ⊠ SU2Irrep]((1, 1 // 2, 1 // 2) => 1) +end + +function single_site_operator(T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}) + V = hubbard_space(particle_symmetry, spin_symmetry) + return TensorMap(zeros, T, V ← V) +end + +function two_site_operator(T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}) + V = hubbard_space(particle_symmetry, spin_symmetry) + return TensorMap(zeros, T, V ⊗ V ← V ⊗ V) +end + +""" + e_plusmin_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that creates a spin-up electron at the first site and annihilates a spin-up electron at the second. +""" +function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}) + t = two_site_operator(T, Trivial, Trivial) + I = sectortype(t) + t[(I(1), I(0), dual(I(0)), dual(I(1)))][1, 1, 1, 1] = 1 + t[(I(1), I(1), dual(I(0)), dual(I(0)))][1, 2, 1, 2] = 1 + t[(I(0), I(0), dual(I(1)), dual(I(1)))][2, 1, 2, 1] = -1 + t[(I(0), I(1), dual(I(1)), dual(I(0)))][2, 2, 2, 2] = -1 + return t +end +function e_plusmin_up(T, ::Type{Trivial}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{Trivial}, ::Type{SU2Irrep}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{Trivial}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}) + t = two_site_operator(T, U1Irrep, U1Irrep) + I = sectortype(t) + t[(I(1, 1, 1 // 2), I(0, 0, 0), dual(I(0, 0, 0)), dual(I(1, 1, 1 // 2)))] .= 1 + t[(I(1, 1, 1 // 2), I(1, 1, -1 // 2), dual(I(0, 0, 0)), dual(I(0, 2, 0)))] .= 1 + t[(I(0, 2, 0), I(0, 0, 0), dual(I(1, 1, -1 // 2)), dual(I(1, 1, 1 // 2)))] .= -1 + t[(I(0, 2, 0), I(1, 1, -1 // 2), dual(I(1, 1, -1 // 2)), dual(I(0, 2, 0)))] .= -1 + return t +end +function e_plusmin_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{SU2Irrep}, ::Type{Trivial}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_plusmin_up(T, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + return error("Not implemented") +end +const e⁺e⁻ꜛ = e_plusmin_up + +""" + e_plusmin_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that creates a spin-down electron at the first site and annihilates a spin-down electron at the second. +""" +function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial}) + t = two_site_operator(T, Trivial, Trivial) + I = sectortype(t) + t[(I(1), I(0), dual(I(0)), dual(I(1)))][2, 1, 1, 2] = 1 + t[(I(1), I(1), dual(I(0)), dual(I(0)))][2, 1, 1, 2] = -1 + t[(I(0), I(0), dual(I(1)), dual(I(1)))][2, 1, 1, 2] = 1 + t[(I(0), I(1), dual(I(1)), dual(I(0)))][2, 1, 1, 2] = -1 + return t +end +function e_plusmin_down(T, ::Type{Trivial}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{Trivial}, ::Type{SU2Irrep}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{Trivial}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}) + t = two_site_operator(T, U1Irrep, U1Irrep) + I = sectortype(t) + t[(I(1, 1, -1 // 2), I(0, 0, 0), dual(I(0, 0, 0)), dual(I(1, 1, -1 // 2)))] .= 1 + t[(I(1, 1, -1 // 2), I(1, 1, 1 // 2), dual(I(0, 0, 0)), dual(I(0, 2, 0)))] .= -1 + t[(I(0, 2, 0), I(0, 0, 0), dual(I(1, 1, 1 // 2)), dual(I(1, 1, -1 // 2)))] .= 1 + t[(I(0, 2, 0), I(1, 1, 1 // 2), dual(I(1, 1, 1 // 2)), dual(I(0, 2, 0)))] .= -1 + return t +end +function e_plusmin_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{SU2Irrep}, ::Type{Trivial}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_plusmin_down(T, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + return error("Not implemented") +end +const e⁺e⁻ꜜ = e_plusmin_down + +""" + e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that annihilates a spin-up electron at the first site and creates a spin-up electron at the second. +""" +function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + return copy(adjoint(e_plusmin_up(T, particle_symmetry, spin_symmetry))) +end +const e⁻⁺ꜛ = e_minplus_up + +""" + e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that annihilates a spin-down electron at the first site and creates a spin-down electron at the second. +""" +function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + return copy(adjoint(e_plusmin_down(T, particle_symmetry, spin_symmetry))) +end +const e⁻e⁺ꜜ = e_minplus_down + +""" + e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that creates a particle at the first site and annihilates a particle at the second. +This is the sum of `e_plusmin_up` and `e_plusmin_down`. +""" +function e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + return e_plusmin_up(T, particle_symmetry, spin_symmetry) + + e_plusmin_down(T, particle_symmetry, spin_symmetry) +end +function e_plusmin(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + t = two_site_operator(T, U1Irrep, SU2Irrep) + I = sectortype(t) + f1 = only(fusiontrees((I(0, 0, 0), I(1, 1, 1 // 2)), I(1, 1, 1 // 2))) + f2 = only(fusiontrees((I(1, 1, 1 // 2), I(0, 0, 0)), I(1, 1, 1 // 2))) + t[f1, f2] .= 1 + f3 = only(fusiontrees((I(1, 1, 1 // 2), I(0, 2, 0)), I(1, 3, 1 // 2))) + f4 = only(fusiontrees((I(0, 2, 0), I(1, 1, 1 // 2)), I(1, 3, 1 // 2))) + t[f3, f4] .= -1 + f5 = only(fusiontrees((I(0, 0, 0), I(0, 2, 0)), I(0, 2, 0))) + f6 = only(fusiontrees((I(1, 1, 1 // 2), I(1, 1, 1 // 2)), I(0, 2, 0))) + t[f5, f6] .= sqrt(2) + f7 = only(fusiontrees((I(1, 1, 1 // 2), I(1, 1, 1 // 2)), I(0, 2, 0))) + f8 = only(fusiontrees((I(0, 2, 0), I(0, 0, 0)), I(0, 2, 0))) + t[f7, f8] .= sqrt(2) + return t +end +const e⁺e⁻ = e_plusmin + +""" + e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the two-body operator that annihilates a particle at the first site and creates a particle at the second. +This is the sum of `e_minplus_up` and `e_minplus_down`. +""" +function e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + return copy(adjoint(e_plusmin(T, particle_symmetry, spin_symmetry))) +end +const e⁻e⁺ = e_minplus + +""" + e_number_up(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the one-body operator that counts the number of spin-up electrons. +""" +function e_number_up(T, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) + t = single_site_operator(T, Trivial, Trivial) + I = sectortype(t) + t[(I(1), I(1))][1, 1] = 1 + t[(I(0), I(0))][2, 2] = 1 + return t +end +function e_number_up(T, ::Type{Trivial}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_number_up(T, ::Type{Trivial}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) +end +function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}=Trivial) + return error("Not implemented") +end +function e_number_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}) + t = single_site_operator(T, U1Irrep, U1Irrep) + I = sectortype(t) + block(t, I(1, 1, 1 // 2)) .= 1 + block(t, I(0, 2, 0)) .= 1 + return t +end +function e_number_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) +end +function e_number_up(T, ::Type{SU2Irrep}, ::Type{Trivial}=Trivial) + return error("Not implemented") +end +function e_number_up(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_number_up(T, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) +end +const nꜛ = e_number_up + +""" + e_number_down(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the one-body operator that counts the number of spin-down electrons. +""" +function e_number_down(T, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) + t = single_site_operator(T, Trivial, Trivial) + I = sectortype(t) + t[(I(1), I(1))][2, 2] = 1 + t[(I(0), I(0))][2, 2] = 1 + return t +end +function e_number_down(T, ::Type{Trivial}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_number_down(T, ::Type{Trivial}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) +end +function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}=Trivial) + return error("Not implemented") +end +function e_number_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}) + t = single_site_operator(T, U1Irrep, U1Irrep) + I = sectortype(t) + block(t, I(1, 1, -1 // 2)) .= 1 + block(t, I(0, 2, 0)) .= 1 + return t +end +function e_number_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) +end +function e_number_down(T, ::Type{SU2Irrep}, ::Type{Trivial}=Trivial) + return error("Not implemented") +end +function e_number_down(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) + return error("Not implemented") +end +function e_number_down(T, ::Type{SU2Irrep}, ::Type{SU2Irrep}) + throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) +end +const nꜜ = e_number_down + +""" + e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the one-body operator that counts the number of particles. +""" +function e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + return e_number_up(T, particle_symmetry, spin_symmetry) + + e_number_down(T, particle_symmetry, spin_symmetry) +end +function e_number(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + t = single_site_operator(T, U1Irrep, SU2Irrep) + I = sectortype(t) + block(t, I(1, 1, 1 // 2)) .= 1 + block(t, I(0, 2, 0)) .= 2 + return t +end +const n = e_number + +""" + e_number_updown(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) + +Return the one-body operator that counts the number of doubly occupied sites. +""" +function e_number_updown(T, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}) + return e_number_up(T, particle_symmetry, spin_symmetry) * + e_number_down(T, particle_symmetry, spin_symmetry) +end +function e_number_updown(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) + t = single_site_operator(T, U1Irrep, SU2Irrep) + I = sectortype(t) + block(t, I(0, 2, 0)) .= 1 + return t +end +const nꜛnꜜ = e_number_updown + +end \ No newline at end of file diff --git a/src/utility.jl b/src/utility.jl index 9535235..a9801a8 100644 --- a/src/utility.jl +++ b/src/utility.jl @@ -28,3 +28,16 @@ function contract_twosite(L::AbstractTensorMap{<:Any,1,2}, R::AbstractTensorMap{ return H end contract_twosite(L::AbstractTensorMap{<:Any,1,1}, R::AbstractTensorMap{<:Any,1,1}) = L ⊗ R + +""" + split_twosite(O) + +Split a two-site operator into two single-site operators with a connecting auxiliary leg. +""" +function split_twosite(O::AbstractTensorMap{<:Any,2,2}) + U, S, V, = tsvd(O, ((3, 1), (4, 2)); trunc=truncbelow(eps(real(scalartype(O))))) + sqrtS = sqrt(S) + @plansor L[p'; p a] := U[p p'; 1] * sqrtS[1; a] + @plansor R[a p'; p] := sqrtS[a; 1] * V[1; p p'] + return L, R +end From fa6f715d009e3a3279c965e6567f14957445d8b7 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 15:07:56 -0500 Subject: [PATCH 2/6] changes to default values --- src/operators/hubbardoperators.jl | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/src/operators/hubbardoperators.jl b/src/operators/hubbardoperators.jl index f9360e2..9e6d373 100644 --- a/src/operators/hubbardoperators.jl +++ b/src/operators/hubbardoperators.jl @@ -28,7 +28,7 @@ end function hubbard_space(::Type{Trivial}, ::Type{SU2Irrep}) return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) end -function hubbard_space(::Type{U1Irrep}, ::Type{Trivial}=Trivial) +function hubbard_space(::Type{U1Irrep}, ::Type{Trivial}) return Vect[FermionParity ⊠ U1Irrep]((0, 0) => 1, (1, 1) => 2, (0, 2) => 1) end function hubbard_space(::Type{U1Irrep}, ::Type{U1Irrep}) @@ -39,7 +39,7 @@ function hubbard_space(::Type{U1Irrep}, ::Type{SU2Irrep}) return Vect[FermionParity ⊠ U1Irrep ⊠ SU2Irrep]((0, 0, 0) => 1, (1, 1, 1 // 2) => 1, (0, 2, 0) => 1) end -function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial}=Trivial) +function hubbard_space(::Type{SU2Irrep}, ::Type{Trivial}) return Vect[FermionParity ⊠ SU2Irrep]((0, 0) => 2, (1, 1 // 2) => 1) end function hubbard_space(::Type{SU2Irrep}, ::Type{U1Irrep}) @@ -66,6 +66,7 @@ end Return the two-body operator that creates a spin-up electron at the first site and annihilates a spin-up electron at the second. """ +e_plusmin_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_up(ComplexF64, P, S) function e_plusmin_up(T, ::Type{Trivial}, ::Type{Trivial}) t = two_site_operator(T, Trivial, Trivial) I = sectortype(t) @@ -112,6 +113,7 @@ const e⁺e⁻ꜛ = e_plusmin_up Return the two-body operator that creates a spin-down electron at the first site and annihilates a spin-down electron at the second. """ +e_plusmin_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin_down(ComplexF64, P, S) function e_plusmin_down(T, ::Type{Trivial}, ::Type{Trivial}) t = two_site_operator(T, Trivial, Trivial) I = sectortype(t) @@ -158,6 +160,7 @@ const e⁺e⁻ꜜ = e_plusmin_down Return the two-body operator that annihilates a spin-up electron at the first site and creates a spin-up electron at the second. """ +e_minplus_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_up(ComplexF64, P, S) function e_minplus_up(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return copy(adjoint(e_plusmin_up(T, particle_symmetry, spin_symmetry))) end @@ -168,6 +171,7 @@ const e⁻⁺ꜛ = e_minplus_up Return the two-body operator that annihilates a spin-down electron at the first site and creates a spin-down electron at the second. """ +e_minplus_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus_down(ComplexF64, P, S) function e_minplus_down(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return copy(adjoint(e_plusmin_down(T, particle_symmetry, spin_symmetry))) end @@ -179,6 +183,7 @@ const e⁻e⁺ꜜ = e_minplus_down Return the two-body operator that creates a particle at the first site and annihilates a particle at the second. This is the sum of `e_plusmin_up` and `e_plusmin_down`. """ +e_plusmin(P::Type{<:Sector}, S::Type{<:Sector}) = e_plusmin(ComplexF64, P, S) function e_plusmin(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return e_plusmin_up(T, particle_symmetry, spin_symmetry) + e_plusmin_down(T, particle_symmetry, spin_symmetry) @@ -208,6 +213,7 @@ const e⁺e⁻ = e_plusmin Return the two-body operator that annihilates a particle at the first site and creates a particle at the second. This is the sum of `e_minplus_up` and `e_minplus_down`. """ +e_minplus(P::Type{<:Sector}, S::Type{<:Sector}) = e_minplus(ComplexF64, P, S) function e_minplus(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return copy(adjoint(e_plusmin(T, particle_symmetry, spin_symmetry))) end @@ -218,7 +224,8 @@ const e⁻e⁺ = e_minplus Return the one-body operator that counts the number of spin-up electrons. """ -function e_number_up(T, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) +e_number_up(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_up(ComplexF64, P, S) +function e_number_up(T::Type{<:Number}, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) t = single_site_operator(T, Trivial, Trivial) I = sectortype(t) t[(I(1), I(1))][1, 1] = 1 @@ -231,7 +238,7 @@ end function e_number_up(T, ::Type{Trivial}, ::Type{SU2Irrep}) throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}=Trivial) +function e_number_up(T, ::Type{U1Irrep}, ::Type{Trivial}) return error("Not implemented") end function e_number_up(T, ::Type{U1Irrep}, ::Type{U1Irrep}) @@ -244,7 +251,7 @@ end function e_number_up(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) throw(ArgumentError("`e_number_up` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_up(T, ::Type{SU2Irrep}, ::Type{Trivial}=Trivial) +function e_number_up(T, ::Type{SU2Irrep}, ::Type{Trivial}) return error("Not implemented") end function e_number_up(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) @@ -260,7 +267,8 @@ const nꜛ = e_number_up Return the one-body operator that counts the number of spin-down electrons. """ -function e_number_down(T, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) +e_number_down(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_down(ComplexF64, P, S) +function e_number_down(T::Type{<:Number}, ::Type{Trivial}=Trivial, ::Type{Trivial}=Trivial) t = single_site_operator(T, Trivial, Trivial) I = sectortype(t) t[(I(1), I(1))][2, 2] = 1 @@ -273,7 +281,7 @@ end function e_number_down(T, ::Type{Trivial}, ::Type{SU2Irrep}) throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}=Trivial) +function e_number_down(T, ::Type{U1Irrep}, ::Type{Trivial}) return error("Not implemented") end function e_number_down(T, ::Type{U1Irrep}, ::Type{U1Irrep}) @@ -286,7 +294,7 @@ end function e_number_down(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) throw(ArgumentError("`e_number_down` is not symmetric under `SU2Irrep` spin symmetry")) end -function e_number_down(T, ::Type{SU2Irrep}, ::Type{Trivial}=Trivial) +function e_number_down(T, ::Type{SU2Irrep}, ::Type{Trivial}) return error("Not implemented") end function e_number_down(T, ::Type{SU2Irrep}, ::Type{U1Irrep}) @@ -302,6 +310,7 @@ const nꜜ = e_number_down Return the one-body operator that counts the number of particles. """ +e_number(P::Type{<:Sector}, S::Type{<:Sector}) = e_number(ComplexF64, P, S) function e_number(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return e_number_up(T, particle_symmetry, spin_symmetry) + e_number_down(T, particle_symmetry, spin_symmetry) @@ -320,6 +329,7 @@ const n = e_number Return the one-body operator that counts the number of doubly occupied sites. """ +e_number_updown(P::Type{<:Sector}, S::Type{<:Sector}) = e_number_updown(ComplexF64, P, S) function e_number_updown(T, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}) return e_number_up(T, particle_symmetry, spin_symmetry) * From baad7b6d4cb70224331602284ba2f3ec08ec81e3 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 15:08:35 -0500 Subject: [PATCH 3/6] hubbard tests --- test/hubbardoperators.jl | 93 ++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 4 ++ 2 files changed, 97 insertions(+) create mode 100644 test/hubbardoperators.jl diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl new file mode 100644 index 0000000..09edbf6 --- /dev/null +++ b/test/hubbardoperators.jl @@ -0,0 +1,93 @@ +using Test +using TensorKit +using MPSKitModels.HubbardOperators +using LinearAlgebra: eigvals + +implemented_symmetries = [(Trivial, Trivial), + (U1Irrep, U1Irrep), + (U1Irrep, SU2Irrep)] +@testset "basic properties" begin + for particle_symmetry in (Trivial, U1Irrep, SU2Irrep), + spin_symmetry in (Trivial, U1Irrep, SU2Irrep) + + if (particle_symmetry, spin_symmetry) in implemented_symmetries + # test hermiticity + @test e_plusmin(particle_symmetry, spin_symmetry)' ≈ + e_minplus(particle_symmetry, spin_symmetry) + if spin_symmetry !== SU2Irrep + @test e_plusmin_down(particle_symmetry, spin_symmetry)' ≈ + e_minplus_down(particle_symmetry, spin_symmetry) + @test e_plusmin_up(particle_symmetry, spin_symmetry)' ≈ + e_minplus_up(particle_symmetry, spin_symmetry) + @test e_plusmin_down(particle_symmetry, spin_symmetry)' ≈ + e_minplus_down(particle_symmetry, spin_symmetry) + @test e_plusmin_up(particle_symmetry, spin_symmetry)' ≈ + e_minplus_up(particle_symmetry, spin_symmetry) + end + + # test number operator + if spin_symmetry !== SU2Irrep + @test e_number(particle_symmetry, spin_symmetry) ≈ + e_number_up(particle_symmetry, spin_symmetry) + + e_number_down(particle_symmetry, spin_symmetry) + @test e_number_updown(particle_symmetry, spin_symmetry) ≈ + e_number_up(particle_symmetry, spin_symmetry) * + e_number_down(particle_symmetry, spin_symmetry) ≈ + e_number_down(particle_symmetry, spin_symmetry) * + e_number_up(particle_symmetry, spin_symmetry) + end + else + @test_broken e_plusmin(particle_symmetry, spin_symmetry) + @test_broken e_minplus(particle_symmetry, spin_symmetry) + end + end +end + +function hamiltonian(particle_symmetry, spin_symmetry; t, U, mu, L) + hopping = t * (e_plusmin(particle_symmetry, spin_symmetry) + + e_minplus(particle_symmetry, spin_symmetry)) + interaction = U * e_number_updown(particle_symmetry, spin_symmetry) + chemical_potential = mu * e_number(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)) + end + sum(1:L) do i + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, interaction)) + end + + sum(1:L) do i + return reduce(⊗, insert!(collect(Any, fill(I, L - 1)), i, chemical_potential)) + end + return H +end + +@testset "spectrum" begin + L = 4 + t = randn() + U = randn() + mu = randn() + + hopping = t * (e_plusmin(Trivial, Trivial) + e_minplus(Trivial, Trivial)) + interaction = U * e_number_updown(Trivial, Trivial) + chemical_potential = mu * e_number(Trivial, Trivial) + I = id(domain(interaction)) + + H_triv = hamiltonian(Trivial, Trivial; t, U, mu, L) + vals_triv = mapreduce(vcat, eigvals(H_triv)) do (c, v) + return repeat(real.(v), dim(c)) + end + sort!(vals_triv) + + H_u1_u1 = hamiltonian(U1Irrep, U1Irrep; t, U, mu, L) + vals_u1_u1 = mapreduce(vcat, eigvals(H_u1_u1)) do (c, v) + return repeat(real.(v), dim(c)) + end + sort!(vals_u1_u1) + @test vals_triv ≈ vals_u1_u1 + + H_u1_su2 = hamiltonian(U1Irrep, SU2Irrep; t, U, mu, L) + vals_u1_su2 = mapreduce(vcat, eigvals(H_u1_su2)) do (c, v) + return repeat(real.(v), dim(c)) + end + sort!(vals_u1_su2) + @test vals_triv ≈ vals_u1_su2 +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index fddb23c..c897426 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,10 @@ end include("fermionoperators.jl") end +@testset "Hubbard operators" begin + include("hubbardoperators.jl") +end + @testset "mpoham" begin include("mpoham.jl") end From 3f8fc15105a4ab6f7e6e45a0ce8f05f3b9a608cd Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 17:42:58 -0500 Subject: [PATCH 4/6] Remove obsolete Hubbard tests --- test/fermionoperators.jl | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/test/fermionoperators.jl b/test/fermionoperators.jl index f3ef3ed..212c48e 100644 --- a/test/fermionoperators.jl +++ b/test/fermionoperators.jl @@ -29,18 +29,3 @@ using MPSKitModels: contract_twosite, contract_onesite @test c_number() ≈ contract_onesite(c⁺(; side=:L), c⁻(; side=:R)) end - -@testset "electrons $particle_symmetry, $spin_symmetry" for (particle_symmetry, - spin_symmetry) in - ((Trivial, Trivial), - (U1Irrep, SU2Irrep)) - ee⁺ = e_minplus(particle_symmetry, spin_symmetry) - e⁺e = e_plusmin(particle_symmetry, spin_symmetry) - @test ee⁺' ≈ e⁺e - @test (e⁺e + ee⁺)' ≈ ee⁺ + e⁺e - @test (e⁺e - ee⁺)' ≈ ee⁺ - e⁺e - - @test e_number(particle_symmetry, spin_symmetry) ≈ - contract_onesite(e⁺(particle_symmetry, spin_symmetry; side=:L), - e⁻(particle_symmetry, spin_symmetry; side=:R)) -end From 93e3a464a80125174733286674e65ba9b31dc562 Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 17:43:58 -0500 Subject: [PATCH 5/6] Remove unused exports --- src/MPSKitModels.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MPSKitModels.jl b/src/MPSKitModels.jl index 7582214..bf61ede 100644 --- a/src/MPSKitModels.jl +++ b/src/MPSKitModels.jl @@ -34,7 +34,7 @@ export c_plus, c_min, c_plusplus, c_minmin, c_plusmin, c_minplus, c_number export c⁺, c⁻, c⁺⁺, c⁻⁻, c⁺⁻, c⁻⁺ export e_plus, e_min, e_plusplus, e_minmin, e_plusmin, e_minplus export e_number, e_number_up, e_number_down, e_number_updown -export e⁺, e⁻, e⁺⁺, e⁻⁻, e⁺⁻, e⁻⁺ +export e⁺⁺, e⁻⁻, e⁺⁻, e⁻⁺ export transverse_field_ising export kitaev_model From 67ad3dcc875b6280078526f4907db0d04266d75f Mon Sep 17 00:00:00 2001 From: Lukas Devos Date: Thu, 14 Nov 2024 17:49:23 -0500 Subject: [PATCH 6/6] Formatter --- src/operators/hubbardoperators.jl | 2 +- test/hubbardoperators.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/operators/hubbardoperators.jl b/src/operators/hubbardoperators.jl index 9e6d373..6b60c9f 100644 --- a/src/operators/hubbardoperators.jl +++ b/src/operators/hubbardoperators.jl @@ -343,4 +343,4 @@ function e_number_updown(T, ::Type{U1Irrep}, ::Type{SU2Irrep}) end const nꜛnꜜ = e_number_updown -end \ No newline at end of file +end diff --git a/test/hubbardoperators.jl b/test/hubbardoperators.jl index 09edbf6..105ac22 100644 --- a/test/hubbardoperators.jl +++ b/test/hubbardoperators.jl @@ -90,4 +90,4 @@ end end sort!(vals_u1_su2) @test vals_triv ≈ vals_u1_su2 -end \ No newline at end of file +end