Export symmetrize_operator for custom operator symmetrization

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

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
 # ---------------------
 """

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})

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,29 +167,46 @@ 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))
     W = reduce(⊗, Gs)
     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)

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
 # ----------------
 """

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)

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
 # -----------------------------
 """

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₂)
 

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