Delete spin_operators_pair

Export instead `to_product_space`

src/Operators/TensorOperators.jl

@@ -78,13 +78,8 @@ function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
     return ret
 end
 
-"""
-    to_product_space(A, B)
 
-Given two lists of local operators acting on sites i and j individually, return
-the corresponding local operators that act on the tensor product space.
-"""
-function to_product_space(A, B)
+function to_product_space_orig(A, B)
     (isempty(A) || isempty(B)) && error("Nonempty lists required")
 
     @assert allequal(size.(A))
@@ -97,3 +92,29 @@ function to_product_space(A, B)
     I2 = Ref(Matrix(I, N2, N2))
     return (kron.(A, I2), kron.(I1, B))
 end
+
+"""
+    to_product_space(A, B, Cs...)
+
+Given lists of operators acting on local Hilbert spaces individually, return the
+corresponding operators that act on the tensor product space. In typical usage,
+the inputs will represent local physical observables and the outputs will be
+used to define quantum couplings.
+"""
+function to_product_space(A, B, Cs...)
+    lists = [A, B, Cs...]
+
+    Ns = map(enumerate(lists)) do (i, list)
+        isempty(list) && error("Empty operator list in argument $i.")
+        allequal(size.(list)) || error("Unequal sized operators in argument $i.")
+        return size(first(list), 1)
+    end
+
+    return map(enumerate(lists)) do (i, list)
+        I1 = I(prod(Ns[begin:i-1]))
+        I2 = I(prod(Ns[i+1:end]))
+        return map(list) do op
+            kron(I1, op, I2)
+        end
+    end
+end

src/Sunny.jl

@@ -39,7 +39,7 @@ include("Operators/Stevens.jl")
 include("Operators/TensorOperators.jl")
 include("Operators/Symbolic.jl")
 include("Operators/Observables.jl")
-export spin_matrices, rotate_operator, print_stevens_expansion
+export spin_matrices, to_product_space, rotate_operator, print_stevens_expansion
 
 include("Symmetry/LatticeUtils.jl")
 include("Symmetry/SymOp.jl")
@@ -65,8 +65,8 @@ include("System/Ewald.jl")
 include("System/Interactions.jl")
 export SpinInfo, System, Site, eachsite, position_to_site, global_position, magnetic_moment, 
     set_coherent!, set_dipole!, polarize_spins!, randomize_spins!, energy, energy_per_site,
-    spin_operators, stevens_operators, large_S_spin_operators, large_S_stevens_operators, set_onsite_coupling!, 
-    spin_operators_pair, set_pair_coupling!, set_external_field!, set_exchange!, dmvec, enable_dipole_dipole!,
+    spin_operators, stevens_operators, large_S_spin_operators, large_S_stevens_operators,
+    set_onsite_coupling!, set_pair_coupling!, set_external_field!, set_exchange!, dmvec, enable_dipole_dipole!,
     to_inhomogeneous, set_external_field_at!, set_vacancy_at!, set_onsite_coupling_at!,
     symmetry_equivalent_bonds, set_exchange_at!, remove_periodicity!
 

src/System/PairExchange.jl

@@ -167,34 +167,22 @@ function set_pair_coupling_aux!(sys::System, bilin::Union{Float64, Mat3}, biquad
 end
 
 
-"""
-    spin_operators_pair(sys::System, i::Int, j::Int)
-
-Returns a pair of spin dipoles `(Si, Sj)` associated with atoms `i` and `j`,
-respectively. The components of these return values are operators that act in
-the tensor product space of the two sites, which makes them useful for defining
-interactions as input to [`set_pair_coupling!`](@ref).
-"""
-function spin_operators_pair(sys::System{N}, i::Int, j::Int) where N
-    Si = spin_matrices(N=sys.Ns[i])
-    Sj = spin_matrices(N=sys.Ns[j])
-    return to_product_space(Si, Sj)
-end
-
 """
     set_pair_coupling!(sys::System, coupling, bond)
 
 Sets an arbitrary `coupling` along `bond`. This coupling will be propagated to
 equivalent bonds in consistency with crystal symmetry. Any previous interactions
 on these bonds will be overwritten. The parameter `bond` has the form `Bond(i,
 j, offset)`, where `i` and `j` are atom indices within the unit cell, and
-`offset` is a displacement in unit cells. The `coupling` may be formed as a
-polynomial of operators obtained from [`spin_operators_pair`](@ref).
+`offset` is a displacement in unit cells. The `coupling` is a represented as a
+matrix acting in the tensor product space of the two sites, and typically
+originates from [`to_product_space`](@ref).
 
 # Examples
 ```julia
 # Add a bilinear and biquadratic exchange
-S = spin_operators_pair(sys, bond.i, bond.j)
+S = spin_matrices(1/2)
+Si, Sj = to_product_space(S, S)
 set_pair_coupling!(sys, Si'*J1*Sj + (Si'*J2*Sj)^2, bond)
 ```
 """

test/shared.jl

@@ -44,8 +44,11 @@ function add_quadratic_interactions!(sys, mode)
                   Γ   J   -D;
                   D   D  J+K]
 
-        Si, Sj = spin_operators_pair(sys, 1, 2)
-        set_pair_coupling!(sys, Si'*J_exch*Sj + 0.01(Si'*Sj)^2, Bond(1, 2, [0, 0, 0]))
+        bond = Bond(1, 2, [0, 0, 0])
+        S1 = spin_matrices(; N=sys.Ns[bond.i])
+        S2 = spin_matrices(; N=sys.Ns[bond.j])
+        Si, Sj = to_product_space(S1, S2)
+        set_pair_coupling!(sys, Si'*J_exch*Sj + 0.01(Si'*Sj)^2, bond)
         =#
     end
 end

test/test_tensors.jl

@@ -48,11 +48,12 @@ end
     E = energy(sys)
     dE_dZ = Sunny.energy_grad_coherents(sys)
 
-    Si, Sj = spin_operators_pair(sys, bond.i, bond.j)
+    S = spin_matrices(; N=5)
+    Si, Sj = to_product_space(S, S)
     set_pair_coupling!(sys, Si'*J_exch*Sj, bond; fast=false)
     E′ = energy(sys)
     dE_dZ′ = Sunny.energy_grad_coherents(sys)
 
     @test E ≈ E′
-    @test dE_dZ ≈ dE_dZ′        
+    @test dE_dZ ≈ dE_dZ′
 end