Progress

src/Reshaping.jl

@@ -46,7 +46,7 @@ function set_interactions_from_origin!(sys::System{N}) where N
         for pc in orig_ints[i].pair
             new_bond = transform_bond(sys.crystal, new_i, origin.crystal, pc.bond)
             isculled = bond_parity(new_bond)
-            push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.biquad, pc.general))
+            push!(new_pc, PairCoupling(isculled, new_bond, pc.scalar, pc.bilin, pc.scalar_biquad, pc.general))
         end
         new_pc = sort!(new_pc, by=c->c.isculled)
         new_ints[new_i].pair = new_pc

src/SpinWaveTheory/SWTCalculations.jl

@@ -146,7 +146,7 @@ function swt_hamiltonian_SUN!(swt::SpinWaveTheory, q_reshaped::Vec3, Hmat::Matri
 
             ### Biquadratic exchange
 
-            J = coupling.biquad
+            J = coupling.scalar_biquad
 
 
             Ti_11 = view(sun_basis_i, 1, 1, :)
@@ -313,7 +313,7 @@ function swt_hamiltonian_dipole!(swt::SpinWaveTheory, q_reshaped::Vec3, Hmat::Ma
 
             ### Biquadratic exchange
 
-            J = coupling.biquad
+            J = coupling.scalar_biquad
             # ⟨Ω₂, Ω₁|(𝐒₁⋅𝐒₂)^2|Ω₁, Ω₂⟩ = (1-1/S+1/(4S^2)) (Ω₁⋅Ω₂)^2 - 1/2 Ω₁⋅Ω₂ + const.
             # The biquadratic part including the biquadratic scaling factor.
             Ri = R_mat[sub_i]

src/System/Interactions.jl

@@ -161,8 +161,8 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N
         ΔE += dot(Δs, J, sⱼ)
 
         # Biquadratic
-        if !iszero(pc.biquad)
-            J = pc.biquad
+        if !iszero(pc.scalar_biquad)
+            J = pc.scalar_biquad
             if sys.mode in (:dipole, :dipole_large_S)
                 ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2)
             else
@@ -272,8 +272,8 @@ function energy_aux(ints::Interactions, sys::System{N}, i::Int, cells) where N
             E += dot(sᵢ, J, sⱼ)
 
             # Biquadratic
-            if !iszero(pc.biquad)
-                J = pc.biquad
+            if !iszero(pc.scalar_biquad)
+                J = pc.scalar_biquad
                 if sys.mode in (:dipole, :dipole_large_S)
                     E += J * (sᵢ⋅sⱼ)^2
                 else
@@ -359,8 +359,8 @@ function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Inter
             ∇E[cellⱼ, bond.j] += J' * sᵢ
 
             # Biquadratic
-            if !iszero(pc.biquad)
-                J = pc.biquad
+            if !iszero(pc.scalar_biquad)
+                J = pc.scalar_biquad
                 if sys.mode in (:dipole, :dipole_large_S)
                     ∇E[cellᵢ, bond.i] += J * 2sⱼ*(sᵢ⋅sⱼ)
                     ∇E[cellⱼ, bond.j] += J * 2sᵢ*(sᵢ⋅sⱼ)

src/System/PairExchange.jl

@@ -29,37 +29,18 @@ function to_float_or_mat3(J)
     return J::Union{Float64, Mat3}
 end
 
-# Perform renormalization derived in https://arxiv.org/abs/2304.03874
-function renormalize_biquad(sys, bond, bilin, biquad)
-    @assert sys.mode == :dipole
-
-    # Multiplicative average is the correct result for two sites with
-    # different S (derivation omitted).
-    S1 = (sys.Ns[bond.i]-1)/2
-    S2 = (sys.Ns[bond.j]-1)/2
-    S = sqrt(S1*S2)
-    r = (1 - 1/S + 1/4S^2)
-
-    # biquad (s⋅sⱼ)^2 -> biquad (r (sᵢ⋅sⱼ)^2 - (sᵢ⋅sⱼ)/2 + S^3 + S^2/4)
-    bilin = bilin - (biquad/2) * one(bilin)
-    biquad = r * biquad
-
-    # Drop the constant shift, `biquad (S^3 + S^2/4)`.
-
-    return (bilin, biquad)
-end
 
 # Internal function only
-function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Float64, tensordec::TensorDecomposition)
+function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Float64, Mat3}, scalar_biquad::Float64, tensordec::TensorDecomposition)
     # Remove previous coupling on this bond
     filter!(c -> c.bond != bond, couplings)
 
     # If the new coupling is exactly zero, return early
