From 05de79db3357390ea5184ab896659cac8499886b Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 14 Aug 2026 18:58:04 +0800 Subject: [PATCH 1/8] Create custom operators from plain Arrays --- docs/src/index.md | 16 +++++++ src/bosonoperators.jl | 18 +++++++- src/fermionoperators.jl | 14 +++++- src/hubbardoperators.jl | 28 +++++++++--- src/quantumgates.jl | 14 +++++- src/spinoperators.jl | 17 +++++++- src/tjoperators.jl | 26 ++++++++++- src/utils.jl | 32 ++++++++++++++ test/custom.jl | 97 +++++++++++++++++++++++++++++++++++++++++ 9 files changed, 252 insertions(+), 10 deletions(-) create mode 100644 test/custom.jl diff --git a/docs/src/index.md b/docs/src/index.md index 8bd6e61..04cfe18 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -22,6 +22,22 @@ eigenbasis for `SpinOperators`), and the projection is available as [`symmetrize symmetrizing custom operators. Operators that are incompatible with a given symmetry throw an `ArgumentError`. +### Custom operators + +Each operator module exports its own `custom` function as the high-level entry point for constructing custom operators from plain arrays. +Import the module for the relevant physical system, then pass the desired symmetry types explicitly: + +```julia +using TensorKit +using TensorKitTensors.SpinOperators +A = [1.0 0.0; 0.0 -1.0] +Z = custom(A, U1Irrep) +``` + +An `N`-site input is a rank-`2N` array ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. +Hubbard and t-J operators take separate particle and spin symmetries. +For t-J, `slave_fermion=true` transforms an array supplied in the normal t-J reference basis into the slave-fermion representation. + ```@docs symmetrize desymmetrize diff --git a/src/bosonoperators.jl b/src/bosonoperators.jl index b8a42d5..5764625 100644 --- a/src/bosonoperators.jl +++ b/src/bosonoperators.jl @@ -3,8 +3,9 @@ module BosonOperators using TensorKit using LinearAlgebra: I import ..TensorKitTensors: symmetrize, desymmetrize, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export boson_space, basis_transform +export boson_space, basis_transform, custom export b_plus, b_min, b_num export b_plus_b_plus, b_plus_b_min, b_min_b_plus, b_min_b_min export b_hopping @@ -62,6 +63,21 @@ end _symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = symmetrize(O, basis_transform(symmetry; kwargs...), boson_space(symmetry; kwargs...)) +""" + custom(A::AbstractArray, symmetry::Type{<:Sector}; cutoff::Integer, tol=nothing) + +Construct a symmetry-aware bosonic operator from a dense rank-`2N` array in the occupation-number basis. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and each axis must have dimension `cutoff + 1`. +""" +function custom( + A::AbstractArray, symmetry::Type{<:Sector}; cutoff::Integer, tol = nothing + ) + O = _check_custom_space(_custom_dense_operator(A), boson_space(Trivial; cutoff)) + return symmetrize( + O, basis_transform(symmetry; cutoff), boson_space(symmetry; cutoff); tol + ) +end + # Single-site operators # --------------------- """ diff --git a/src/fermionoperators.jl b/src/fermionoperators.jl index 2269f96..f147931 100644 --- a/src/fermionoperators.jl +++ b/src/fermionoperators.jl @@ -3,8 +3,9 @@ module FermionOperators using TensorKit using LinearAlgebra: I import ..TensorKitTensors: symmetrize, desymmetrize, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export fermion_space, basis_transform +export fermion_space, basis_transform, custom export f_num export f_plus_f_min, f_min_f_plus, f_plus_f_plus, f_min_f_min export f_hopping @@ -59,6 +60,17 @@ end _symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}) = symmetrize(O, basis_transform(symmetry), fermion_space(symmetry)) +""" + custom(A::AbstractArray, symmetry::Type{<:Sector}; tol=nothing) + +Construct a symmetry-aware spinless-fermion operator from a dense rank-`2N` array in the basis ``|0⟩, |1⟩``. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. Fermion-parity grading is enforced during the construction. +""" +function custom(A::AbstractArray, symmetry::Type{<:Sector}; tol = nothing) + O = _check_custom_space(_custom_dense_operator(A), fermion_space(Trivial)) + return symmetrize(O, basis_transform(symmetry), fermion_space(symmetry); tol) +end + # Single-site operators # --------------------- function single_site_operator(elt::Type{<:Number}, symmetry::Type{<:Sector}) diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index 437cc46..c124cc6 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -2,8 +2,9 @@ module HubbardOperators using TensorKit import ..TensorKitTensors: symmetrize, desymmetrize, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export hubbard_space, basis_transform +export hubbard_space, basis_transform, custom export e_num, u_num, d_num, ud_num, half_ud_num, h_num export S_x, S_y, S_z, S_plus, S_min export u_plus_u_min, d_plus_d_min @@ -166,11 +167,11 @@ const _PARTICLE_GAUGE = Complex{Int}[1 0 0 0; 0 -1 0 0; 0 0 im 0; 0 0 0 im] # If the staggered gauge commutes we don't incorporate it to retain the option for real tensors. function _symmetrize_hubbard( O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector} + spin_symmetry::Type{<:Sector}; tol = nothing ) U = basis_transform(particle_symmetry, spin_symmetry) V = hubbard_space(particle_symmetry, spin_symmetry) - particle_symmetry === SU2Irrep || return symmetrize(O, U, V) + particle_symmetry === SU2Irrep || return symmetrize(O, U, V; tol) Vref = domain(U)[1] Gs = ntuple(k -> TensorMap(_PARTICLE_GAUGE^(k - 1), Vref ← Vref), numout(O)) @@ -178,17 +179,34 @@ function _symmetrize_hubbard( Od = desymmetrize(O) # check if commutes with the staggered gauge - W * Od ≈ Od * W && return symmetrize(O, U, V) + W * Od ≈ Od * W && return symmetrize(O, U, V; tol) scalartype(O) <: Real && throw(ArgumentError("operator with `SU2Irrep` particle symmetry that does not commute with the staggered gauge requires a complex scalar type")) - return symmetrize(O, map(g -> U * g, Gs), V) + return symmetrize(O, map(g -> U * g, Gs), V; tol) end # Symmetrize a Hubbard operator through its (staggered-gauge) basis transformation _symmetrize_operator(O::AbstractTensorMap, P::Type{<:Sector}, S::Type{<:Sector}) = _symmetrize_hubbard(O, P, S) +""" + custom(A::AbstractArray, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; tol=nothing) + +Construct a symmetry-aware Hubbard operator from a dense rank-`2N` array in the reference basis ``|0⟩, |↑↓⟩, |↑⟩, |↓⟩``. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. +Fermion parity and the staggered gauge required by `SU2Irrep` particle symmetry are applied automatically. +""" +function custom( + A::AbstractArray, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}; tol = nothing + ) + O = _check_custom_space( + _custom_dense_operator(A), hubbard_space(Trivial, Trivial) + ) + return _symmetrize_hubbard(O, particle_symmetry, spin_symmetry; tol) +end + function n_site_operator(::Val{N}, elt::Type{<:Number}) where {N} V = hubbard_space(Trivial, Trivial) return zeros(elt, V^N ← V^N) diff --git a/src/quantumgates.jl b/src/quantumgates.jl index 2b0feef..da82ae7 100644 --- a/src/quantumgates.jl +++ b/src/quantumgates.jl @@ -4,8 +4,9 @@ using TensorKit using LinearAlgebra: I using RationalRoots: signedroot import ..TensorKitTensors: symmetrize, desymmetrize, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export qubit_space, basis_transform, controlled +export qubit_space, basis_transform, custom, controlled export pauli_x, pauli_y, pauli_z, proj_0, proj_1, hadamard, s_gate, t_gate export phase_shift, rotation_x, rotation_y, rotation_z, u3 export cnot, cy, cz, cphase, ch, cs, swap, iswap, dcx, ecr @@ -61,6 +62,17 @@ end _symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = symmetrize(O, basis_transform(symmetry), qubit_space(symmetry)) +""" + custom(A::AbstractArray, symmetry::Type{<:Sector}; tol=nothing) + +Construct a symmetry-aware qubit operator from a dense rank-`2N` array in the computational basis. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and every axis must have dimension two. +""" +function custom(A::AbstractArray, symmetry::Type{<:Sector}; tol = nothing) + O = _check_custom_space(_custom_dense_operator(A), qubit_space(Trivial)) + return symmetrize(O, basis_transform(symmetry), qubit_space(symmetry); tol) +end + # Gate combinators # ---------------- """ diff --git a/src/spinoperators.jl b/src/spinoperators.jl index 4530f94..7600186 100644 --- a/src/spinoperators.jl +++ b/src/spinoperators.jl @@ -4,8 +4,9 @@ using TensorKit using LinearAlgebra: I using RationalRoots: signedroot import ..TensorKitTensors: symmetrize, desymmetrize, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export spin_space, basis_transform, casimir +export spin_space, basis_transform, custom, casimir export S_x, S_y, S_z, S_plus, S_min export S_x_S_x, S_y_S_y, S_z_S_z, S_plus_S_min, S_min_S_plus, S_exchange export σˣ, σʸ, σᶻ, σ⁺, σ⁻ @@ -89,6 +90,20 @@ end _symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = symmetrize(O, basis_transform(symmetry; kwargs...), spin_space(symmetry; kwargs...)) +""" + custom(A::AbstractArray, symmetry::Type{<:Sector}; spin=1 // 2, tol=nothing) + +Construct a symmetry-aware spin operator from a dense rank-`2N` array in the standard spin basis. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and each axis must have dimension `2spin + 1`. +""" +function custom( + A::AbstractArray, symmetry::Type{<:Sector}; spin = 1 // 2, tol = nothing + ) + V = spin_space(symmetry; spin) + O = _check_custom_space(_custom_dense_operator(A), spin_space(Trivial; spin)) + return symmetrize(O, basis_transform(symmetry; spin), V; tol) +end + # Pauli matrices # -------------- function _pauliterm(spin, i, j) diff --git a/src/tjoperators.jl b/src/tjoperators.jl index 3d6297a..4d4fc68 100644 --- a/src/tjoperators.jl +++ b/src/tjoperators.jl @@ -4,8 +4,9 @@ using LinearAlgebra: diagind using TensorKit import ..HubbardOperators import ..TensorKitTensors: symmetrize, desymmetrize, fuse_local_operators, @operator +import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export tj_space, basis_transform, tj_projector +export tj_space, basis_transform, custom, tj_projector export transform_slave_fermion # the operator names and their aliases are exported by the generation loop at the bottom of # this file, from the `_OPERATORS` registry @@ -166,6 +167,29 @@ function _symmetrize_operator( ) end +""" + custom(A::AbstractArray, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false, tol=nothing) + +Construct a symmetry-aware t-J operator from a dense rank-`2N` array in the normal reference basis ``|0⟩, |↑⟩, |↓⟩``. +The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. +Fermion parity is enforced, and `slave_fermion=true` applies the full slave-fermion transformation, before imposing the requested symmetries. +""" +function custom( + A::AbstractArray, particle_symmetry::Type{<:Sector}, + spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false, tol = nothing + ) + O = _check_custom_space(_custom_dense_operator(A), tj_space(Trivial, Trivial)) + if slave_fermion + V = tj_space(Trivial, Trivial) + U = basis_transform(Trivial, Trivial) + O = transform_slave_fermion(symmetrize(O, U, V; tol)) + end + return symmetrize( + O, basis_transform(particle_symmetry, spin_symmetry; slave_fermion), + tj_space(particle_symmetry, spin_symmetry; slave_fermion); tol + ) +end + # Relation to the Hubbard model # ----------------------------- """ diff --git a/src/utils.jl b/src/utils.jl index e10e0e1..f84a616 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -129,6 +129,38 @@ function _default_tol(::Type{T}, I::Type{<:Sector}) where {T <: Number} return sqrt(Tsector <: AbstractFloat ? max(ε, eps(Tsector)) : ε) end +# Convert the common dense-array representation of a square N-site operator to a TensorMap +# over ComplexSpace. The individual operator modules then check the inferred local dimension +# and project this map onto their own symmetric (and, where applicable, fermion-graded) spaces. +function _custom_dense_operator(A::AbstractArray) + eltype(A) <: Number || throw(ArgumentError("operator array must have a numeric element type")) + rank = ndims(A) + rank >= 2 && iseven(rank) || + throw(ArgumentError("operator array must have positive even rank `2N`, got rank $rank")) + local_dimension = size(A, 1) + expected_size = ntuple(Returns(local_dimension), rank) + size(A) == expected_size || throw( + ArgumentError( + "all operator-array legs must have the same size, got $(size(A))" + ) + ) + sites = rank ÷ 2 + V = ComplexSpace(local_dimension) + return TensorMap(A, V^sites ← V^sites) +end + +function _check_custom_space(O::AbstractTensorMap, V::ElementarySpace) + n = numin(O) + Vdense = desymmetrize(V) + expected_space = Vdense^n ← Vdense^n + space(O) == expected_space || throw( + ArgumentError( + "operator array has space `$(space(O))`, expected `$expected_space`" + ) + ) + return O +end + """ fuse_local_operators(O₁, O₂) diff --git a/test/custom.jl b/test/custom.jl new file mode 100644 index 0000000..2d613e0 --- /dev/null +++ b/test/custom.jl @@ -0,0 +1,97 @@ +using TensorKit +using Test +import TensorKitTensors +using TensorKitTensors: desymmetrize +import TensorKitTensors.SpinOperators as SO +import TensorKitTensors.BosonOperators as BO +import TensorKitTensors.FermionOperators as FO +import TensorKitTensors.HubbardOperators as HO +import TensorKitTensors.TJOperators as TJ +import TensorKitTensors.QuantumGates as QG + +@testset "SpinOperators" begin + Sz = convert(Array, SO.S_z()) + @test SO.custom(Sz, U1Irrep) ≈ SO.S_z(U1Irrep) + + Sz_spin1 = SO.S_z(Float64, Trivial; spin = 1) + @test SO.custom(convert(Array, Sz_spin1), U1Irrep; spin = 1) ≈ + SO.S_z(Float64, U1Irrep; spin = 1) + + # rank determines the number of sites and the scalar type is preserved + SS_big = SO.S_exchange(Complex{BigFloat}) + SS_custom = SO.custom(convert(Array, SS_big), SU2Irrep) + @test numout(SS_custom) == numin(SS_custom) == 2 + @test scalartype(SS_custom) === Complex{BigFloat} + @test SS_custom ≈ SO.S_exchange(Complex{BigFloat}, SU2Irrep) + + # tol is used only for symmetry projection + X_pert = copy(convert(Array, SO.S_x())) + X_pert[1, 1] += 1.0e-6 + @test_throws ArgumentError SO.custom(X_pert, Z2Irrep) + @test SO.custom(X_pert, Z2Irrep; tol = 1.0e-3) ≈ + SO.S_x(Z2Irrep) atol = 1.0e-5 + + # common dense-array validation + @test_throws ArgumentError SO.custom(fill("x", 2, 2), Trivial) + @test_throws ArgumentError SO.custom(zeros(2), Trivial) + @test_throws ArgumentError SO.custom(fill(1.0), Trivial) + @test_throws ArgumentError SO.custom(zeros(2, 3), Trivial) + @test_throws ArgumentError SO.custom(zeros(4, 4), Trivial) + @test_throws ArgumentError SO.custom(Sz, FermionParity) +end + +@testset "BosonOperators" begin + cutoff = 3 + N = BO.b_num(Float64, Trivial; cutoff) + @test BO.custom(convert(Array, N), U1Irrep; cutoff) ≈ + BO.b_num(Float64, U1Irrep; cutoff) + + hopping = BO.b_hopping(Float64, Trivial; cutoff) + @test BO.custom(convert(Array, hopping), U1Irrep; cutoff) ≈ + BO.b_hopping(Float64, U1Irrep; cutoff) +end + +@testset "FermionOperators" begin + hopping = FO.f_hopping(ComplexF64, Trivial) + A = convert(Array, desymmetrize(hopping)) + @test FO.custom(A, Trivial) ≈ hopping + @test FO.custom(A, U1Irrep) ≈ + FO.f_hopping(ComplexF64, U1Irrep) + + # `Trivial` still enforces the mandatory fermion-parity grading + parity_odd = ComplexF64[0 1; 0 0] + @test_throws ArgumentError FO.custom(parity_odd, Trivial) +end + +@testset "HubbardOperators" begin + hopping = HO.e_hopping(ComplexF64, Trivial, Trivial) + A = convert(Array, desymmetrize(hopping)) + @test HO.custom(A, SU2Irrep, SU2Irrep) ≈ + HO.e_hopping(ComplexF64, SU2Irrep, SU2Irrep) + + hopping_real = HO.e_hopping(Float64, Trivial, Trivial) + @test_throws ArgumentError HO.custom( + convert(Array, desymmetrize(hopping_real)), SU2Irrep, SU2Irrep + ) + + # test an operator with more than two-sites + paircor = HO.Δ⁺ij_Δkl(ComplexF64, Trivial, Trivial) + A = convert(Array, desymmetrize(paircor)) + @test HO.custom(A, U1Irrep, SU2Irrep) ≈ + HO.Δ⁺ij_Δkl(ComplexF64, U1Irrep, SU2Irrep) +end + +@testset "TJOperators" begin + hopping = TJ.e_hopping(ComplexF64, Trivial, Trivial) + A = convert(Array, desymmetrize(hopping)) + @test TJ.custom(A, U1Irrep, SU2Irrep; slave_fermion = true) ≈ + TJ.e_hopping(ComplexF64, U1Irrep, SU2Irrep; slave_fermion = true) +end + +@testset "QuantumGates" begin + Z = convert(Array, QG.pauli_z()) + @test QG.custom(Z, U1Irrep) ≈ QG.pauli_z(U1Irrep) + + CZ = convert(Array, QG.cz()) + @test QG.custom(CZ, U1Irrep) ≈ QG.cz(U1Irrep) +end From 1d82549b6920e0e2756f9f9bd1cc07ee48ea7bf6 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Fri, 14 Aug 2026 22:07:11 +0800 Subject: [PATCH 2/8] Directly export `symmetrize_operator` --- docs/src/index.md | 23 ++++---- src/bosonoperators.jl | 18 ++---- src/fermionoperators.jl | 18 +++--- src/hubbardoperators.jl | 36 +++--------- src/quantumgates.jl | 19 +++--- src/spinoperators.jl | 17 ++---- src/tjoperators.jl | 43 ++++++-------- src/utils.jl | 40 ++----------- test/custom.jl | 97 ------------------------------- test/symmetrize_operator.jl | 111 ++++++++++++++++++++++++++++++++++++ 10 files changed, 177 insertions(+), 245 deletions(-) delete mode 100644 test/custom.jl create mode 100644 test/symmetrize_operator.jl diff --git a/docs/src/index.md b/docs/src/index.md index 04cfe18..e7db5db 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -14,29 +14,28 @@ Pages = Main.operatorpages ## Symmetric operators through basis transformations Each operator module defines its operators only once, in a non-symmetric reference basis. -The symmetric versions are generated automatically by rotating the reference operator with a -documented unitary basis transformation and projecting the result onto the symmetric tensor -structure. The basis transformations are exposed through each module's `basis_transform` -function (e.g. the Hadamard matrix that maps the ``S^z`` basis onto the ``ℤ₂`` spin-flip -eigenbasis for `SpinOperators`), and the projection is available as [`symmetrize`](@ref) for -symmetrizing custom operators. Operators that are incompatible with a given symmetry throw -an `ArgumentError`. +The symmetric versions are generated automatically by rotating the reference operator with a documented unitary basis transformation and projecting the result onto the symmetric tensor structure. +The basis transformations are exposed through each module's `basis_transform` function, and each module exports `symmetrize_operator` to apply the appropriate transformation and symmetry projection to an operator on its desymmetrized form. +Operators that are incompatible with a given symmetry throw an `ArgumentError`. ### Custom operators -Each operator module exports its own `custom` function as the high-level entry point for constructing custom operators from plain arrays. -Import the module for the relevant physical system, then pass the desired symmetry types explicitly: +Construct a `TensorMap` on the relevant module's `Trivial` reference space, then pass it to that module's `symmetrize_operator` with the desired symmetry types: ```julia using TensorKit using TensorKitTensors.SpinOperators A = [1.0 0.0; 0.0 -1.0] -Z = custom(A, U1Irrep) +V = spin_space(Trivial) +sites = ndims(A) ÷ 2 +O = TensorMap(A, V^sites ← V^sites) +Z = symmetrize_operator(O, U1Irrep) ``` -An `N`-site input is a rank-`2N` array ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. +The caller is responsible for constructing a square homogeneous-site operator with the intended local space and array-axis ordering. +For fermionic systems, constructing the input on the module's `Trivial` space preserves the mandatory fermion-parity grading. Hubbard and t-J operators take separate particle and spin symmetries. -For t-J, `slave_fermion=true` transforms an array supplied in the normal t-J reference basis into the slave-fermion representation. +For t-J, `slave_fermion=true` transforms an operator supplied in the normal t-J representation into the slave-fermion representation. ```@docs symmetrize diff --git a/src/bosonoperators.jl b/src/bosonoperators.jl index 5764625..cef2b93 100644 --- a/src/bosonoperators.jl +++ b/src/bosonoperators.jl @@ -3,9 +3,8 @@ module BosonOperators using TensorKit using LinearAlgebra: I import ..TensorKitTensors: symmetrize, desymmetrize, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export boson_space, basis_transform, custom +export boson_space, basis_transform, symmetrize_operator export b_plus, b_min, b_num export b_plus_b_plus, b_plus_b_min, b_min_b_plus, b_min_b_min export b_hopping @@ -59,20 +58,15 @@ function basis_transform(symmetry::Type{<:Sector}; cutoff::Integer) return TensorMap(Matrix{Int}(I, cutoff + 1, cutoff + 1), V ← boson_space(Trivial; cutoff)) end -# Symmetrize a boson operator through its basis transformation -_symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = - symmetrize(O, basis_transform(symmetry; kwargs...), boson_space(symmetry; kwargs...)) - """ - custom(A::AbstractArray, symmetry::Type{<:Sector}; cutoff::Integer, tol=nothing) + symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; cutoff::Integer, tol=nothing) -Construct a symmetry-aware bosonic operator from a dense rank-`2N` array in the occupation-number basis. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and each axis must have dimension `cutoff + 1`. +Symmetrize a bosonic operator defined on `boson_space(Trivial; cutoff)` through the basis transformation for `symmetry`. """ -function custom( - A::AbstractArray, symmetry::Type{<:Sector}; cutoff::Integer, tol = nothing +function symmetrize_operator( + O::AbstractTensorMap, symmetry::Type{<:Sector}; cutoff::Integer, + tol = nothing ) - O = _check_custom_space(_custom_dense_operator(A), boson_space(Trivial; cutoff)) return symmetrize( O, basis_transform(symmetry; cutoff), boson_space(symmetry; cutoff); tol ) diff --git a/src/fermionoperators.jl b/src/fermionoperators.jl index f147931..bf0e12b 100644 --- a/src/fermionoperators.jl +++ b/src/fermionoperators.jl @@ -3,9 +3,8 @@ module FermionOperators using TensorKit using LinearAlgebra: I import ..TensorKitTensors: symmetrize, desymmetrize, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export fermion_space, basis_transform, custom +export fermion_space, basis_transform, symmetrize_operator export f_num export f_plus_f_min, f_min_f_plus, f_plus_f_plus, f_min_f_min export f_hopping @@ -56,18 +55,15 @@ function basis_transform(symmetry::Type{<:Sector}) return TensorMap(Matrix{Int}(I, 2, 2), V ← desymmetrize(fermion_space(Trivial))) end -# Symmetrize a fermion operator through its basis transformation -_symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}) = - symmetrize(O, basis_transform(symmetry), fermion_space(symmetry)) - """ - custom(A::AbstractArray, symmetry::Type{<:Sector}; tol=nothing) + symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; tol=nothing) -Construct a symmetry-aware spinless-fermion operator from a dense rank-`2N` array in the basis ``|0⟩, |1⟩``. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. Fermion-parity grading is enforced during the construction. +Symmetrize a spinless-fermion operator defined on `fermion_space(Trivial)` through the basis transformation for `symmetry`. +The input space must retain the mandatory fermion-parity grading. """ -function custom(A::AbstractArray, symmetry::Type{<:Sector}; tol = nothing) - O = _check_custom_space(_custom_dense_operator(A), fermion_space(Trivial)) +function symmetrize_operator( + O::AbstractTensorMap, symmetry::Type{<:Sector}; tol = nothing + ) return symmetrize(O, basis_transform(symmetry), fermion_space(symmetry); tol) end diff --git a/src/hubbardoperators.jl b/src/hubbardoperators.jl index c124cc6..c9f5426 100644 --- a/src/hubbardoperators.jl +++ b/src/hubbardoperators.jl @@ -2,9 +2,8 @@ module HubbardOperators using TensorKit import ..TensorKitTensors: symmetrize, desymmetrize, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export hubbard_space, basis_transform, custom +export hubbard_space, basis_transform, symmetrize_operator export e_num, u_num, d_num, ud_num, half_ud_num, h_num export S_x, S_y, S_z, S_plus, S_min export u_plus_u_min, d_plus_d_min @@ -163,9 +162,13 @@ end # exact (Gaussian) integer entries promote to any scalar type without loss of precision const _PARTICLE_GAUGE = Complex{Int}[1 0 0 0; 0 -1 0 0; 0 0 im 0; 0 0 0 im] -# Symmetrize a Hubbard operator, inserting the staggered gauge for `SU2Irrep` particle symmetry -# If the staggered gauge commutes we don't incorporate it to retain the option for real tensors. -function _symmetrize_hubbard( +""" + symmetrize_operator(O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; tol=nothing) + +Symmetrize a Hubbard operator defined on `hubbard_space(Trivial, Trivial)` through the basis transformation for the requested particle and spin symmetries. +The staggered gauge required by `SU2Irrep` particle symmetry is applied automatically. +""" +function symmetrize_operator( O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; tol = nothing ) @@ -178,7 +181,7 @@ function _symmetrize_hubbard( W = reduce(⊗, Gs) Od = desymmetrize(O) - # check if commutes with the staggered gauge + # Do not incorporate the staggered gauge when it commutes, retaining real tensors. W * Od ≈ Od * W && return symmetrize(O, U, V; tol) scalartype(O) <: Real && @@ -186,27 +189,6 @@ function _symmetrize_hubbard( return symmetrize(O, map(g -> U * g, Gs), V; tol) end -# Symmetrize a Hubbard operator through its (staggered-gauge) basis transformation -_symmetrize_operator(O::AbstractTensorMap, P::Type{<:Sector}, S::Type{<:Sector}) = - _symmetrize_hubbard(O, P, S) - -""" - custom(A::AbstractArray, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; tol=nothing) - -Construct a symmetry-aware Hubbard operator from a dense rank-`2N` array in the reference basis ``|0⟩, |↑↓⟩, |↑⟩, |↓⟩``. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. -Fermion parity and the staggered gauge required by `SU2Irrep` particle symmetry are applied automatically. -""" -function custom( - A::AbstractArray, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}; tol = nothing - ) - O = _check_custom_space( - _custom_dense_operator(A), hubbard_space(Trivial, Trivial) - ) - return _symmetrize_hubbard(O, particle_symmetry, spin_symmetry; tol) -end - function n_site_operator(::Val{N}, elt::Type{<:Number}) where {N} V = hubbard_space(Trivial, Trivial) return zeros(elt, V^N ← V^N) diff --git a/src/quantumgates.jl b/src/quantumgates.jl index da82ae7..228527c 100644 --- a/src/quantumgates.jl +++ b/src/quantumgates.jl @@ -4,9 +4,8 @@ using TensorKit using LinearAlgebra: I using RationalRoots: signedroot import ..TensorKitTensors: symmetrize, desymmetrize, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export qubit_space, basis_transform, custom, controlled +export qubit_space, basis_transform, symmetrize_operator, controlled export pauli_x, pauli_y, pauli_z, proj_0, proj_1, hadamard, s_gate, t_gate export phase_shift, rotation_x, rotation_y, rotation_z, u3 export cnot, cy, cz, cphase, ch, cs, swap, iswap, dcx, ecr @@ -57,19 +56,15 @@ function basis_transform(symmetry::Type{<:Sector}) return TensorMap(Matrix{Int}(I, 2, 2), V ← qubit_space(Trivial)) end -# Symmetrize a gate through its basis transformation. The keywords of the parametrized gates -# are accepted and discarded: `qubit_space` and `basis_transform` do not depend on them. -_symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = - symmetrize(O, basis_transform(symmetry), qubit_space(symmetry)) - """ - custom(A::AbstractArray, symmetry::Type{<:Sector}; tol=nothing) + symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; tol=nothing, kwargs...) -Construct a symmetry-aware qubit operator from a dense rank-`2N` array in the computational basis. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and every axis must have dimension two. +Symmetrize a quantum gate defined on `qubit_space(Trivial)` through the basis transformation for `symmetry`. +Additional keywords used by parametrized gates are accepted but do not affect the transformation. """ -function custom(A::AbstractArray, symmetry::Type{<:Sector}; tol = nothing) - O = _check_custom_space(_custom_dense_operator(A), qubit_space(Trivial)) +function symmetrize_operator( + O::AbstractTensorMap, symmetry::Type{<:Sector}; tol = nothing, kwargs... + ) return symmetrize(O, basis_transform(symmetry), qubit_space(symmetry); tol) end diff --git a/src/spinoperators.jl b/src/spinoperators.jl index 7600186..17c6f75 100644 --- a/src/spinoperators.jl +++ b/src/spinoperators.jl @@ -4,9 +4,8 @@ using TensorKit using LinearAlgebra: I using RationalRoots: signedroot import ..TensorKitTensors: symmetrize, desymmetrize, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export spin_space, basis_transform, custom, casimir +export spin_space, basis_transform, symmetrize_operator, casimir export S_x, S_y, S_z, S_plus, S_min export S_x_S_x, S_y_S_y, S_z_S_z, S_plus_S_min, S_min_S_plus, S_exchange export σˣ, σʸ, σᶻ, σ⁺, σ⁻ @@ -86,21 +85,15 @@ function basis_transform(::Type{SU2Irrep}; spin = 1 // 2) return TensorMap(Matrix{Int}(I, d, d), desymmetrize(spin_space(SU2Irrep; spin)) ← V) end -# Symmetrize a spin operator through its basis transformation -_symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; kwargs...) = - symmetrize(O, basis_transform(symmetry; kwargs...), spin_space(symmetry; kwargs...)) - """ - custom(A::AbstractArray, symmetry::Type{<:Sector}; spin=1 // 2, tol=nothing) + symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; spin=1 // 2, tol=nothing) -Construct a symmetry-aware spin operator from a dense rank-`2N` array in the standard spin basis. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)` and each axis must have dimension `2spin + 1`. +Symmetrize a spin operator defined on `spin_space(Trivial; spin)` through the basis transformation for `symmetry`. """ -function custom( - A::AbstractArray, symmetry::Type{<:Sector}; spin = 1 // 2, tol = nothing +function symmetrize_operator( + O::AbstractTensorMap, symmetry::Type{<:Sector}; spin = 1 // 2, tol = nothing ) V = spin_space(symmetry; spin) - O = _check_custom_space(_custom_dense_operator(A), spin_space(Trivial; spin)) return symmetrize(O, basis_transform(symmetry; spin), V; tol) end diff --git a/src/tjoperators.jl b/src/tjoperators.jl index 4d4fc68..57a0bbb 100644 --- a/src/tjoperators.jl +++ b/src/tjoperators.jl @@ -4,9 +4,8 @@ using LinearAlgebra: diagind using TensorKit import ..HubbardOperators import ..TensorKitTensors: symmetrize, desymmetrize, fuse_local_operators, @operator -import ..TensorKitTensors: _custom_dense_operator, _check_custom_space -export tj_space, basis_transform, custom, tj_projector +export tj_space, basis_transform, symmetrize_operator, tj_projector export transform_slave_fermion # the operator names and their aliases are exported by the generation loop at the bottom of # this file, from the `_OPERATORS` registry @@ -154,35 +153,20 @@ function _state_indices(c, particle_symmetry::Type{<:Sector}, spin_symmetry::Typ end end -# Symmetrize a t-J operator through its basis transformation, in the requested basis. Note -# that `slave_fermion` is forwarded from the reference operator, which is already expressed -# in that basis (see `_maybe_slave_fermion`), such that only a permutation remains here. -function _symmetrize_operator( - O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, - spin_symmetry::Type{<:Sector}; kwargs... - ) - return symmetrize( - O, basis_transform(particle_symmetry, spin_symmetry; kwargs...), - tj_space(particle_symmetry, spin_symmetry; kwargs...) - ) -end - """ - custom(A::AbstractArray, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false, tol=nothing) + symmetrize_operator(O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false, tol=nothing) -Construct a symmetry-aware t-J operator from a dense rank-`2N` array in the normal reference basis ``|0⟩, |↑⟩, |↓⟩``. -The axes must be ordered as `(out₁, …, outₙ, in₁, …, inₙ)`. -Fermion parity is enforced, and `slave_fermion=true` applies the full slave-fermion transformation, before imposing the requested symmetries. +Symmetrize a t-J operator defined on `tj_space(Trivial, Trivial)` through the basis transformation for the requested particle and spin symmetries. +With `slave_fermion=true`, input in the normal t-J representation is transformed to the slave-fermion representation before imposing the requested symmetries. """ -function custom( - A::AbstractArray, particle_symmetry::Type{<:Sector}, +function symmetrize_operator( + O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false, tol = nothing ) - O = _check_custom_space(_custom_dense_operator(A), tj_space(Trivial, Trivial)) - if slave_fermion - V = tj_space(Trivial, Trivial) - U = basis_transform(Trivial, Trivial) - O = transform_slave_fermion(symmetrize(O, U, V; tol)) + # Public callers supply the normal reference representation, whereas `@operator` + # supplies an operator already transformed by `_maybe_slave_fermion`. + if slave_fermion && space(O) == _reference_operator_space(O) + O = transform_slave_fermion(O) end return symmetrize( O, basis_transform(particle_symmetry, spin_symmetry; slave_fermion), @@ -190,6 +174,13 @@ function custom( ) end +function _reference_operator_space(O::AbstractTensorMap) + N = numout(O) + N == numin(O) || return nothing + V = tj_space(Trivial, Trivial) + return V^N ← V^N +end + # Relation to the Hubbard model # ----------------------------- """ diff --git a/src/utils.jl b/src/utils.jl index f84a616..1be787d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -129,38 +129,6 @@ function _default_tol(::Type{T}, I::Type{<:Sector}) where {T <: Number} return sqrt(Tsector <: AbstractFloat ? max(ε, eps(Tsector)) : ε) end -# Convert the common dense-array representation of a square N-site operator to a TensorMap -# over ComplexSpace. The individual operator modules then check the inferred local dimension -# and project this map onto their own symmetric (and, where applicable, fermion-graded) spaces. -function _custom_dense_operator(A::AbstractArray) - eltype(A) <: Number || throw(ArgumentError("operator array must have a numeric element type")) - rank = ndims(A) - rank >= 2 && iseven(rank) || - throw(ArgumentError("operator array must have positive even rank `2N`, got rank $rank")) - local_dimension = size(A, 1) - expected_size = ntuple(Returns(local_dimension), rank) - size(A) == expected_size || throw( - ArgumentError( - "all operator-array legs must have the same size, got $(size(A))" - ) - ) - sites = rank ÷ 2 - V = ComplexSpace(local_dimension) - return TensorMap(A, V^sites ← V^sites) -end - -function _check_custom_space(O::AbstractTensorMap, V::ElementarySpace) - n = numin(O) - Vdense = desymmetrize(V) - expected_space = Vdense^n ← Vdense^n - space(O) == expected_space || throw( - ArgumentError( - "operator array has space `$(space(O))`, expected `$expected_space`" - ) - ) - return O -end - """ fuse_local_operators(O₁, O₂) @@ -250,7 +218,7 @@ function _operator_defs(alias, def, source) # symmetric terminal: # `op(elt::, s₁::Type{<:Sector}, …; kwargs...) = - # _symmetrize_operator(op(elt, Trivial, …; kwargs...), s₁, …; kwargs...)` + # symmetrize_operator(op(elt, Trivial, …; kwargs...), s₁, …; kwargs...)` tname = gensym(:elt) tsyms = [gensym(:S) for _ in 1:N] tsig = call_expr( @@ -258,7 +226,7 @@ function _operator_defs(alias, def, source) (Expr(:(::), tsyms[i], :(Type{<:Sector})) for i in 1:N)... ) refcall = call_expr(fname, tname, ntuple(Returns(:Trivial), N)...) - hookcall = call_expr(:_symmetrize_operator, refcall, tsyms...) + hookcall = call_expr(:symmetrize_operator, refcall, tsyms...) method_terminal = method_def(tsig, hookcall) (method_allsyms, method_elt, method_terminal) @@ -300,11 +268,11 @@ defaulting the unspecified symmetries of a multi-symmetry operator. If the reference method takes no symmetry arguments, it is itself the terminal: only the `op()` element-type default is generated (the symmetry-first, elt-only, and symmetric-terminal -methods, including the `_symmetrize_operator` hook, are omitted). +methods, including the `symmetrize_operator` hook, are omitted). Otherwise, the last generated method, the *symmetric terminal*, delegates to the module-local hook -`_symmetrize_operator(op(elt, Trivial, …; kwargs...), symmetry₁, …; kwargs...)`, which each +`symmetrize_operator(op(elt, Trivial, …; kwargs...), symmetry₁, …; kwargs...)`, which each operator module defines once (typically wrapping [`symmetrize`](@ref) with the module's `basis_transform` and local space). Every generated method simply accepts `kwargs...` and forwards them verbatim to both the reference call and the hook; the wrapped reference method diff --git a/test/custom.jl b/test/custom.jl deleted file mode 100644 index 2d613e0..0000000 --- a/test/custom.jl +++ /dev/null @@ -1,97 +0,0 @@ -using TensorKit -using Test -import TensorKitTensors -using TensorKitTensors: desymmetrize -import TensorKitTensors.SpinOperators as SO -import TensorKitTensors.BosonOperators as BO -import TensorKitTensors.FermionOperators as FO -import TensorKitTensors.HubbardOperators as HO -import TensorKitTensors.TJOperators as TJ -import TensorKitTensors.QuantumGates as QG - -@testset "SpinOperators" begin - Sz = convert(Array, SO.S_z()) - @test SO.custom(Sz, U1Irrep) ≈ SO.S_z(U1Irrep) - - Sz_spin1 = SO.S_z(Float64, Trivial; spin = 1) - @test SO.custom(convert(Array, Sz_spin1), U1Irrep; spin = 1) ≈ - SO.S_z(Float64, U1Irrep; spin = 1) - - # rank determines the number of sites and the scalar type is preserved - SS_big = SO.S_exchange(Complex{BigFloat}) - SS_custom = SO.custom(convert(Array, SS_big), SU2Irrep) - @test numout(SS_custom) == numin(SS_custom) == 2 - @test scalartype(SS_custom) === Complex{BigFloat} - @test SS_custom ≈ SO.S_exchange(Complex{BigFloat}, SU2Irrep) - - # tol is used only for symmetry projection - X_pert = copy(convert(Array, SO.S_x())) - X_pert[1, 1] += 1.0e-6 - @test_throws ArgumentError SO.custom(X_pert, Z2Irrep) - @test SO.custom(X_pert, Z2Irrep; tol = 1.0e-3) ≈ - SO.S_x(Z2Irrep) atol = 1.0e-5 - - # common dense-array validation - @test_throws ArgumentError SO.custom(fill("x", 2, 2), Trivial) - @test_throws ArgumentError SO.custom(zeros(2), Trivial) - @test_throws ArgumentError SO.custom(fill(1.0), Trivial) - @test_throws ArgumentError SO.custom(zeros(2, 3), Trivial) - @test_throws ArgumentError SO.custom(zeros(4, 4), Trivial) - @test_throws ArgumentError SO.custom(Sz, FermionParity) -end - -@testset "BosonOperators" begin - cutoff = 3 - N = BO.b_num(Float64, Trivial; cutoff) - @test BO.custom(convert(Array, N), U1Irrep; cutoff) ≈ - BO.b_num(Float64, U1Irrep; cutoff) - - hopping = BO.b_hopping(Float64, Trivial; cutoff) - @test BO.custom(convert(Array, hopping), U1Irrep; cutoff) ≈ - BO.b_hopping(Float64, U1Irrep; cutoff) -end - -@testset "FermionOperators" begin - hopping = FO.f_hopping(ComplexF64, Trivial) - A = convert(Array, desymmetrize(hopping)) - @test FO.custom(A, Trivial) ≈ hopping - @test FO.custom(A, U1Irrep) ≈ - FO.f_hopping(ComplexF64, U1Irrep) - - # `Trivial` still enforces the mandatory fermion-parity grading - parity_odd = ComplexF64[0 1; 0 0] - @test_throws ArgumentError FO.custom(parity_odd, Trivial) -end - -@testset "HubbardOperators" begin - hopping = HO.e_hopping(ComplexF64, Trivial, Trivial) - A = convert(Array, desymmetrize(hopping)) - @test HO.custom(A, SU2Irrep, SU2Irrep) ≈ - HO.e_hopping(ComplexF64, SU2Irrep, SU2Irrep) - - hopping_real = HO.e_hopping(Float64, Trivial, Trivial) - @test_throws ArgumentError HO.custom( - convert(Array, desymmetrize(hopping_real)), SU2Irrep, SU2Irrep - ) - - # test an operator with more than two-sites - paircor = HO.Δ⁺ij_Δkl(ComplexF64, Trivial, Trivial) - A = convert(Array, desymmetrize(paircor)) - @test HO.custom(A, U1Irrep, SU2Irrep) ≈ - HO.Δ⁺ij_Δkl(ComplexF64, U1Irrep, SU2Irrep) -end - -@testset "TJOperators" begin - hopping = TJ.e_hopping(ComplexF64, Trivial, Trivial) - A = convert(Array, desymmetrize(hopping)) - @test TJ.custom(A, U1Irrep, SU2Irrep; slave_fermion = true) ≈ - TJ.e_hopping(ComplexF64, U1Irrep, SU2Irrep; slave_fermion = true) -end - -@testset "QuantumGates" begin - Z = convert(Array, QG.pauli_z()) - @test QG.custom(Z, U1Irrep) ≈ QG.pauli_z(U1Irrep) - - CZ = convert(Array, QG.cz()) - @test QG.custom(CZ, U1Irrep) ≈ QG.cz(U1Irrep) -end diff --git a/test/symmetrize_operator.jl b/test/symmetrize_operator.jl new file mode 100644 index 0000000..76586a7 --- /dev/null +++ b/test/symmetrize_operator.jl @@ -0,0 +1,111 @@ +using TensorKit +using Test +import TensorKitTensors.SpinOperators as SO +import TensorKitTensors.BosonOperators as BO +import TensorKitTensors.FermionOperators as FO +import TensorKitTensors.HubbardOperators as HO +import TensorKitTensors.TJOperators as TJ +import TensorKitTensors.QuantumGates as QG + +function reference_operator(A::AbstractArray, V::ElementarySpace) + sites = ndims(A) ÷ 2 + return TensorMap(A, V^sites ← V^sites) +end + +@testset "exports" begin + for operators in (SO, BO, FO, HO, TJ, QG) + @test :symmetrize_operator in names(operators) + end +end + +@testset "SpinOperators" begin + Sz = convert(Array, SO.S_z()) + O = reference_operator(Sz, SO.spin_space(Trivial)) + @test SO.symmetrize_operator(O, U1Irrep) ≈ SO.S_z(U1Irrep) + + Sz_spin1 = SO.S_z(Float64, Trivial; spin = 1) + @test SO.symmetrize_operator(Sz_spin1, U1Irrep; spin = 1) ≈ + SO.S_z(Float64, U1Irrep; spin = 1) + + # the number of sites and scalar type are inherited from the reference TensorMap + SS_big = SO.S_exchange(Complex{BigFloat}) + SS_symmetric = SO.symmetrize_operator(SS_big, SU2Irrep) + @test numout(SS_symmetric) == numin(SS_symmetric) == 2 + @test scalartype(SS_symmetric) === Complex{BigFloat} + @test SS_symmetric ≈ SO.S_exchange(Complex{BigFloat}, SU2Irrep) + + # tol is used only for symmetry projection + X_pert = copy(convert(Array, SO.S_x())) + X_pert[1, 1] += 1.0e-6 + O_pert = reference_operator(X_pert, SO.spin_space(Trivial)) + @test_throws ArgumentError SO.symmetrize_operator(O_pert, Z2Irrep) + @test SO.symmetrize_operator(O_pert, Z2Irrep; tol = 1.0e-3) ≈ + SO.S_x(Z2Irrep) atol = 1.0e-5 + + @test_throws ArgumentError SO.symmetrize_operator(O, FermionParity) +end + +@testset "BosonOperators" begin + cutoff = 3 + N = BO.b_num(Float64, Trivial; cutoff) + @test BO.symmetrize_operator(N, U1Irrep; cutoff) ≈ + BO.b_num(Float64, U1Irrep; cutoff) + + hopping = BO.b_hopping(Float64, Trivial; cutoff) + @test BO.symmetrize_operator(hopping, U1Irrep; cutoff) ≈ + BO.b_hopping(Float64, U1Irrep; cutoff) +end + +@testset "FermionOperators" begin + hopping = FO.f_hopping(ComplexF64, Trivial) + @test FO.symmetrize_operator(hopping, Trivial) ≈ hopping + @test FO.symmetrize_operator(hopping, U1Irrep) ≈ + FO.f_hopping(ComplexF64, U1Irrep) + + # the reference TensorMap itself enforces the mandatory fermion-parity grading + parity_odd = ComplexF64[0 1; 0 0] + V = FO.fermion_space(Trivial) + @test_throws ArgumentError TensorMap(parity_odd, V ← V) +end + +@testset "HubbardOperators" begin + hopping = HO.e_hopping(ComplexF64, Trivial, Trivial) + @test HO.symmetrize_operator(hopping, SU2Irrep, SU2Irrep) ≈ + HO.e_hopping(ComplexF64, SU2Irrep, SU2Irrep) + + hopping_real = HO.e_hopping(Float64, Trivial, Trivial) + @test_throws ArgumentError HO.symmetrize_operator( + hopping_real, SU2Irrep, SU2Irrep + ) + + # test an operator with more than two sites + paircor = HO.Δ⁺ij_Δkl(ComplexF64, Trivial, Trivial) + @test HO.symmetrize_operator(paircor, U1Irrep, SU2Irrep) ≈ + HO.Δ⁺ij_Δkl(ComplexF64, U1Irrep, SU2Irrep) +end + +@testset "TJOperators" begin + hopping = TJ.e_hopping(ComplexF64, Trivial, Trivial) + hopping_slave = TJ.symmetrize_operator( + hopping, U1Irrep, SU2Irrep; slave_fermion = true + ) + generated_slave = TJ.e_hopping( + ComplexF64, U1Irrep, SU2Irrep; slave_fermion = true + ) + @test hopping_slave ≈ generated_slave + @test generated_slave ≈ TJ.transform_slave_fermion( + TJ.e_hopping(ComplexF64, U1Irrep, SU2Irrep) + ) +end + +@testset "QuantumGates" begin + Z = QG.pauli_z(Trivial) + @test QG.symmetrize_operator(Z, U1Irrep) ≈ QG.pauli_z(U1Irrep) + + CZ = QG.cz(Trivial) + @test QG.symmetrize_operator(CZ, U1Irrep) ≈ QG.cz(U1Irrep) + + theta = 0.37 + @test QG.rotation_z(U1Irrep; theta) ≈ + QG.symmetrize_operator(QG.rotation_z(Trivial; theta), U1Irrep) +end From dcc431033ac21ea017cadbf5fa83694407a1dcc6 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 15 Aug 2026 18:07:30 +0800 Subject: [PATCH 3/8] Minor docs update --- docs/src/index.md | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index e7db5db..ebc52f7 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -27,15 +27,11 @@ using TensorKit using TensorKitTensors.SpinOperators A = [1.0 0.0; 0.0 -1.0] V = spin_space(Trivial) -sites = ndims(A) ÷ 2 -O = TensorMap(A, V^sites ← V^sites) +O = TensorMap(A, V ← V) Z = symmetrize_operator(O, U1Irrep) ``` The caller is responsible for constructing a square homogeneous-site operator with the intended local space and array-axis ordering. -For fermionic systems, constructing the input on the module's `Trivial` space preserves the mandatory fermion-parity grading. -Hubbard and t-J operators take separate particle and spin symmetries. -For t-J, `slave_fermion=true` transforms an operator supplied in the normal t-J representation into the slave-fermion representation. ```@docs symmetrize From c48b33cb6e2c5f456c4d24ab6dfa2b9906bc60a5 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 15 Aug 2026 22:34:47 +0800 Subject: [PATCH 4/8] Drop dedicated tests for symmetrize_operator --- test/symmetrize_operator.jl | 111 ------------------------------------ 1 file changed, 111 deletions(-) delete mode 100644 test/symmetrize_operator.jl diff --git a/test/symmetrize_operator.jl b/test/symmetrize_operator.jl deleted file mode 100644 index 76586a7..0000000 --- a/test/symmetrize_operator.jl +++ /dev/null @@ -1,111 +0,0 @@ -using TensorKit -using Test -import TensorKitTensors.SpinOperators as SO -import TensorKitTensors.BosonOperators as BO -import TensorKitTensors.FermionOperators as FO -import TensorKitTensors.HubbardOperators as HO -import TensorKitTensors.TJOperators as TJ -import TensorKitTensors.QuantumGates as QG - -function reference_operator(A::AbstractArray, V::ElementarySpace) - sites = ndims(A) ÷ 2 - return TensorMap(A, V^sites ← V^sites) -end - -@testset "exports" begin - for operators in (SO, BO, FO, HO, TJ, QG) - @test :symmetrize_operator in names(operators) - end -end - -@testset "SpinOperators" begin - Sz = convert(Array, SO.S_z()) - O = reference_operator(Sz, SO.spin_space(Trivial)) - @test SO.symmetrize_operator(O, U1Irrep) ≈ SO.S_z(U1Irrep) - - Sz_spin1 = SO.S_z(Float64, Trivial; spin = 1) - @test SO.symmetrize_operator(Sz_spin1, U1Irrep; spin = 1) ≈ - SO.S_z(Float64, U1Irrep; spin = 1) - - # the number of sites and scalar type are inherited from the reference TensorMap - SS_big = SO.S_exchange(Complex{BigFloat}) - SS_symmetric = SO.symmetrize_operator(SS_big, SU2Irrep) - @test numout(SS_symmetric) == numin(SS_symmetric) == 2 - @test scalartype(SS_symmetric) === Complex{BigFloat} - @test SS_symmetric ≈ SO.S_exchange(Complex{BigFloat}, SU2Irrep) - - # tol is used only for symmetry projection - X_pert = copy(convert(Array, SO.S_x())) - X_pert[1, 1] += 1.0e-6 - O_pert = reference_operator(X_pert, SO.spin_space(Trivial)) - @test_throws ArgumentError SO.symmetrize_operator(O_pert, Z2Irrep) - @test SO.symmetrize_operator(O_pert, Z2Irrep; tol = 1.0e-3) ≈ - SO.S_x(Z2Irrep) atol = 1.0e-5 - - @test_throws ArgumentError SO.symmetrize_operator(O, FermionParity) -end - -@testset "BosonOperators" begin - cutoff = 3 - N = BO.b_num(Float64, Trivial; cutoff) - @test BO.symmetrize_operator(N, U1Irrep; cutoff) ≈ - BO.b_num(Float64, U1Irrep; cutoff) - - hopping = BO.b_hopping(Float64, Trivial; cutoff) - @test BO.symmetrize_operator(hopping, U1Irrep; cutoff) ≈ - BO.b_hopping(Float64, U1Irrep; cutoff) -end - -@testset "FermionOperators" begin - hopping = FO.f_hopping(ComplexF64, Trivial) - @test FO.symmetrize_operator(hopping, Trivial) ≈ hopping - @test FO.symmetrize_operator(hopping, U1Irrep) ≈ - FO.f_hopping(ComplexF64, U1Irrep) - - # the reference TensorMap itself enforces the mandatory fermion-parity grading - parity_odd = ComplexF64[0 1; 0 0] - V = FO.fermion_space(Trivial) - @test_throws ArgumentError TensorMap(parity_odd, V ← V) -end - -@testset "HubbardOperators" begin - hopping = HO.e_hopping(ComplexF64, Trivial, Trivial) - @test HO.symmetrize_operator(hopping, SU2Irrep, SU2Irrep) ≈ - HO.e_hopping(ComplexF64, SU2Irrep, SU2Irrep) - - hopping_real = HO.e_hopping(Float64, Trivial, Trivial) - @test_throws ArgumentError HO.symmetrize_operator( - hopping_real, SU2Irrep, SU2Irrep - ) - - # test an operator with more than two sites - paircor = HO.Δ⁺ij_Δkl(ComplexF64, Trivial, Trivial) - @test HO.symmetrize_operator(paircor, U1Irrep, SU2Irrep) ≈ - HO.Δ⁺ij_Δkl(ComplexF64, U1Irrep, SU2Irrep) -end - -@testset "TJOperators" begin - hopping = TJ.e_hopping(ComplexF64, Trivial, Trivial) - hopping_slave = TJ.symmetrize_operator( - hopping, U1Irrep, SU2Irrep; slave_fermion = true - ) - generated_slave = TJ.e_hopping( - ComplexF64, U1Irrep, SU2Irrep; slave_fermion = true - ) - @test hopping_slave ≈ generated_slave - @test generated_slave ≈ TJ.transform_slave_fermion( - TJ.e_hopping(ComplexF64, U1Irrep, SU2Irrep) - ) -end - -@testset "QuantumGates" begin - Z = QG.pauli_z(Trivial) - @test QG.symmetrize_operator(Z, U1Irrep) ≈ QG.pauli_z(U1Irrep) - - CZ = QG.cz(Trivial) - @test QG.symmetrize_operator(CZ, U1Irrep) ≈ QG.cz(U1Irrep) - - theta = 0.37 - @test QG.rotation_z(U1Irrep; theta) ≈ - QG.symmetrize_operator(QG.rotation_z(Trivial; theta), U1Irrep) -end From 4ea2605bf5e553116d28210a879cc6797355387a Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 15 Aug 2026 22:39:54 +0800 Subject: [PATCH 5/8] Docs update again --- docs/src/index.md | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index ebc52f7..ddd0608 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -20,19 +20,16 @@ Operators that are incompatible with a given symmetry throw an `ArgumentError`. ### Custom operators -Construct a `TensorMap` on the relevant module's `Trivial` reference space, then pass it to that module's `symmetrize_operator` with the desired symmetry types: +For each operator module, one can define a custom operator in `TensorMap` format in the non-symmetric reference basis, and then impose the symmetries with that module's `symmetrize_operator` function: ```julia using TensorKit using TensorKitTensors.SpinOperators -A = [1.0 0.0; 0.0 -1.0] V = spin_space(Trivial) -O = TensorMap(A, V ← V) +O = TensorMap([1.0 0.0; 0.0 -1.0], V ← V) Z = symmetrize_operator(O, U1Irrep) ``` -The caller is responsible for constructing a square homogeneous-site operator with the intended local space and array-axis ordering. - ```@docs symmetrize desymmetrize From 0a77820d7a0f2bf22f0f86b2c5539d885addfb05 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sat, 15 Aug 2026 22:40:01 +0800 Subject: [PATCH 6/8] Bump patch version --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6a768aa..37d11e4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "TensorKitTensors" uuid = "41b62e7d-e9d1-4e23-942c-79a97adf954b" -version = "0.3.0" +version = "0.3.1" authors = ["QuantumKitHub"] [workspace] From 75a32bb68e4feecd16767e75fe1424b9142cd9f9 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Sun, 16 Aug 2026 10:40:34 +0800 Subject: [PATCH 7/8] Streamline t-J construction in slave-fermion basis --- docs/src/index.md | 2 ++ src/tjoperators.jl | 42 +++++++++++------------------------------- 2 files changed, 13 insertions(+), 31 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index ddd0608..c2b4a38 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -30,6 +30,8 @@ O = TensorMap([1.0 0.0; 0.0 -1.0], V ← V) Z = symmetrize_operator(O, U1Irrep) ``` +For `TJOperators`, construct the input on `tj_space(Trivial, Trivial; slave_fermion)` in the representation selected by the same `slave_fermion` keyword passed to `symmetrize_operator`. + ```@docs symmetrize desymmetrize diff --git a/src/tjoperators.jl b/src/tjoperators.jl index 57a0bbb..96da519 100644 --- a/src/tjoperators.jl +++ b/src/tjoperators.jl @@ -156,31 +156,19 @@ end """ symmetrize_operator(O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool=false, tol=nothing) -Symmetrize a t-J operator defined on `tj_space(Trivial, Trivial)` through the basis transformation for the requested particle and spin symmetries. -With `slave_fermion=true`, input in the normal t-J representation is transformed to the slave-fermion representation before imposing the requested symmetries. +Symmetrize a t-J operator defined on `tj_space(Trivial, Trivial; slave_fermion)` through the basis transformation for the requested particle and spin symmetries. +The input must already use the normal or slave-fermion representation selected by `slave_fermion`. """ function symmetrize_operator( O::AbstractTensorMap, particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:Sector}; slave_fermion::Bool = false, tol = nothing ) - # Public callers supply the normal reference representation, whereas `@operator` - # supplies an operator already transformed by `_maybe_slave_fermion`. - if slave_fermion && space(O) == _reference_operator_space(O) - O = transform_slave_fermion(O) - end return symmetrize( O, basis_transform(particle_symmetry, spin_symmetry; slave_fermion), tj_space(particle_symmetry, spin_symmetry; slave_fermion); tol ) end -function _reference_operator_space(O::AbstractTensorMap) - N = numout(O) - N == numin(O) || return nothing - V = tj_space(Trivial, Trivial) - return V^N ← V^N -end - # Relation to the Hubbard model # ----------------------------- """ @@ -221,16 +209,13 @@ function tj_projector(particle_symmetry::Type{<:Sector}, spin_symmetry::Type{<:S return proj end -# Project a Hubbard reference operator onto the t-J space, and express it in the requested -# basis. The projector is parity-even, so it distributes over `⊗` and projecting the reference -# operator is equivalent to projecting the symmetric one. Since it has integer entries, the -# scalar type of the Hubbard operator is preserved. -function _project_hubbard(O::AbstractTensorMap, slave_fermion::Bool) - Pⁿ = _reference_projector(Val(numout(O))) - return _maybe_slave_fermion(Pⁿ * O * Pⁿ', slave_fermion) -end -function _reference_projector(::Val{N}) where {N} - return reduce(⊗, ntuple(Returns(tj_projector(Trivial, Trivial)), Val(N))) +"Project a Hubbard reference operator onto the normal t-J space." +function _project_hubbard(O::AbstractTensorMap) + # The projector is parity-even, so it distributes over `⊗` + n = numout(O) + Pⁿ = reduce(⊗, ntuple(Returns(tj_projector(Trivial, Trivial)), Val(n))) + # integer entries in `Pⁿ` preserve the scalar type of `O` + return Pⁿ * O * Pⁿ' end # Slave-fermion basis @@ -271,12 +256,6 @@ function transform_slave_fermion(V::ElementarySpace) return fuse(V, V_aux) end -# Express a reference operator in the requested basis. Reference operators are always obtained -# in the plain t-J basis, and transformed only once, at the very end: the slave-fermion -# transformation does not distribute over `⊗`. -_maybe_slave_fermion(O::AbstractTensorMap, slave_fermion::Bool) = - slave_fermion ? transform_slave_fermion(O) : O - # Operators # --------- @@ -360,7 +339,8 @@ for (name, alias) in [ elt::Type{<:Number}, ::Type{Trivial}, ::Type{Trivial}; slave_fermion::Bool = false ) - return _project_hubbard(HubbardOperators.$name(elt, Trivial, Trivial), slave_fermion) + O = _project_hubbard(HubbardOperators.$name(elt, Trivial, Trivial)) + return slave_fermion ? transform_slave_fermion(O) : O end end From 13ab67e1bea68aea1f7a699fc53d20f74e6d8a56 Mon Sep 17 00:00:00 2001 From: Yue Zhengyuan Date: Mon, 17 Aug 2026 20:43:04 +0800 Subject: [PATCH 8/8] Removed instructions for constructing input on TJOperators from index.md [skip ci]. --- docs/src/index.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/docs/src/index.md b/docs/src/index.md index c2b4a38..ddd0608 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -30,8 +30,6 @@ O = TensorMap([1.0 0.0; 0.0 -1.0], V ← V) Z = symmetrize_operator(O, U1Irrep) ``` -For `TJOperators`, construct the input on `tj_space(Trivial, Trivial; slave_fermion)` in the representation selected by the same `slave_fermion` keyword passed to `symmetrize_operator`. - ```@docs symmetrize desymmetrize