Computing the expectation value of J^2 fails for non-relativistic spin-orbitals

[jagot]

Setup:

julia> using AtomicLevels

julia> using EnergyExpressions

julia> import EnergyExpressions: OrbitalMatrixElement

julia> using AngularMomentumAlgebra

julia> import AngularMomentumAlgebra: integrate_spinor

julia> 𝐉 = TotalAngularMomentum()
𝐉̂⁽¹⁾

julia> 𝐉² = 𝐉⋅𝐉
(𝐉̂⁽¹⁾⋅𝐉̂⁽¹⁾)

The one-body part $j_i^2$ of $J^2 = j_i^2 + \mathbf{j}_i\cdot\mathbf{j}_j$ works fine:

julia> integrate_spinor(OrbitalMatrixElement([so"ls₀α"], TensorOperator{1}(𝐉²), [so"1s₀α"]))
0.75⟨ls₀α|1s₀α⟩ᵣ

whereas the two-body part $\mathbf{j}_i\cdot\mathbf{j}_j$ does not:

julia> integrate_spinor(OrbitalMatrixElement([so"1s₀β", so"ls₀α"], TensorOperator{2}(𝐉²), [so"1s₀α", so"ks₀β"]))
ERROR: MethodError: no method matching -(::Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}}, ::Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}})
Closest candidates are:
  -(::LC, ::T) where {T, N, LC<:(AngularMomentumAlgebra.LinearCombination{<:T, N})} at ~/.julia/dev/AngularMomentumAlgebra/src/linear_combinations.jl:61
  -(::T, ::LC) where {T, N, LC<:(AngularMomentumAlgebra.LinearCombination{<:T, N})} at ~/.julia/dev/AngularMomentumAlgebra/src/linear_combinations.jl:64
Stacktrace:
 [1] matrix_element2(γjm₁′::Tuple{Tuple{Int64, Int64}, Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}}}, γjm₂′::Tuple{Tuple{Int64, Int64}, Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}}}, X::TensorScalarProduct{TotalAngularMomentum, TotalAngularMomentum}, γjm₁::Tuple{Tuple{Int64, Int64}, Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}}}, γjm₂::Tuple{Tuple{Int64, Int64}, Tuple{HalfIntegers.Half{Int64}, HalfIntegers.Half{Int64}}})
   @ AngularMomentumAlgebra ~/.julia/dev/AngularMomentumAlgebra/src/tensor_matrix_elements.jl:451
 [2] matrix_element(::Tuple{TotalAngularMomentumSubSystem, TotalAngularMomentumSubSystem}, ::Tuple{SpinOrbital{Orbital{Int64}, Tuple{Int64, HalfIntegers.Half{Int64}}}, SpinOrbital{Orbital{Symbol}, Tuple{Int64, HalfIntegers.Half{Int64}}}}, X::TensorScalarProduct{TotalAngularMomentum, TotalAngularMomentum}, ::Tuple{SpinOrbital{Orbital{Int64}, Tuple{Int64, HalfIntegers.Half{Int64}}}, SpinOrbital{Orbital{Symbol}, Tuple{Int64, HalfIntegers.Half{Int64}}}})
   @ AngularMomentumAlgebra ~/.julia/dev/AngularMomentumAlgebra/src/tensor_matrix_elements.jl:979
 [3] dot(::Tuple{SpinOrbital{Orbital{Int64}, Tuple{Int64, HalfIntegers.Half{Int64}}}, SpinOrbital{Orbital{Symbol}, Tuple{Int64, HalfIntegers.Half{Int64}}}}, X::TensorScalarProduct{TotalAngularMomentum, TotalAngularMomentum}, ::Tuple{SpinOrbital{Orbital{Int64}, Tuple{Int64, HalfIntegers.Half{Int64}}}, SpinOrbital{Orbital{Symbol}, Tuple{Int64, HalfIntegers.Half{Int64}}}})
   @ AngularMomentumAlgebra ~/.julia/dev/AngularMomentumAlgebra/src/tensor_matrix_elements.jl:687
 [4] spin_ang_coeff(me::OrbitalMatrixElement{2, SpinOrbital{<:Orbital, Tuple{Int64, HalfIntegers.Half{Int64}}}, TensorOperator{2, TensorScalarProduct{TotalAngularMomentum, TotalAngularMomentum}}, SpinOrbital{<:Orbital, Tuple{Int64, HalfIntegers.Half{Int64}}}})
   @ AngularMomentumAlgebra ~/.julia/dev/AngularMomentumAlgebra/src/tensors.jl:304
 [5] integrate_spinor(me::OrbitalMatrixElement{2, SpinOrbital{<:Orbital, Tuple{Int64, HalfIntegers.Half{Int64}}}, TensorOperator{2, TensorScalarProduct{TotalAngularMomentum, TotalAngularMomentum}}, SpinOrbital{<:Orbital, Tuple{Int64, HalfIntegers.Half{Int64}}}})
   @ AngularMomentumAlgebra ~/.julia/dev/AngularMomentumAlgebra/src/tensors.jl:311
 [6] top-level scope
   @ REPL[11]:1