-    iszero(bilin) && iszero(biquad) && isempty(tensordec.data) && return
+    iszero(bilin) && iszero(scalar_biquad) && isempty(tensordec.data) && return
 
     # Otherwise, add the new coupling to the list
     isculled = bond_parity(bond)
-    push!(couplings, PairCoupling(isculled, bond, scalar, bilin, biquad, tensordec))
+    push!(couplings, PairCoupling(isculled, bond, scalar, bilin, scalar_biquad, tensordec))
 
     # Sorting after each insertion will introduce quadratic scaling in length of
     # `couplings`. In typical usage, the `couplings` list will be short.
@@ -69,16 +50,17 @@ function push_coupling!(couplings, bond::Bond, scalar::Float64, bilin::Union{Flo
 end
 
 
-function decompose_general_coupling(op, gen1, gen2; fast)
-    N1 = size(gen1[1], 1)
-    N2 = size(gen2[1], 1)
+function decompose_general_coupling(op, N1, N2; fast)
+    @assert size(op) == (N1*N2, N1*N2)
 
     if fast
         # Remove scalar part
         scalar = real(tr(op))
         op = op - (scalar/size(op, 1))*I
 
         # Remove bilinear part
+        gen1 = spin_matrices_of_dim(; N=N1)
+        gen2 = spin_matrices_of_dim(; N=N2)
         bilin = zeros(3, 3)
         for α in 1:3, β in 1:3
             v = kron(gen1[α], gen2[β])
@@ -89,16 +71,24 @@ function decompose_general_coupling(op, gen1, gen2; fast)
         end
         bilin = Mat3(bilin)
 
-        # TODO: Remove biquadratic part
-        u = sum(kron(gen1[α], gen2[α]) for α in 1:3)
-        if norm(op) > 1e-12 && normalize(op) ≈ normalize(u^2 + u/2)
-            @info "Detected scalar biquadratic. Not yet optimized."
+        # Remove biquadratic part
+        biquad = zeros(5, 5)
+        Oi = stevens_matrices_of_dim(2; N=N1)
+        Oj = stevens_matrices_of_dim(2; N=N2)
+        for α in 1:5, β in 1:5
+            v = kron(Oi[α], Oj[β])
+            J = tr(v' * op) / tr(v' * v)
+            @assert imag(J) < 1e-12
+            biquad[α, β] = real(J)
+            op = op - v * bilin[α, β]
         end
+        biquad = SMatrix{5, 5}(biquad)
     else
         scalar = bilin = 0.0
+        biquad = nothing
     end
 
-    return scalar, bilin, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2))
+    return scalar, bilin, biquad, TensorDecomposition(gen1, gen2, svd_tensor_expansion(op, N1, N2))
 end
 
 function Base.zero(::Type{TensorDecomposition})
@@ -133,11 +123,11 @@ function Base.isapprox(op1::TensorDecomposition, op2::TensorDecomposition; kwarg
     return isapprox(op1′, op2′; kwargs...)
 end
 
-function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float64, Mat3}, biquad::Float64, tensordec::TensorDecomposition, bond::Bond)
+function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float64, Mat3}, tensordec::TensorDecomposition, bond::Bond)
     # If `sys` has been reshaped, then operate first on `sys.origin`, which
     # contains full symmetry information.
     if !isnothing(sys.origin)
-        set_pair_coupling_aux!(sys.origin, scalar, bilin, biquad, tensordec, bond)
+        set_pair_coupling_aux!(sys.origin, scalar, bilin, tensordec, bond)
         set_interactions_from_origin!(sys)
         return
     end
