begin refactors

src/Operators/Spin.jl

@@ -20,7 +20,6 @@ function spin_matrices(; N::Int)
     return [Sx, Sy, Sz]
 end
 
-
 # Returns ⟨Z|Sᵅ|Z⟩
 @generated function expected_spin(Z::CVec{N}) where N
     S = spin_matrices(; N)
@@ -39,6 +38,41 @@ end
     end
 end
 
+function expected_multipolar_moments(Z::CVec{N}) where N
+    M = Z * Z'
+    Q = zeros(Float64,N,N)
+    for i = 1:N, j = 1:N
+        if i <= j
+            Q[i,j] = real(M[i,j] + M[j,i])
+        else
+            Q[i,j] = imag(M[j,i] - M[i,j])
+        end
+    end
+    Q[1:(N^2 - 1)]
+end
+
+function multipolar_generators_times_Z(dE_dT,Z::CVec{N}) where N
+    out = MVector{N,ComplexF64}(undef)
+    out .= 0
+    li = LinearIndices(zeros(N,N))
+    for i = 1:N, j = 1:N
+        k = li[i,j]
+        if k == N^2
+          break
+        end
+        if i <= j
+            out[i] += Z[j] * dE_dT[k]
+            out[j] += Z[i] * dE_dT[k]
+        else
+            out[i] += im * Z[j] * dE_dT[k]
+            out[j] -= im * Z[i] * dE_dT[k]
+        end
+    end
+    out
+end
+
+
+
 # Find a ket (up to an irrelevant phase) that corresponds to a pure dipole.
 # TODO, we can do this faster by using the exponential map of spin operators,
 # expressed as a polynomial expansion,

src/System/Interactions.jl

@@ -131,9 +131,15 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N
         (; onsite, pair) = interactions_inhomog(sys)[site]
     end
 
-    s₀ = dipoles[site]
     Z₀ = coherents[site]
+
+    s₀ = dipoles[site]
     Δs = s - s₀
+
+    q = expected_multipolar_moments(Z)
+    q₀ = expected_multipolar_moments(Z₀)
+    Δq = q - q₀
+
     ΔE = 0.0
 
     # Zeeman coupling to external field
@@ -150,25 +156,14 @@ function local_energy_change(sys::System{N}, site, state::SpinState) where N
         ΔE += real(dot(Z, Λ, Z) - dot(Z₀, Λ, Z₀))
     end
 
-    # Quadratic exchange matrix
+    # Quadratic (in the SU(N) generators) exchange matrix
     for coupling in pair
         (; bond) = coupling
         cellⱼ = offsetc(to_cell(site), bond.n, latsize)
-        sⱼ = dipoles[cellⱼ, bond.j]
-
-        # Bilinear
-        J = coupling.bilin
-        ΔE += dot(Δs, J, sⱼ)
-
-        # Biquadratic
-        if !iszero(coupling.biquad)
-            J = coupling.biquad
-            if sys.mode == :dipole
-                ΔE += J * ((s⋅sⱼ)^2 - (s₀⋅sⱼ)^2)
-            elseif sys.mode == :SUN
-                error("Biquadratic currently unsupported in SU(N) mode.") 
-            end
-        end
+        qⱼ = expected_multipolar_moments(coherents[cellⱼ, bond.j])
+
+        J = coupling.matrix
+        ΔE += dot(Δq, J, qⱼ)
     end
 
     # Long-range dipole-dipole
@@ -247,22 +242,11 @@ function energy_aux(sys::System{N}, ints::Interactions, i::Int, cells, foreachbo
     end
 
     foreachbond(ints.pair) do coupling, site1, site2
-        sᵢ = dipoles[site1]
-        sⱼ = dipoles[site2]
-
-        # Bilinear
-        J = coupling.bilin
-        E += dot(sᵢ, J, sⱼ)
-
-        # Biquadratic
-        if !iszero(coupling.biquad)
-            J = coupling.biquad
-            if sys.mode == :dipole
-                E += J * (sᵢ⋅sⱼ)^2
-            elseif sys.mode == :SUN
-                error("Biquadratic currently unsupported in SU(N) mode.")
-            end
-        end
+        qᵢ = expected_multipolar_moments(coherents[site1])
+        qⱼ = expected_multipolar_moments(coherents[site2])
+
+        J = coupling.matrix
+        E += dot(qᵢ, J, qⱼ)
     end
 
     return E
@@ -273,7 +257,7 @@ end
 # contributions from Zeeman coupling, bilinear exchange, and long-range
 # dipole-dipole. Excluded terms include onsite coupling, and general pair
 # coupling (biquadratic and beyond).
-function set_energy_grad_dipoles!(∇E, dipoles::Array{Vec3, 4}, sys::System{N}) where N
+function set_energy_grad_dipoles!(∇E, coherents, dipoles::Array{Vec3, 4}, sys::System{N}) where N
     (; crystal, latsize, extfield, ewald) = sys
 
     fill!(∇E, zero(Vec3))
@@ -288,12 +272,12 @@ function set_energy_grad_dipoles!(∇E, dipoles::Array{Vec3, 4}, sys::System{N})
         if is_homogeneous(sys)
             # Interaction is the same at every cell
             interactions = sys.interactions_union[i]
-            set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, eachcell(sys), homog_bond_iterator(latsize))
+            set_energy_grad_dipoles_aux!(∇E, coherents, dipoles, interactions, sys, i, eachcell(sys), homog_bond_iterator(latsize))
         else
             for cell in eachcell(sys)
                 # There is a different interaction at every cell
                 interactions = sys.interactions_union[cell,i]
