With c0 shift

ext/PlottingExt/ViewCrystal.jl

@@ -60,6 +60,57 @@ function propagate_reference_bond_for_cell(cryst, b_ref)
     return reduce(vcat, found)
 end
 
+# Get the quadratic anisotropy as a 3×3 exchange matrix for atom `i` in the
+# chemical cell.
+function quadratic_anisotropy(sys, i)
+    # Extract quadratic Stevens coefficients
+    interactions = isnothing(sys) ? nothing : Sunny.interactions_homog(something(sys.origin, sys))
+    onsite = interactions[i].onsite
+    (c0, c2) = if onsite isa Sunny.HermitianC64
+        Sunny.matrix_to_stevens_coefficients(onsite)[[0, 2]]
+    else
+        @assert onsite isa Sunny.StevensExpansion
+        (onsite.c0, onsite.c2)
+    end
+
+    # Undo RCS renormalization for quadrupolar anisotropy for spin-s
+    if sys.mode == :dipole
+        s = (sys.Ns[i] - 1) / 2
+        c2 /= Sunny.rcs_factors(s)[2]
+    end
+
+    # Stevens quadrupoles expressed as 3×3 bilinears
+    quadrupole_basis = [
+        [1 0 0; 0 -1 0; 0 0 0],    # 𝒪₂₂  = SˣSˣ - SʸSʸ
+        [0 0 1; 0 0 0; 1 0 0] / 2, # 𝒪₂₁  = (SˣSᶻ + SᶻSˣ)/2
+        [-1 0 0; 0 -1 0; 0 0 2],   # 𝒪₂₀  = 2SᶻSᶻ - SˣSˣ - SʸSʸ
+        [0 0 0; 0 0 1; 0 1 0] / 2, # 𝒪₂₋₁ = (SʸSᶻ + SᶻSʸ)/2
+        [0 1 0; 1 0 0; 0 0 0],     # 𝒪₂₋₂ = SˣSʸ + SʸSˣ
+    ]
+
+    # Check consistency with Stevens operators as symbolic polynomials
+    S = spin_matrices(Inf)
+    O = stevens_matrices(Inf)
+    for (b, q) in zip(quadrupole_basis, 2:-1:-2)
+        @assert iszero(S' * b * S - O[2, q])
+    end
+
+    # The c0 coefficient incorporates a factor of S². For quantum spin
+    # operators, S² = s(s+1) I. For the large-s classical limit, S² = s² is a
+    # scalar.
+    S² = if sys.mode == :dipole_uncorrected
+        # Undoes extraction in `operator_to_stevens_coefficients`. Note that
+        # spin magnitude s² is set to κ², as originates from `onsite_coupling`
+        # for p::AbstractPolynomialLike.
+        sys.κs[i]^2
+    else
+        # Undoes extraction in `matrix_to_stevens_coefficients` where 𝒪₀₀ = I.
+        s = (sys.Ns[i]-1) / 2
+        s * (s+1)
+    end
+
+    return c2' * quadrupole_basis + only(c0) * I / S²
+end
 
 # Get the 3×3 exchange matrix for bond `b`
 function exchange_on_bond(interactions, b)
@@ -274,7 +325,7 @@ function is_type_degenerate(cryst, i)
 end
 
 # Construct atom labels for use in DataInspector
-function label_atoms(cryst; ismagnetic)
+function label_atoms(cryst; ismagnetic, sys)
     return map(1:natoms(cryst)) do i
         typ = cryst.types[i]
         rstr = Sunny.fractional_vec3_to_string(cryst.positions[i])
@@ -287,11 +338,17 @@ function label_atoms(cryst; ismagnetic)
             push!(ret, isempty(typ) ? "Position $rstr" : "'$typ' at $rstr")
         end
         if ismagnetic
-            # See similar logic in print_site()
-            refatoms = [b.i for b in Sunny.reference_bonds(cryst, 0.0)]
-            i_ref = Sunny.findfirstval(i_ref -> Sunny.is_related_by_symmetry(cryst, i, i_ref), refatoms)
-            R_site = Sunny.rotation_between_sites(cryst, i, i_ref)
-            push!(ret, Sunny.allowed_g_tensor_string(cryst, i_ref; R_site, prefix="Aniso: ", digits=8, atol=1e-12))
+            if isnothing(sys)
+                # See similar logic in print_site()
+                refatoms = [b.i for b in Sunny.reference_bonds(cryst, 0.0)]
+                i_ref = Sunny.findfirstval(i_ref -> Sunny.is_related_by_symmetry(cryst, i, i_ref), refatoms)
+                R_site = Sunny.rotation_between_sites(cryst, i, i_ref)
+                push!(ret, Sunny.allowed_g_tensor_string(cryst, i_ref; R_site, prefix="Aniso: ", digits=8, atol=1e-12))
+            else
+                ansio = quadratic_anisotropy(sys, i)
+                basis_strs = Sunny.number_to_simple_string.(ansio; digits=3)
+                push!(ret, Sunny.formatted_matrix(basis_strs; prefix="Ansiso: "))
+            end
         end
         join(ret, "\n")
     end
@@ -340,7 +397,7 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy
             # Labels for non-ghost atoms
             inspector_label = nothing
             if !isghost
-                labels = label_atoms(xtal; ismagnetic)[idxs]
+                labels = label_atoms(xtal; ismagnetic, sys)[idxs]
                 inspector_label = (_plot, index, _position) -> labels[index]
             end
 

src/Operators/Stevens.jl

@@ -72,7 +72,8 @@ function stevens_abstract_polynomials(; J, k::Int)
 end
 
 
-# Construct Stevens operators as polynomials in the spin operators.
+# Construct Stevens operators as polynomials in the spin operators. Listed in
+# descending order q = k,..-k.
 function stevens_matrices_of_dim(k::Int; N::Int)
     if k >= N
         return fill(Hermitian(zeros(ComplexF64, N, N)), 2k+1)
@@ -134,6 +135,9 @@ end
 const stevens_αinv = map(inv, stevens_α)
 
 
+# Expands matrix A in Stevens operators. The coefficients are returned as an
+# OffsetArray c[k] with indices k = 0..6. Elements of c[k][:] are the Stevens
+# coefficients in descending order q = k..-k.
 function matrix_to_stevens_coefficients(A::HermitianC64)
     N = size(A,1)
     @assert N == size(A,2)

src/Operators/Symbolic.jl

@@ -91,9 +91,9 @@ end
 
 
 # Extract Stevens operator coefficients from spin polynomial
-function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S)
+function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S²)
     p = expand_in_stevens_operators(p)
