From 325a30b74a1810d390035b2b07e2173f0d500fab Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 30 Jan 2026 06:31:36 -0500 Subject: [PATCH] Very basic CUDA support --- Project.toml | 13 ++- ext/MPSKitModelsCUDAExt.jl | 183 ++++++++++++++++++++++++++++++ src/models/hamiltonians.jl | 68 +++++------ src/operators/bosonoperators.jl | 130 ++++++++++++--------- src/operators/fermionoperators.jl | 30 ++--- src/operators/localoperators.jl | 2 + src/operators/spinoperators.jl | 178 ++++++++++++++--------------- 7 files changed, 410 insertions(+), 194 deletions(-) create mode 100644 ext/MPSKitModelsCUDAExt.jl diff --git a/Project.toml b/Project.toml index 7cf0447..c7ac825 100644 --- a/Project.toml +++ b/Project.toml @@ -12,7 +12,16 @@ TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" +[weakdeps] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +cuTENSOR = "011b41b2-24ef-40a8-b3eb-fa098493e9e1" + +[extensions] +MPSKitModelsCUDAExt = ["CUDA", "cuTENSOR"] + [compat] +CUDA = "5" +cuTENSOR = "2" MPSKit = "0.13.1" MacroTools = "0.5" PrecompileTools = "1" @@ -23,9 +32,11 @@ TupleTools = "1" julia = "1.10" [extras] +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +cuTENSOR = "011b41b2-24ef-40a8-b3eb-fa098493e9e1" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TestExtras = "5ed8adda-3752-4e41-b88a-e8b09835ee3a" [targets] -test = ["Test", "TestExtras", "QuadGK"] +test = ["Test", "TestExtras", "QuadGK", "CUDA", "cuTENSOR"] diff --git a/ext/MPSKitModelsCUDAExt.jl b/ext/MPSKitModelsCUDAExt.jl new file mode 100644 index 0000000..5b8919b --- /dev/null +++ b/ext/MPSKitModelsCUDAExt.jl @@ -0,0 +1,183 @@ +module MPSKitModelsCUDAExt + +using MPSKitModels, CUDA, cuTENSOR + +import MPSKitModels: build_a_plus_left!, build_a_plus_right!, build_a_min_left!, build_a_min_right! +import MPSKitModels: e_plusmin_up, e_plusmin_down, e_number_up, e_number_down, spinmatrices +using MPSKitModels: SU2Irrep, U1Irrep, Trivial, fusiontrees, sectortype, block, dual + +function e_number_down(::Type{TA}, ::Type{Trivial} = Trivial, ::Type{Trivial} = Trivial) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, Trivial, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1), I(1))][2, 2] = 1 + t[(I(0), I(0))][2, 2] = 1 + end + return t +end +function e_number_down(::Type{TA}, ::Type{Trivial}, ::Type{U1Irrep}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, Trivial, U1Irrep) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1, -1 // 2), dual(I(1, -1 // 2)))][1, 1] = 1 + t[(I(0, 0), I(0, 0))][2, 2] = 1 + end + return t +end +function e_number_down(::Type{TA}, ::Type{U1Irrep}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, U1Irrep, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + block(t, I(1, 1))[2, 2] = 1 # expected to be [1,2] + block(t, I(0, 2))[1, 1] = 1 + end + return t +end +function e_number_down(::Type{TA}, ::Type{U1Irrep}, ::Type{U1Irrep}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, U1Irrep, U1Irrep) + I = sectortype(t) + CUDA.@allowscalar begin + block(t, I(1, 1, -1 // 2)) .= 1 + block(t, I(0, 2, 0)) .= 1 + end + return t +end +function e_number_up(::Type{TA}, ::Type{Trivial}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, Trivial, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1), I(1))][1, 1] = 1 + t[(I(0), I(0))][2, 2] = 1 + end + return t +end +function e_number_up(::Type{TA}, ::Type{U1Irrep}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, U1Irrep, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + block(t, I(1, 1))[1, 1] = 1 + block(t, I(0, 2))[1, 1] = 1 + end + return t +end +function e_number_up(::Type{TA}, ::Type{U1Irrep}, ::Type{U1Irrep}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.single_site_operator(TA, U1Irrep, U1Irrep) + I = sectortype(t) + CUDA.@allowscalar begin + block(t, I(1, 1, 1 // 2)) .= 1 + block(t, I(0, 2, 0)) .= 1 + end + return t +end +function e_plusmin_up(::Type{TA}, ::Type{Trivial}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.two_site_operator(TA, Trivial, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + 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 + end + return t +end +function e_plusmin_up(::Type{TA}, ::Type{U1Irrep}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.two_site_operator(TA, U1Irrep, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1, 1), I(0, 0), dual(I(0, 0)), dual(I(1, 1)))][1, 1, 1, 1] = 1 + t[(I(1, 1), I(1, 1), dual(I(0, 0)), dual(I(0, 2)))][1, 2, 1, 1] = 1 + t[(I(0, 2), I(0, 0), dual(I(1, 1)), dual(I(1, 1)))][1, 1, 2, 1] = -1 + t[(I(0, 2), I(1, 1), dual(I(1, 1)), dual(I(0, 2)))][1, 2, 2, 1] = -1 + end + return t +end + +function e_plusmin_down(::Type{TA}, ::Type{Trivial}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.two_site_operator(TA, Trivial, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + 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 + end + return t +end +function e_plusmin_down(::Type{TA}, ::Type{Trivial}, ::Type{U1Irrep}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.two_site_operator(T, Trivial, U1Irrep) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1, -1 // 2), I(0, 0), dual(I(0, 0)), dual(I(1, -1 // 2)))][1, 1, 1, 1] = 1 + t[(I(1, -1 // 2), I(1, 1 // 2), dual(I(0, 0)), dual(I(0, 0)))][1, 1, 1, 2] = -1 + t[(I(0, 0), I(0, 0), dual(I(1, 1 // 2)), dual(I(1, -1 // 2)))][2, 1, 1, 1] = 1 + t[(I(0, 0), I(1, 1 // 2), dual(I(1, 1 // 2)), dual(I(0, 0)))][2, 1, 1, 2] = -1 + end + return t +end +function e_plusmin_down(::Type{TA}, ::Type{U1Irrep}, ::Type{Trivial}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.two_site_operator(TA, U1Irrep, Trivial) + I = sectortype(t) + CUDA.@allowscalar begin + t[(I(1, 1), I(0, 0), dual(I(0, 0)), dual(I(1, 1)))][2, 1, 1, 2] = 1 + t[(I(1, 1), I(1, 1), dual(I(0, 0)), dual(I(0, 2)))][2, 1, 1, 1] = -1 + t[(I(0, 2), I(0, 0), dual(I(1, 1)), dual(I(1, 1)))][1, 1, 1, 2] = 1 + t[(I(0, 2), I(1, 1), dual(I(1, 1)), dual(I(0, 2)))][1, 1, 1, 1] = -1 + end + return t +end +function e_plusmin_down(::Type{TA}, ::Type{U1Irrep}, ::Type{U1Irrep}) where {TA <: CuArray} + t = MPSKitModels.HubbardOperators.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 build_a_plus_left!(::Type{U1Irrep}, ::Type{TA}, a⁺) where {TA <: CuArray} + for (f1, f2) in fusiontrees(a⁺) + c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] + if c₁.charge == c₂.charge + 1 + CUDA.@allowscalar a⁺[f1, f2] .= -sqrt(c₁.charge) + end + end + return +end + +function build_a_plus_right!(::Type{U1Irrep}, ::Type{TA}, a⁺) where {TA <: CuArray} + for (f1, f2) in fusiontrees(a⁺) + c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] + if c₁.charge == c₂.charge + 1 + CUDA.@allowscalar a⁺[f1, f2] .= -sqrt(c₁.charge) + end + end + return +end + +function build_a_min_left!(::Type{U1Irrep}, ::Type{TA}, a⁻) where {TA <: CuArray} + for (f1, f2) in fusiontrees(a⁻) + c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] + if c₁.charge + 1 == c₂.charge + CUDA.@allowscalar a⁻[f1, f2] .= -sqrt(c₂.charge) + end + end + return +end + +function build_a_min_right!(::Type{U1Irrep}, ::Type{TA}, a⁻) where {TA <: CuArray} + for (f1, f2) in fusiontrees(a⁻) + c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] + if c₁.charge + 1 == c₂.charge + CUDA.@allowscalar a⁻[f1, f2] .= -sqrt(c₂.charge) + end + end + return +end + +function spinmatrices(s::Union{Rational{Int}, Int}, ::Type{TorA}) where {T, TorA <: CuArray{T}} + hSx, hSy, hSz, hI = spinmatrices(s, T) + return CuArray(hSx), CuArray(hSy), CuArray(hSz), CuArray(hI) +end + +end diff --git a/src/models/hamiltonians.jl b/src/models/hamiltonians.jl index 7d36c28..fe6b04a 100644 --- a/src/models/hamiltonians.jl +++ b/src/models/hamiltonians.jl @@ -205,12 +205,12 @@ function heisenberg_XYZ(lattice::AbstractLattice; kwargs...) return heisenberg_XYZ(ComplexF64, lattice; kwargs...) end function heisenberg_XYZ( - T::Type{<:Number} = ComplexF64, lattice::AbstractLattice = InfiniteChain(1); + ::Type{TorA}, lattice::AbstractLattice = InfiniteChain(1); Jx = 1.0, Jy = 1.0, Jz = 1.0, spin = 1 - ) - term = rmul!(S_xx(T, Trivial; spin = spin), Jx) + - rmul!(S_yy(T, Trivial; spin = spin), Jy) + - rmul!(S_zz(T, Trivial; spin = spin), Jz) + ) where {TorA} + term = rmul!(S_xx(TorA, Trivial; spin = spin), Jx) + + rmul!(S_yy(TorA, Trivial; spin = spin), Jy) + + rmul!(S_zz(TorA, Trivial; spin = spin), Jz) return @mpoham sum(nearest_neighbours(lattice)) do (i, j) return term{i, j} end @@ -239,19 +239,19 @@ function bilinear_biquadratic_model( return bilinear_biquadratic_model(ComplexF64, symmetry, lattice; kwargs...) end function bilinear_biquadratic_model( - elt::Type{<:Number}, lattice::AbstractLattice; + ::Type{TorA}, lattice::AbstractLattice; kwargs... - ) - return bilinear_biquadratic_model(elt, Trivial, lattice; kwargs...) + ) where {TorA} + return bilinear_biquadratic_model(TorA, Trivial, lattice; kwargs...) end function bilinear_biquadratic_model( - elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + ::Type{TorA}, symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); spin = 1, J = 1.0, θ = 0.0 - ) + ) where {TorA} return @mpoham sum(nearest_neighbours(lattice)) do (i, j) - return J * cos(θ) * S_exchange(elt, symmetry; spin = spin){i, j} + - J * sin(θ) * (S_exchange(elt, symmetry; spin = spin)^2){i, j} + return J * cos(θ) * S_exchange(TorA, symmetry; spin = spin){i, j} + + J * sin(θ) * (S_exchange(TorA, symmetry; spin = spin)^2){i, j} end end @@ -328,18 +328,18 @@ function hubbard_model( ) return hubbard_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) end -function hubbard_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) - return hubbard_model(elt, Trivial, Trivial, lattice; kwargs...) +function hubbard_model(::Type{TorA}, lattice::AbstractLattice; kwargs...) where {TorA} + return hubbard_model(TorA, Trivial, Trivial, lattice; kwargs...) end function hubbard_model( - T::Type{<:Number} = ComplexF64, particle_symmetry::Type{<:Sector} = Trivial, + ::Type{TorA}, particle_symmetry::Type{<:Sector} = Trivial, spin_symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); t = 1.0, U = 1.0, mu = 0.0, n::Integer = 0 - ) - hopping = e⁺e⁻(T, particle_symmetry, spin_symmetry) + - e⁻e⁺(T, particle_symmetry, spin_symmetry) - interaction_term = nꜛnꜜ(T, particle_symmetry, spin_symmetry) - N = e_number(T, particle_symmetry, spin_symmetry) + ) where {TorA} + hopping = e⁺e⁻(TorA, particle_symmetry, spin_symmetry) + + e⁻e⁺(TorA, particle_symmetry, spin_symmetry) + interaction_term = nꜛnꜜ(TorA, particle_symmetry, spin_symmetry) + N = e_number(TorA, particle_symmetry, spin_symmetry) return @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) return -t * hopping{i, j} @@ -377,14 +377,14 @@ function bose_hubbard_model( return bose_hubbard_model(ComplexF64, symmetry, lattice; kwargs...) end function bose_hubbard_model( - elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + ::Type{TorA}, symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, n = 0 - ) - hopping_term = a_plusmin(elt, symmetry; cutoff = cutoff) + - a_minplus(elt, symmetry; cutoff = cutoff) - N = a_number(elt, symmetry; cutoff = cutoff) - interaction_term = contract_onesite(N, N - id(domain(N))) + ) where {TorA} + hopping_term = a_plusmin(TorA, symmetry; cutoff = cutoff) + + a_minplus(TorA, symmetry; cutoff = cutoff) + N = a_number(TorA, symmetry; cutoff = cutoff) + interaction_term = contract_onesite(N, N - id(TorA, domain(N))) H = @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) @@ -436,19 +436,19 @@ function tj_model( ) return tj_model(ComplexF64, particle_symmetry, spin_symmetry; kwargs...) end -function tj_model(elt::Type{<:Number}, lattice::AbstractLattice; kwargs...) - return tj_model(elt, Trivial, Trivial, lattice; kwargs...) +function tj_model(::Type{TorA}, lattice::AbstractLattice; kwargs...) where {TorA} + return tj_model(TorA, Trivial, Trivial, lattice; kwargs...) end function tj_model( - T::Type{<:Number} = ComplexF64, particle_symmetry::Type{<:Sector} = Trivial, + ::Type{TorA}, particle_symmetry::Type{<:Sector} = Trivial, spin_symmetry::Type{<:Sector} = Trivial, lattice::AbstractLattice = InfiniteChain(1); t = 2.5, J = 1.0, mu = 0.0, slave_fermion::Bool = false - ) - hopping = TJOperators.e_plusmin(T, particle_symmetry, spin_symmetry; slave_fermion) + - TJOperators.e_minplus(T, particle_symmetry, spin_symmetry; slave_fermion) - num = TJOperators.e_number(T, particle_symmetry, spin_symmetry; slave_fermion) + ) where {TorA} + hopping = TJOperators.e_plusmin(TorA, particle_symmetry, spin_symmetry; slave_fermion) + + TJOperators.e_minplus(TorA, particle_symmetry, spin_symmetry; slave_fermion) + num = TJOperators.e_number(TorA, particle_symmetry, spin_symmetry; slave_fermion) heisenberg = TJOperators.S_exchange( - T, particle_symmetry, spin_symmetry; + TorA, particle_symmetry, spin_symmetry; slave_fermion ) - (1 / 4) * (num ⊗ num) diff --git a/src/operators/bosonoperators.jl b/src/operators/bosonoperators.jl index 0f6c4ac..507545f 100644 --- a/src/operators/bosonoperators.jl +++ b/src/operators/bosonoperators.jl @@ -5,38 +5,48 @@ The truncated bosonic creation operator, with a maximum of `cutoff` bosons per site. """ function a_plus end -a_plus(; kwargs...) = a_plus(ComplexF64, Trivial; kwargs...) -a_plus(elt::Type{<:Number}; kwargs...) = a_plus(elt, Trivial; kwargs...) -a_plus(symm::Type{<:Sector}; kwargs...) = a_plus(ComplexF64, symm; kwargs...) +a_plus(; kwargs...) = a_plus(ComplexF64, Vector{ComplexF64}, Trivial; kwargs...) +a_plus(::Type{T}, ::Type{TA}; kwargs...) where {T <: Number, TA <: AbstractArray{T}} = a_plus(T, TA, Trivial; kwargs...) +a_plus(symm::Type{<:Sector}; kwargs...) = a_plus(ComplexF64, Vector{ComplexF64}, symm; kwargs...) -function a_plus(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) - a⁺ = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) +function a_plus(::Type{TorA}, ::Type{Trivial}; cutoff::Integer = 5) where {TorA} + a⁺ = zeros(TorA, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 1:cutoff a⁺[n + 1, n] = sqrt(n) end return a⁺ end -function a_plus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) +function build_a_plus_left!(::Type{U1Irrep}, ::Type{TA}, a⁺) where {TA} + for (f1, f2) in fusiontrees(a⁺) + c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] + if c₁.charge == c₂.charge + 1 + a⁺[f1, f2] .= -sqrt(c₁.charge) + end + end + return +end + +function build_a_plus_right!(::Type{U1Irrep}, ::Type{TA}, a⁺) where {TA} + for (f1, f2) in fusiontrees(a⁺) + c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] + if c₁.charge == c₂.charge + 1 + a⁺[f1, f2] .= -sqrt(c₁.charge) + end + end + return +end + +function a_plus(::Type{TorA}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) where {TorA} pspace = U1Space(n => 1 for n in 0:cutoff) if side === :L vspace = U1Space(1 => 1) - a⁺ = zeros(elt, pspace ← pspace ⊗ vspace) - for (f1, f2) in fusiontrees(a⁺) - c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] - if c₁.charge == c₂.charge + 1 - a⁺[f1, f2] .= -sqrt(c₁.charge) - end - end + a⁺ = zeros(TorA, pspace ← pspace ⊗ vspace) + build_a_plus_left!(U1Irrep, TorA, a⁺) elseif side === :R vspace = U1Space(-1 => 1) - a⁺ = zeros(elt, vspace ⊗ pspace ← pspace) - for (f1, f2) in fusiontrees(a⁺) - c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] - if c₁.charge == c₂.charge + 1 - a⁺[f1, f2] .= -sqrt(c₁.charge) - end - end + a⁺ = zeros(TorA, vspace ⊗ pspace ← pspace) + build_a_plus_right!(U1Irrep, TorA, a⁺) else throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) end @@ -53,37 +63,47 @@ The truncated bosonic annihilation operator, with a maximum of `cutoff` bosons p """ function a_min end a_min(; kwargs...) = a_min(ComplexF64, Trivial; kwargs...) -a_min(elt::Type{<:Number}; kwargs...) = a_min(elt, Trivial; kwargs...) +a_min(::Type{TorA}; kwargs...) where {TorA} = a_min(TorA, Trivial; kwargs...) a_min(symm::Type{<:Sector}; kwargs...) = a_min(ComplexF64, symm; kwargs...) -function a_min(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) - a⁻ = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) +function a_min(::Type{TorA}, ::Type{Trivial}; cutoff::Integer = 5) where {TorA} + a⁻ = zeros(TorA, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 1:cutoff a⁻[n, n + 1] = sqrt(n) end return a⁻ end -function a_min(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) +function build_a_min_left!(::Type{U1Irrep}, ::Type{TA}, a⁻) where {TA} + for (f1, f2) in fusiontrees(a⁻) + c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] + if c₁.charge + 1 == c₂.charge + a⁻[f1, f2] .= -sqrt(c₂.charge) + end + end + return +end + +function build_a_min_right!(::Type{U1Irrep}, ::Type{TA}, a⁻) where {TA} + for (f1, f2) in fusiontrees(a⁻) + c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] + if c₁.charge + 1 == c₂.charge + a⁻[f1, f2] .= -sqrt(c₂.charge) + end + end + return +end + +function a_min(::Type{TorA}, ::Type{U1Irrep}; cutoff::Integer = 5, side = :L) where {TorA} pspace = U1Space(n => 1 for n in 0:cutoff) if side === :L vspace = U1Space(-1 => 1) - a⁻ = zeros(elt, pspace ← pspace ⊗ vspace) - for (f1, f2) in fusiontrees(a⁻) - c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] - if c₁.charge + 1 == c₂.charge - a⁻[f1, f2] .= -sqrt(c₂.charge) - end - end + a⁻ = zeros(TorA, pspace ← pspace ⊗ vspace) + build_a_min_left!(U1Irrep, TorA, a⁻) elseif side === :R vspace = U1Space(1 => 1) - a⁻ = zeros(elt, vspace ⊗ pspace ← pspace) - for (f1, f2) in fusiontrees(a⁻) - c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] - if c₁.charge + 1 == c₂.charge - a⁻[f1, f2] .= -sqrt(c₂.charge) - end - end + a⁻ = zeros(TorA, vspace ⊗ pspace ← pspace) + build_a_min_right!(U1Irrep, TorA, a⁻) else throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) end @@ -94,29 +114,29 @@ const a⁻ = a_min a_plusmin(symmetry::Type{<:Sector}; kwargs...) = a_plusmin(ComplexF64, symmetry; kwargs...) function a_plusmin( - elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + ::Type{TorA}, ::Type{Trivial} = Trivial; cutoff::Integer = 5 - ) - return contract_twosite(a⁺(elt; cutoff = cutoff), a⁻(elt; cutoff = cutoff)) + ) where {TorA} + return contract_twosite(a⁺(TorA; cutoff = cutoff), a⁻(elt; cutoff = cutoff)) end -function a_plusmin(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) +function a_plusmin(::Type{TorA}, ::Type{U1Irrep}; cutoff::Integer = 5) where {TorA} return contract_twosite( - a⁺(elt, U1Irrep; cutoff = cutoff, side = :L), - a⁻(elt, U1Irrep; cutoff = cutoff, side = :R) + a⁺(TorA, U1Irrep; cutoff = cutoff, side = :L), + a⁻(TorA, U1Irrep; cutoff = cutoff, side = :R) ) end a_minplus(symmetry::Type{<:Sector}; kwargs...) = a_minplus(ComplexF64, symmetry; kwargs...) function a_minplus( - elt::Type{<:Number} = ComplexF64, ::Type{Trivial} = Trivial; + ::Type{TorA}, ::Type{Trivial} = Trivial; cutoff::Integer = 5 - ) - return contract_twosite(a⁻(elt; cutoff = cutoff), a⁺(elt; cutoff = cutoff)) + ) where {TorA} + return contract_twosite(a⁻(TorA; cutoff = cutoff), a⁺(elt; cutoff = cutoff)) end -function a_minplus(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) +function a_minplus(::Type{TorA}, ::Type{U1Irrep}; cutoff::Integer = 5) where {TorA} return contract_twosite( - a⁻(elt, U1Irrep; cutoff = cutoff, side = :L), - a⁺(elt, U1Irrep; cutoff = cutoff, side = :R) + a⁻(TorA, U1Irrep; cutoff = cutoff, side = :L), + a⁺(TorA, U1Irrep; cutoff = cutoff, side = :R) ) end @@ -127,20 +147,20 @@ The truncated bosonic number operator, with a maximum of `cutoff` bosons per sit """ function a_number end a_number(; kwargs...) = a_number(ComplexF64, Trivial; kwargs...) -a_number(elt::Type{<:Number}; kwargs...) = a_number(elt, Trivial; kwargs...) +a_number(::Type{TorA}; kwargs...) where {TorA} = a_number(TorA, Trivial; kwargs...) a_number(symm::Type{<:Sector}; kwargs...) = a_number(ComplexF64, symm; kwargs...) -function a_number(elt::Type{<:Number}, ::Type{Trivial}; cutoff::Integer = 5) - N = zeros(elt, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) +function a_number(::Type{TorA}, ::Type{Trivial}; cutoff::Integer = 5) where {TorA} + N = zeros(TorA, ComplexSpace(cutoff + 1), ComplexSpace(cutoff + 1)) for n in 0:cutoff N[n + 1, n + 1] = n end return N end -function a_number(elt::Type{<:Number}, ::Type{U1Irrep}; cutoff::Integer = 5) +function a_number(::Type{TorA}, ::Type{U1Irrep}; cutoff::Integer = 5) where {TorA} pspace = U1Space(n => 1 for n in 0:cutoff) - N = zeros(elt, pspace, pspace) + N = zeros(TorA, pspace, pspace) for (c, b) in blocks(N) b .= c.charge end diff --git a/src/operators/fermionoperators.jl b/src/operators/fermionoperators.jl index 8f148f9..6ff364c 100644 --- a/src/operators/fermionoperators.jl +++ b/src/operators/fermionoperators.jl @@ -8,15 +8,15 @@ Fermionic creation operator. """ -function c_plus(elt::Type{<:Number} = ComplexF64; side = :L) +function c_plus(::Type{TorA}; side = :L) where {TorA} vspace = Vect[fℤ₂](1 => 1) if side === :L pspace = Vect[fℤ₂](0 => 1, 1 => 1) - c⁺ = zeros(elt, pspace ← pspace ⊗ vspace) + c⁺ = zeros(TorA, pspace ← pspace ⊗ vspace) block(c⁺, fℤ₂(1)) .= one(elt) elseif side === :R - C = c_plus(elt; side = :L) - F = isomorphism(storagetype(C), vspace, flip(vspace)) + C = c_plus(TorA; side = :L) + F = isomorphism(TorA, vspace, flip(vspace)) @planar c⁺[-1 -2; -3] := C[-2; 1 2] * τ[1 2; 3 -3] * F[3; -1] else throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) @@ -31,13 +31,13 @@ const c⁺ = c_plus Fermionic annihilation operator. """ -function c_min(elt::Type{<:Number} = ComplexF64; side = :L) +function c_min(::Type{TorA}; side = :L) where {TorA} if side === :L - C = c_plus(elt; side = :L)' - F = isomorphism(flip(space(C, 2)), space(C, 2)) + C = c_plus(TorA; side = :L)' + F = isomorphism(TorA, flip(space(C, 2)), space(C, 2)) @planar c⁻[-1; -2 -3] := C[-1 1; -2] * F[-3; 1] elseif side === :R - c⁻ = permute(c_plus(elt; side = :L)', ((2, 1), (3,))) + c⁻ = permute(c_plus(TorA; side = :L)', ((2, 1), (3,))) else throw(ArgumentError("invalid side `:$side`, expected `:L` or `:R`")) end @@ -46,13 +46,13 @@ end const c⁻ = c_min -c_plusmin(elt = ComplexF64) = contract_twosite(c⁺(elt; side = :L), c⁻(elt; side = :R)) +c_plusmin(::Type{TorA}) where {TorA} = contract_twosite(c⁺(TorA; side = :L), c⁻(TorA; side = :R)) const c⁺c⁻ = c_plusmin -c_minplus(elt = ComplexF64) = contract_twosite(c⁻(elt; side = :L), c⁺(elt; side = :R)) +c_minplus(::Type{TorA}) where {TorA} = contract_twosite(c⁻(TorA; side = :L), c⁺(TorA; side = :R)) const c⁻c⁺ = c_minplus -c_plusplus(elt = ComplexF64) = contract_twosite(c⁺(elt; side = :L), c⁺(elt; side = :R)) +c_plusplus(::Type{TorA}) where {TorA} = contract_twosite(c⁺(TorA; side = :L), c⁺(TorA; side = :R)) const c⁺c⁺ = c_plusplus -c_minmin(elt = ComplexF64) = contract_twosite(c⁻(elt; side = :L), c⁻(elt; side = :R)) +c_minmin(::Type{TorA}) where {TorA} = contract_twosite(c⁻(TorA; side = :L), c⁻(TorA; side = :R)) const c⁻c⁻ = c_minmin """ @@ -60,9 +60,9 @@ const c⁻c⁻ = c_minmin Fermionic number operator. """ -function c_number(elt::Type{<:Number} = ComplexF64) +function c_number(::Type{TorA}) where {TorA} pspace = Vect[fℤ₂](0 => 1, 1 => 1) - n = zeros(elt, pspace ← pspace) - block(n, fℤ₂(1)) .= one(elt) + n = zeros(TorA, pspace ← pspace) + block(n, fℤ₂(1)) .= one(eltype(TorA)) return n end diff --git a/src/operators/localoperators.jl b/src/operators/localoperators.jl index a52e701..712813c 100644 --- a/src/operators/localoperators.jl +++ b/src/operators/localoperators.jl @@ -23,6 +23,8 @@ struct LocalOperator{T <: AbstractTensorMap{<:Number, <:Any, 2, 2}, G <: Lattice return new{T, G}(O, inds) end end +TensorKit.storagetype(lo::LocalOperator{T, G}) where {T, G} = storagetype(T) +TensorKit.storagetype(::Type{LocalOperator{T, G}}) where {T, G} = storagetype(T) function LocalOperator( t::AbstractTensorMap{<:Number, <:Any, N, N}, inds::Vararg{G, N} diff --git a/src/operators/spinoperators.jl b/src/operators/spinoperators.jl index e2e1fe5..2523d9e 100644 --- a/src/operators/spinoperators.jl +++ b/src/operators/spinoperators.jl @@ -17,12 +17,12 @@ end the spinmatrices according to [Wikipedia](https://en.wikipedia.org/wiki/Spin_(physics)#Higher_spins). """ -function spinmatrices(s::Union{Rational{Int}, Int}, elt = ComplexF64) +function spinmatrices(s::Union{Rational{Int}, Int}, ::Type{TorA}) where {TorA <: Number} N = Int(2 * s) - Sx = zeros(elt, N + 1, N + 1) - Sy = zeros(complex(elt), N + 1, N + 1) - Sz = zeros(elt, N + 1, N + 1) + Sx = zeros(TorA, N + 1, N + 1) + Sy = zeros(complex(TorA), N + 1, N + 1) + Sz = zeros(TorA, N + 1, N + 1) for row in 1:(N + 1) for col in 1:(N + 1) @@ -56,29 +56,29 @@ See also [`σˣ`](@ref). """ function S_x end S_x(; kwargs...) = S_x(ComplexF64, Trivial; kwargs...) -S_x(elt::Type{<:Number}; kwargs...) = S_x(elt, Trivial; kwargs...) +S_x(::Type{TorA}; kwargs...) where {TorA} = S_x(TorA, Trivial; kwargs...) S_x(symm::Type{<:Sector}; kwargs...) = S_x(ComplexF64, symm; kwargs...) -function S_x(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) - S_x_mat, _, _ = spinmatrices(spin, elt) +function S_x(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} + S_x_mat, _, _ = spinmatrices(spin, TorA) pspace = ComplexSpace(size(S_x_mat, 1)) return TensorMap(S_x_mat, pspace ← pspace) end -function S_x(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2) +function S_x(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2) where {TorA} spin == 1 // 2 || error("not implemented") pspace = Z2Space(0 => 1, 1 => 1) - X = zeros(elt, pspace, pspace) - block(X, Z2Irrep(0)) .= one(elt) / 2 - block(X, Z2Irrep(1)) .= -one(elt) / 2 + X = zeros(TorA, pspace, pspace) + block(X, Z2Irrep(0)) .= one(eltype(TorA)) / 2 + block(X, Z2Irrep(1)) .= -one(eltype(TorA)) / 2 return X end -function S_x(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) +function S_x(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) where {TorA} pspace = U1Space(i => 1 for i in (-spin):spin) vspace = U1Space(1 => 1, -1 => 1) if side == :L - X = zeros(elt, pspace ← pspace ⊗ vspace) + X = zeros(TorA, pspace ← pspace ⊗ vspace) for (f1, f2) in fusiontrees(X) c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] if c₁.charge == c₂.charge + 1 || c₁.charge + 1 == c₂.charge @@ -86,7 +86,7 @@ function S_x(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) end end elseif side == :R - X = zeros(elt, vspace ⊗ pspace ← pspace) + X = zeros(TorA, vspace ⊗ pspace ← pspace) for (f1, f2) in fusiontrees(X) c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] if c₁.charge == c₂.charge + 1 || c₁.charge + 1 == c₂.charge @@ -118,25 +118,25 @@ See also [`σʸ`](@ref). """ function S_y end S_y(; kwargs...) = S_y(ComplexF64, Trivial; kwargs...) -S_y(elt::Type{<:Complex}; kwargs...) = S_y(elt, Trivial; kwargs...) +S_y(::Type{TorA}; kwargs...) where {TorA} = S_y(elt, Trivial; kwargs...) S_y(symm::Type{<:Sector}; kwargs...) = S_y(ComplexF64, symm; kwargs...) -function S_y(elt::Type{<:Complex}, ::Type{Trivial}; spin = 1 // 2) - _, Y, _, _ = spinmatrices(spin, elt) +function S_y(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} + _, Y, _, _ = spinmatrices(spin, TorA) pspace = ComplexSpace(size(Y, 1)) return TensorMap(Y, pspace ← pspace) end -function S_y(elt::Type{<:Complex}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) +function S_y(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) where {TorA} spin == 1 // 2 || error("not implemented") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(1 => 1) if side == :L - Y = zeros(elt, pspace ← pspace ⊗ vspace) + Y = zeros(TorA, pspace ← pspace ⊗ vspace) block(Y, Z2Irrep(0)) .= one(elt)im / 2 block(Y, Z2Irrep(1)) .= -one(elt)im / 2 elseif side == :R - Y = zeros(elt, vspace ⊗ pspace ← pspace) + Y = zeros(TorA, vspace ⊗ pspace ← pspace) block(Y, Z2Irrep(0)) .= -one(elt)im / 2 block(Y, Z2Irrep(1)) .= one(elt)im / 2 else @@ -145,11 +145,11 @@ function S_y(elt::Type{<:Complex}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) return Y end -function S_y(elt::Type{<:Complex}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) +function S_y(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) where {TorA} pspace = U1Space(i => 1 for i in (-spin):spin) vspace = U1Space(1 => 1, -1 => 1) if side == :L - Y = zeros(elt, pspace ← pspace ⊗ vspace) + Y = zeros(TorA, pspace ← pspace ⊗ vspace) for (f1, f2) in fusiontrees(Y) c₁, c₂ = f1.uncoupled[1], f2.uncoupled[1] if c₁.charge == c₂.charge + 1 @@ -159,7 +159,7 @@ function S_y(elt::Type{<:Complex}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) end end elseif side == :R - Y = zeros(elt, vspace ⊗ pspace ← pspace) + Y = zeros(TorA, vspace ⊗ pspace ← pspace) for (f1, f2) in fusiontrees(Y) c₁, c₂ = f1.uncoupled[2], f2.uncoupled[1] if c₁.charge == c₂.charge + 1 @@ -194,25 +194,25 @@ See also [`σᶻ`](@ref). """ function S_z end S_z(; kwargs...) = S_z(ComplexF64, Trivial; kwargs...) -S_z(elt::Type{<:Number}; kwargs...) = S_z(elt, Trivial; kwargs...) +S_z(::Type{TorA}; kwargs...) where {TorA} = S_z(TorA, Trivial; kwargs...) S_z(symm::Type{<:Sector}; kwargs...) = S_z(ComplexF64, symm; kwargs...) -function S_z(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) - _, _, S_z_mat = spinmatrices(spin, elt) +function S_z(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} + _, _, S_z_mat = spinmatrices(spin, TorA) pspace = ComplexSpace(size(S_z_mat, 1)) return TensorMap(S_z_mat, pspace ← pspace) end -function S_z(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) +function S_z(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) where {TorA} spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(1 => 1) if side == :L - Z = zeros(elt, pspace ← pspace ⊗ vspace) + Z = zeros(TorA, pspace ← pspace ⊗ vspace) block(Z, Z2Irrep(0)) .= one(elt) / 2 block(Z, Z2Irrep(1)) .= one(elt) / 2 elseif side == :R - Z = zeros(elt, vspace ⊗ pspace ← pspace) + Z = zeros(TorA, vspace ⊗ pspace ← pspace) block(Z, Z2Irrep(0)) .= one(elt) / 2 block(Z, Z2Irrep(1)) .= one(elt) / 2 else @@ -221,10 +221,10 @@ function S_z(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) return Z end -function S_z(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2) +function S_z(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2) where {TorA} charges = U1Irrep.((-spin):spin) pspace = U1Space((v => 1 for v in charges)) - Z = zeros(elt, pspace ← pspace) + Z = zeros(TorA, pspace ← pspace) for (c, b) in blocks(Z) b .= c.charge end @@ -250,24 +250,24 @@ See also [`σ⁺`](@ref). """ function S_plus end S_plus(; kwargs...) = S_plus(ComplexF64, Trivial; kwargs...) -S_plus(elt::Type{<:Number}; kwargs...) = S_plus(elt, Trivial; kwargs...) +S_plus(::Type{TorA}; kwargs...) where {TorA} = S_plus(TorA, Trivial; kwargs...) S_plus(symm::Type{<:Sector}; kwargs...) = S_plus(ComplexF64, symm; kwargs...) -function S_plus(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) - S⁺ = S_x(elt, Trivial; spin = spin) + 1im * S_y(complex(elt), Trivial; spin = spin) +function S_plus(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} + S⁺ = S_x(TorA, Trivial; spin = spin) + 1im * S_y(complex(TorA), Trivial; spin = spin) return elt <: Real ? real(S⁺) : S⁺ end -function S_plus(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) +function S_plus(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) where {TorA} spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(0 => 1, 1 => 1) if side == :L - S⁺ = zeros(elt, pspace ← pspace ⊗ vspace) + S⁺ = zeros(TorA, pspace ← pspace ⊗ vspace) block(S⁺, Z2Irrep(0)) .= [1 -1] / 2 block(S⁺, Z2Irrep(1)) .= [-1 1] / 2 elseif side == :R - S⁺ = zeros(elt, vspace ⊗ pspace ← pspace) + S⁺ = zeros(TorA, vspace ⊗ pspace ← pspace) block(S⁺, Z2Irrep(0)) .= [1 1]' / 2 block(S⁺, Z2Irrep(1)) .= [-1 -1]' / 2 else @@ -276,17 +276,17 @@ function S_plus(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) return S⁺ end -function S_plus(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) +function S_plus(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) where {TorA} pspace = U1Space(i => 1 for i in (-spin):spin) if side == :L vspace = U1Space(1 => 1) - S⁺ = zeros(elt, pspace ← pspace ⊗ vspace) + S⁺ = zeros(TorA, pspace ← pspace ⊗ vspace) for (c, b) in blocks(S⁺) b .= 2 * _pauliterm(spin, c, only(c ⊗ U1Irrep(-1))) end elseif side == :R vspace = U1Space(-1 => 1) - S⁺ = zeros(elt, vspace ⊗ pspace ← pspace) + S⁺ = zeros(TorA, vspace ⊗ pspace ← pspace) for (c, b) in blocks(S⁺) b .= 2 * _pauliterm(spin, only(c ⊗ U1Irrep(+1)), c) end @@ -315,24 +315,24 @@ See also [`σ⁻`](@ref). """ function S_min end S_min(; kwargs...) = S_min(ComplexF64, Trivial; kwargs...) -S_min(elt::Type{<:Number}; kwargs...) = S_min(elt, Trivial; kwargs...) +S_min(::Type{TorA}; kwargs...) where {TorA} = S_min(TorA, Trivial; kwargs...) S_min(symm::Type{<:Sector}; kwargs...) = S_min(ComplexF64, symm; kwargs...) -function S_min(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) - S⁻ = S_x(elt, Trivial; spin = spin) - 1im * S_y(complex(elt), Trivial; spin = spin) - return elt <: Real ? real(S⁻) : S⁻ +function S_min(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} + S⁻ = S_x(TorA, Trivial; spin = spin) - 1im * S_y(complex(TorA), Trivial; spin = spin) + return eltype(TorA) <: Real ? real(S⁻) : S⁻ end -function S_min(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) +function S_min(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) where {TorA} spin == 1 // 2 || error("Z2 symmetry only implemented for spin 1 // 2") pspace = Z2Space(0 => 1, 1 => 1) vspace = Z2Space(0 => 1, 1 => 1) if side == :L - S⁻ = zeros(elt, pspace ← pspace ⊗ vspace) + S⁻ = zeros(TorA, pspace ← pspace ⊗ vspace) block(S⁻, Z2Irrep(0)) .= [1 1] / 2 block(S⁻, Z2Irrep(1)) .= [-1 -1] / 2 elseif side == :R - S⁻ = zeros(elt, vspace ⊗ pspace ← pspace) + S⁻ = zeros(TorA, vspace ⊗ pspace ← pspace) block(S⁻, Z2Irrep(0)) .= [1 -1]' / 2 block(S⁻, Z2Irrep(1)) .= [1 -1]' / 2 else @@ -341,17 +341,17 @@ function S_min(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2, side = :L) return S⁻ end -function S_min(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) +function S_min(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2, side = :L) where {TorA} pspace = U1Space(i => 1 for i in (-spin):spin) if side == :L vspace = U1Space(-1 => 1) - S⁻ = zeros(elt, pspace ← pspace ⊗ vspace) + S⁻ = zeros(TorA, pspace ← pspace ⊗ vspace) for (c, b) in blocks(S⁻) b .= 2 * _pauliterm(spin, only(c ⊗ U1Irrep(+1)), c) end elseif side == :R vspace = U1Space(+1 => 1) - S⁻ = zeros(elt, vspace ⊗ pspace ← pspace) + S⁻ = zeros(TorA, vspace ⊗ pspace ← pspace) for (c, b) in blocks(S⁻) b .= 2 * _pauliterm(spin, c, only(c ⊗ U1Irrep(-1))) end @@ -404,17 +404,17 @@ for (L, R) in ((:x, :x), (:y, :y), (:z, :z), (:plus, :min), (:min, :plus)) ($f)(elt::Type{<:Number}; kwargs...) = ($f)(elt, Trivial; kwargs...) ($f)(symmetry::Type{<:Sector}; kwargs...) = ($f)(ComplexF64, symmetry; kwargs...) - function ($f)(elt::Type{<:Number}, ::Type{Trivial}; spin = 1 // 2) + function ($f)(::Type{TorA}, ::Type{Trivial}; spin = 1 // 2) where {TorA} return contract_twosite( - $(fₗ)(elt, Trivial; spin = spin), - $(fᵣ)(elt, Trivial; spin = spin) + $(fₗ)(TorA, Trivial; spin = spin), + $(fᵣ)(TorA, Trivial; spin = spin) ) end - function ($f)(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin = 1 // 2) + function ($f)(::Type{TorA}, symmetry::Type{<:Sector}; spin = 1 // 2) where {TorA} return contract_twosite( - $(fₗ)(elt, symmetry; spin = spin, side = :L), - $(fᵣ)(elt, symmetry; spin = spin, side = :R) + $(fₗ)(TorA, symmetry; spin = spin, side = :L), + $(fᵣ)(TorA, symmetry; spin = spin, side = :R) ) end @@ -426,11 +426,11 @@ for (L, R) in ((:x, :x), (:y, :y), (:z, :z), (:plus, :min), (:min, :plus)) end end -function S_xx(elt::Type{<:Number}, ::Type{Z2Irrep}; spin = 1 // 2) - return contract_twosite(S_x(elt, Z2Irrep; spin = spin), S_x(elt, Z2Irrep; spin = spin)) +function S_xx(::Type{TorA}, ::Type{Z2Irrep}; spin = 1 // 2) where {TorA} + return contract_twosite(S_x(TorA, Z2Irrep; spin = spin), S_x(TorA, Z2Irrep; spin = spin)) end -function S_zz(elt::Type{<:Number}, ::Type{U1Irrep}; spin = 1 // 2) - return contract_twosite(S_z(elt, U1Irrep; spin = spin), S_z(elt, U1Irrep; spin = spin)) +function S_zz(::Type{TorA}, ::Type{U1Irrep}; spin = 1 // 2) where {TorA} + return contract_twosite(S_z(TorA, U1Irrep; spin = spin), S_z(TorA, U1Irrep; spin = spin)) end """ @@ -443,25 +443,25 @@ See also [`σσ`](@ref). """ function S_exchange end S_exchange(; kwargs...) = S_exchange(ComplexF64, Trivial; kwargs...) -S_exchange(elt::Type{<:Number}; kwargs...) = S_exchange(elt, Trivial; kwargs...) +S_exchange(::Type{TorA}; kwargs...) where {TorA} = S_exchange(TorA, Trivial; kwargs...) function S_exchange(symmetry::Type{<:Sector}; kwargs...) return S_exchange(ComplexF64, symmetry; kwargs...) end -function S_exchange(elt::Type{<:Number}, symmetry::Type{<:Sector}; spin = 1 // 2) - elt_complex = complex(elt) +function S_exchange(::Type{TorA}, symmetry::Type{<:Sector}; spin = 1 // 2) where {TorA} + TorA_complex = TensorKit.similarstoragetype(TorA, complex(eltype(TorA))) SS = ( - S_plusmin(elt_complex, symmetry; spin = spin) + - S_minplus(elt_complex, symmetry; spin = spin) + S_plusmin(TorA_complex, symmetry; spin = spin) + + S_minplus(TorA_complex, symmetry; spin = spin) ) / 2 + - S_zz(elt_complex, symmetry; spin = spin) - return elt <: Real ? real(SS) : SS + S_zz(TorA_complex, symmetry; spin = spin) + return eltype(TorA) <: Real ? real(SS) : SS end -function S_exchange(elt::Type{<:Number}, ::Type{SU2Irrep}; spin = 1 // 2) +function S_exchange(::Type{TorA}, ::Type{SU2Irrep}; spin = 1 // 2) where {TorA} pspace = SU2Space(spin => 1) aspace = SU2Space(1 => 1) - Sleft = ones(elt, pspace ← pspace ⊗ aspace) - Sright = -ones(elt, aspace ⊗ pspace ← pspace) + Sleft = ones(TorA, pspace ← pspace ⊗ aspace) + Sright = -ones(TorA, aspace ⊗ pspace ← pspace) @tensor SS[-1 -2; -3 -4] := Sleft[-1; -3 1] * Sright[1 -2; -4] * (spin^2 + spin) return SS @@ -483,23 +483,23 @@ The Potts operator ``Z ⊗ Z'``, where ``Z^q = 1``. """ function potts_ZZ end potts_ZZ(; kwargs...) = potts_ZZ(ComplexF64, Trivial; kwargs...) -potts_ZZ(elt::Type{<:Number}; kwargs...) = potts_ZZ(elt, Trivial; kwargs...) +potts_ZZ(::Type{TorA}; kwargs...) where {TorA} = potts_ZZ(TorA, Trivial; kwargs...) function potts_ZZ(symmetry::Type{<:Sector}; kwargs...) return potts_ZZ(ComplexF64, symmetry; kwargs...) end -function potts_ZZ(elt::Type{<:Number}, ::Type{Trivial}; q = 3) - Z = potts_Z(elt, Trivial; q = q) +function potts_ZZ(::Type{TorA}, ::Type{Trivial}; q = 3) where {TorA} + Z = potts_Z(TorA, Trivial; q = q) return Z ⊗ Z' end -function potts_ZZ(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q = Q) where {Q} +function potts_ZZ(::Type{TorA}, ::Type{ZNIrrep{Q}}; q = Q) where {Q, TorA} @assert q == Q "q must match the irrep charge" pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1)) - ZZ = zeros(elt, pspace ⊗ pspace ← pspace ⊗ pspace) + ZZ = zeros(TorA, pspace ⊗ pspace ← pspace ⊗ pspace) for charge in 0:(Q - 1) for i in 1:Q - block(ZZ, ZNIrrep{Q}(charge))[i, mod1(i + 1, Q)] = one(elt) + block(ZZ, ZNIrrep{Q}(charge))[i, mod1(i + 1, Q)] .= one(eltype(TorA)) end end return ZZ @@ -510,15 +510,15 @@ end the Weyl-Heisenberg matrices according to [Wikipedia](https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices#Sylvester's_generalized_Pauli_matrices_(non-Hermitian)). """ -function weyl_heisenberg_matrices(Q::Int, elt = ComplexF64) - U = zeros(elt, Q, Q) # clock matrix - V = zeros(elt, Q, Q) # shift matrix - W = zeros(elt, Q, Q) # DFT +function weyl_heisenberg_matrices(Q::Int, ::Type{TorA}) where {TorA} + U = zeros(TorA, Q, Q) # clock matrix + V = zeros(TorA, Q, Q) # shift matrix + W = zeros(TorA, Q, Q) # DFT ω = cis(2 * pi / Q) for row in 1:Q U[row, row] = ω^(row - 1) - V[row, mod1(row - 1, Q)] = one(elt) + V[row, mod1(row - 1, Q)] .= one(eltype(TorA)) for col in 1:Q W[row, col] = ω^((row - 1) * (col - 1)) end @@ -533,10 +533,10 @@ The Potts ``Z`` operator, also known as the clock operator, where ``Z^q = 1``. """ function potts_Z end potts_Z(; kwargs...) = potts_Z(ComplexF64, Trivial; kwargs...) -potts_Z(elt::Type{<:Complex}; kwargs...) = potts_Z(elt, Trivial; kwargs...) +potts_Z(::Type{TorA}; kwargs...) where {TorA} = potts_Z(TorA, Trivial; kwargs...) potts_Z(symm::Type{<:Sector}; kwargs...) = potts_Z(ComplexF64, symm; kwargs...) -function potts_Z(elt::Type{<:Number}, ::Type{Trivial}; q = 3) - U, _, _ = weyl_heisenberg_matrices(q, elt) +function potts_Z(::Type{TorA}, ::Type{Trivial}; q = 3) where {TorA} + U, _, _ = weyl_heisenberg_matrices(q, TorA) Z = TensorMap(U, ComplexSpace(q) ← ComplexSpace(q)) return Z end @@ -549,18 +549,18 @@ The Potts ``X`` operator, also known as the shift operator, where ``X^q = 1``. """ function potts_X end potts_X(; kwargs...) = potts_X(ComplexF64, Trivial; kwargs...) -potts_X(elt::Type{<:Complex}; kwargs...) = potts_X(elt, Trivial; kwargs...) +potts_X(::Type{TorA}; kwargs...) where {TorA} = potts_X(TorA, Trivial; kwargs...) potts_X(symm::Type{<:Sector}; kwargs...) = potts_X(ComplexF64, symm; kwargs...) -function potts_X(elt::Type{<:Number}, ::Type{Trivial}; q = 3) - _, V, _ = weyl_heisenberg_matrices(q, elt) +function potts_X(::Type{TorA}, ::Type{Trivial}; q = 3) where {TorA} + _, V, _ = weyl_heisenberg_matrices(q, TorA) X = TensorMap(V, ComplexSpace(q) ← ComplexSpace(q)) return X end -function potts_X(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q = Q) where {Q} +function potts_X(::Type{TorA}, ::Type{ZNIrrep{Q}}; q = Q) where {Q, TorA} @assert q == Q "q must match the irrep charge" pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1)) - X = zeros(elt, pspace ← pspace) + X = zeros(TorA, pspace ← pspace) ω = cis(2 * pi / Q) for i in 1:Q block(X, ZNIrrep{Q}(i)) .= ω^i