-                set_energy_grad_dipoles_aux!(∇E, dipoles, interactions, sys, i, (cell,), inhomog_bond_iterator(latsize, cell))
+                set_energy_grad_dipoles_aux!(∇E, coherents, dipoles, interactions, sys, i, (cell,), inhomog_bond_iterator(latsize, cell))
             end
         end
     end
@@ -306,7 +290,7 @@ end
 # Calculate the energy gradient `∇E' for the sublattice `i' at all elements of
 # `cells`. The function `foreachbond` enables efficient iteration over
 # neighboring cell pairs (without double counting).
-function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells, foreachbond) where N
+function set_energy_grad_dipoles_aux!(∇E, coherents, dipoles::Array{Vec3, 4}, ints::Interactions, sys::System{N}, i::Int, cells, foreachbond) where N
     # Single-ion anisotropy only contributes in dipole mode. In SU(N) mode, the
     # anisotropy matrix will be incorporated directly into ℌ.
     if N == 0
@@ -318,24 +302,12 @@ function set_energy_grad_dipoles_aux!(∇E, dipoles::Array{Vec3, 4}, ints::Inter
     end
 
     foreachbond(ints.pair) do coupling, site1, site2
-        sᵢ = dipoles[site1]
-        sⱼ = dipoles[site2]
-
-        # Bilinear
-        J = coupling.bilin
-        ∇E[site1] += J  * sⱼ
-        ∇E[site2] += J' * sᵢ
-
-        # Biquadratic
-        if !iszero(coupling.biquad)
-            J = coupling.biquad
-            if sys.mode == :dipole
-                ∇E[site1] += J * 2sⱼ*(sᵢ⋅sⱼ)
-                ∇E[site2] += J * 2sᵢ*(sᵢ⋅sⱼ)
-            elseif sys.mode == :SUN
-                error("Biquadratic currently unsupported in SU(N) mode.")
-            end
-        end
+        qᵢ = expected_multipolar_moments(coherents[site1])
+        qⱼ = expected_multipolar_moments(coherents[site2])
+
+        J = coupling.matrix
+        ∇E[site1] += J  * qⱼ
+        ∇E[site2] += J' * qᵢ
     end
 end
 
@@ -350,20 +322,22 @@ function set_energy_grad_coherents!(HZ, Z, sys::System{N}) where N
     # or the general pair couplings, which must be handled in a more general
     # way.
     dE_ds, dipoles = get_dipole_buffers(sys, 2)
-    @. dipoles = expected_spin(Z)
-    set_energy_grad_dipoles!(dE_ds, dipoles, sys)
+    @. multipoles = expected_spin(Z)
+    set_energy_grad_multipoles!(dE_ds, multipoles, sys)
 
     if is_homogeneous(sys)
         ints = interactions_homog(sys)
         for site in eachsite(sys)
             Λ = ints[to_atom(site)].onsite :: Matrix
-            HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site])
+            #HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site])
+            HZ[site] = Λ * Z[site] + multipolar_generators_times_Z(dE_dT[site],Z[site])
         end
     else
         ints = interactions_inhomog(sys)
         for site in eachsite(sys)
             Λ = ints[site].onsite :: Matrix
-            HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site])
+            #HZ[site] = mul_spin_matrices(Λ, dE_ds[site], Z[site])
+            HZ[site] = Λ * Z[site] + multipolar_generators_times_Z(dE_dT[site],Z[site])
         end 
     end
 
@@ -373,7 +347,7 @@ end
 # Internal testing functions
 function energy_grad_dipoles(sys::System{N}) where N
     ∇E = zero(sys.dipoles)
-    set_energy_grad_dipoles!(∇E, sys.dipoles, sys)
+    set_energy_grad_dipoles!(∇E, sys.coherents, sys.dipoles, sys)
     return ∇E
 end
 function energy_grad_coherents(sys::System{N}) where N

src/System/Types.jl

@@ -23,11 +23,8 @@ struct PairCoupling
 
     # In :dipole mode, these will be renormalized couplings following
     # the procedure in https://arxiv.org/abs/2304.03874
-    bilin    :: Union{Float64, Mat3} # Bilinear exchange
-    biquad   :: Float64              # Scalar biquadratic
+    matrix :: Matrix{Float64}
 
-    # General pair interactions, only valid in SU(N) mode
-    # general  :: Vector{Tuple{Hermitian{ComplexF64}, Hermitian{ComplexF64}}}
     # TODO: update clone_interactions(), set_interactions_from_origin!
 end