@@ -167,7 +157,8 @@ function set_pair_coupling_aux!(sys::System, scalar::Float64, bilin::Union{Float
         for bond′ in all_symmetry_related_bonds_for_atom(sys.crystal, i, bond)
             bilin′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, bilin)
             tensordec′ = transform_coupling_for_bonds(sys.crystal, bond′, bond, tensordec)
-            push_coupling!(ints[i].pair, bond′, scalar, bilin′, biquad, tensordec′)
+            scalar_biquad = 0.0
+            push_coupling!(ints[i].pair, bond′, scalar, bilin′, scalar_biquad, tensordec′)
         end
     end
 end
@@ -199,12 +190,12 @@ function set_pair_coupling!(sys::System{N}, op::AbstractMatrix, bond; fast=true)
         error("System with mode `:dipole_large_S` requires a symbolic operator.")
     end
 
-    gen1 = spin_matrices(spin_irrep_label(sys, bond.i))
-    gen2 = spin_matrices(spin_irrep_label(sys, bond.j))
-    scalar, bilin, tensordec = decompose_general_coupling(op, gen1, gen2; fast)
-    biquad = 0.0
+    N1 = Int(2spin_irrep_label(sys, bond.i)+1)
+    N2 = Int(2spin_irrep_label(sys, bond.j)+1)
+    scalar, bilin, biquad, tensordec = decompose_general_coupling(op, N1, N2; fast)
+    scalar_biquad = 0.0
 
-    set_pair_coupling_aux!(sys, scalar, bilin, biquad, tensordec, bond)
+    set_pair_coupling_aux!(sys, scalar, bilin, scalar_biquad, tensordec, bond)
 end
 
 function set_pair_coupling!(sys::System{N}, fn::Function, bond) where N
@@ -250,15 +241,12 @@ J2 = 2*I + dmvec([0,0,3])
 set_exchange!(sys, J2, bond)
 ```
 """
-function set_exchange!(sys::System{N}, J, bond::Bond; biquad=0, large_S=false) where N
-    large_S && @assert sys.mode == :dipole_large_S
+function set_exchange!(sys::System{N}, J, bond::Bond; scalar_biquad=nothing) where N
+    !isnothing(scalar_biquad) && error("Use `set_pair_coupling!` to assign a biquadratic interaction.")
 
     is_homogeneous(sys) || error("Use `set_exchange_at!` for an inhomogeneous system.")
     bilin = to_float_or_mat3(J)
-    if sys.mode == :dipole
-        bilin, biquad = renormalize_biquad(sys, bond, bilin, biquad)
-    end
-    set_pair_coupling_aux!(sys, 0.0, bilin, Float64(biquad), zero(TensorDecomposition), bond)
+    set_pair_coupling_aux!(sys, 0.0, bilin, zero(TensorDecomposition), bond)
 end
 
 
@@ -308,7 +296,7 @@ end
 
 
 """
-    set_exchange_at!(sys::System, J, site1::Site, site2::Site; biquad=0, large_S=false, offset=nothing)
+    set_exchange_at!(sys::System, J, site1::Site, site2::Site; offset=nothing)
 
 Sets the exchange interaction along the single bond connecting two
 [`Site`](@ref)s, ignoring crystal symmetry. The system must support
@@ -320,21 +308,18 @@ due to possible periodic wrapping. Resolve this ambiguity by passing an explicit
 
 See also [`set_exchange!`](@ref).
 """
-function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=0, large_S=false, offset=nothing) where N
+function set_exchange_at!(sys::System{N}, J, site1::Site, site2::Site; biquad=nothing, offset=nothing) where N
+    !isnothing(biquad) && error("Use `set_pair_coupling_at!` to assign a biquadratic interaction.")
+
     site1 = to_cartesian(site1)
     site2 = to_cartesian(site2)
     bond = sites_to_internal_bond(sys, site1, site2, offset)
 
     is_homogeneous(sys) && error("Use `to_inhomogeneous` first.")
     ints = interactions_inhomog(sys)
 
-    bilin = to_float_or_mat3(J)
-    if sys.mode == :dipole && !large_S
-        bilin, biquad = renormalize_biquad(sys, bond, bilin, biquad)
-    end
-
-    push_coupling!(ints[site1].pair, bond, 0.0, bilin, Float64(biquad), zero(TensorDecomposition))
-    push_coupling!(ints[site2].pair, reverse(bond), 0.0, bilin', Float64(biquad), zero(TensorDecomposition))
+    push_coupling!(ints[site1].pair, bond, 0.0, bilin, zero(TensorDecomposition))
+    push_coupling!(ints[site2].pair, reverse(bond), 0.0, bilin', zero(TensorDecomposition))
 
     return
 end

src/System/Types.jl

@@ -34,7 +34,7 @@ struct PairCoupling
     # procedure in https://arxiv.org/abs/2304.03874
     scalar   :: Float64              # Constant shift
     bilin    :: Union{Float64, Mat3} # Bilinear exchange
-    biquad   :: Float64              # Scalar biquadratic
+    scalar_biquad   :: Float64       # Scalar biquadratic
 
     # General pair interactions, only valid in SU(N) mode
     general  :: TensorDecomposition