Export symmetrize_operator for custom operator symmetrization

docs/src/index.md

@@ -14,13 +14,24 @@ 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
+
+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]
+V = spin_space(Trivial)
+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.
 
 ```@docs
 symmetrize

src/bosonoperators.jl

@@ -4,7 +4,7 @@ using TensorKit
 using LinearAlgebra: I
 import ..TensorKitTensors: symmetrize, desymmetrize, @operator
 
-export boson_space, basis_transform
+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
@@ -58,9 +58,19 @@ 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...))
+"""
+    symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; cutoff::Integer, tol=nothing)
+
+Symmetrize a bosonic operator defined on `boson_space(Trivial; cutoff)` through the basis transformation for `symmetry`.
+"""
+function symmetrize_operator(
+        O::AbstractTensorMap, symmetry::Type{<:Sector}; cutoff::Integer,
+        tol = nothing
+    )
+    return symmetrize(
+        O, basis_transform(symmetry; cutoff), boson_space(symmetry; cutoff); tol
+    )
+end
 
 # Single-site operators
 # ---------------------

src/fermionoperators.jl

@@ -4,7 +4,7 @@ using TensorKit
 using LinearAlgebra: I
 import ..TensorKitTensors: symmetrize, desymmetrize, @operator
 
-export fermion_space, basis_transform
+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
@@ -55,9 +55,17 @@ 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))
+"""
+    symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; tol=nothing)
+
+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 symmetrize_operator(
+        O::AbstractTensorMap, symmetry::Type{<:Sector}; tol = nothing
+    )
+    return symmetrize(O, basis_transform(symmetry), fermion_space(symmetry); tol)
+end
 
 # Single-site operators
 # ---------------------

src/hubbardoperators.jl

@@ -3,7 +3,7 @@ module HubbardOperators
 using TensorKit
 import ..TensorKitTensors: symmetrize, desymmetrize, @operator
 
-export hubbard_space, basis_transform
+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
@@ -162,33 +162,33 @@ 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}
+        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)
+    # 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 &&
         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)
-
 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

@@ -5,7 +5,7 @@ using LinearAlgebra: I
 using RationalRoots: signedroot
 import ..TensorKitTensors: symmetrize, desymmetrize, @operator
 
-export qubit_space, basis_transform, 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
@@ -56,10 +56,17 @@ 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))
+"""
+    symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; tol=nothing, kwargs...)
+
+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 symmetrize_operator(
+        O::AbstractTensorMap, symmetry::Type{<:Sector}; tol = nothing, kwargs...
+    )
+    return symmetrize(O, basis_transform(symmetry), qubit_space(symmetry); tol)
+end
 
 # Gate combinators
 # ----------------

src/spinoperators.jl

@@ -5,7 +5,7 @@ using LinearAlgebra: I
 using RationalRoots: signedroot
 import ..TensorKitTensors: symmetrize, desymmetrize, @operator
 
-export spin_space, basis_transform, 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 σˣ, σʸ, σᶻ, σ⁺, σ⁻
@@ -85,9 +85,17 @@ 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...))
+"""
+    symmetrize_operator(O::AbstractTensorMap, symmetry::Type{<:Sector}; spin=1 // 2, tol=nothing)
+
+Symmetrize a spin operator defined on `spin_space(Trivial; spin)` through the basis transformation for `symmetry`.
+"""
+function symmetrize_operator(
+        O::AbstractTensorMap, symmetry::Type{<:Sector}; spin = 1 // 2, tol = nothing
+    )
+    V = spin_space(symmetry; spin)
+    return symmetrize(O, basis_transform(symmetry; spin), V; tol)
+end
 
 # Pauli matrices
 # --------------

src/tjoperators.jl

@@ -5,7 +5,7 @@ using TensorKit
 import ..HubbardOperators
 import ..TensorKitTensors: symmetrize, desymmetrize, fuse_local_operators, @operator
 
-export tj_space, basis_transform, 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
@@ -153,19 +153,34 @@ 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(
+"""
+    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.
+"""
+function symmetrize_operator(
         O::AbstractTensorMap, particle_symmetry::Type{<:Sector},
-        spin_symmetry::Type{<:Sector}; kwargs...
+        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; kwargs...),
-        tj_space(particle_symmetry, spin_symmetry; kwargs...)
+        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
 # -----------------------------
 """

src/utils.jl

@@ -218,15 +218,15 @@ function _operator_defs(alias, def, source)
 
         # symmetric terminal:
         # `op(elt::<c>, 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(
             fname, Expr(:(::), tname, eltconstraint),
             (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)
@@ -268,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