-    p = DP.subs(p, spin_squared_symbol => S^2)
+    p = DP.subs(p, spin_squared_symbol => S²)
     return map(stevens_symbols) do 𝒪ₖ
         map(𝒪ₖ) do 𝒪kq
             j = findfirst(==(𝒪kq), DP.monomials(p))

src/System/OnsiteCoupling.jl

@@ -24,8 +24,8 @@ function onsite_coupling(sys, site, p::DP.AbstractPolynomialLike)
         error("Symbolic operator only valid for system with mode `:dipole_uncorrected`.")
     end
 
-    s = sys.κs[site]
-    c = operator_to_stevens_coefficients(p, s)
+    S² = sys.κs[site]^2
+    c = operator_to_stevens_coefficients(p, S²)
 
     # No renormalization here because symbolic polynomials `p` are associated
     # with the large-s limit.

src/System/Types.jl

@@ -1,4 +1,6 @@
-# Stevens function expansion, renormalized for dipole projection
+# Coefficients for an expansion in Stevens functions. Each component ck lists
+# coefficients in descending order q = k,..-k. In :dipole mode, the
+# renormalization factors (`rcs_factor`) are included in these coefficients.
 struct StevensExpansion
     kmax :: Int
     c0 :: SVector{1, Float64}

test/test_operators.jl

@@ -330,11 +330,11 @@ end
     # Test Stevens coefficients extraction
     S = spin_matrices(Inf)
     O = stevens_matrices(Inf)
-    S_mag = π
+    S² = π
     p = S'*S * O[4, 2]
-    c = Sunny.operator_to_stevens_coefficients(p, S_mag)
+    c = Sunny.operator_to_stevens_coefficients(p, S²)
     @test iszero(c[1]) && iszero(c[2]) && iszero(c[3]) && iszero(c[5]) && iszero(c[6])
-    @test c[4] ≈ [0.0, 0.0, S_mag^2, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    @test c[4] ≈ [0.0, 0.0, S², 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
     # Test round trip Stevens -> spin -> Stevens
     c_ref = map(OffsetArrays.OffsetArray(0:6, 0:6)) do k