Consistency in expansion of anisotropy

Project.toml

@@ -1,7 +1,7 @@
 name = "Sunny"
 uuid = "2b4a2ac8-8f8b-43e8-abf4-3cb0c45e8736"
 authors = ["The Sunny team"]
-version = "0.7.4"
+version = "0.7.5"
 
 [deps]
 Brillouin = "23470ee3-d0df-4052-8b1a-8cbd6363e7f0"

docs/src/versions.md

@@ -1,5 +1,13 @@
 # Version History
 
+## v0.7.5
+(In development)
+
+* [`view_crystal`](@ref) shows allowed quadratic anisotropy.
+  [`print_site`](@ref) accepts an optional reference atom `i_ref`, with default
+  of `i`. The optional reference bond `b_ref` of [`print_bond`](@ref) now
+  defaults to `b`.
+
 ## v0.7.4
 (Dec 6, 2024)
 

ext/PlottingExt/ViewCrystal.jl

@@ -36,7 +36,7 @@ function reference_bonds_upto(cryst, nbonds, ndims)
     refbonds = filter(reference_bonds(cryst, max_dist)) do b
         return !(b.i == b.j && iszero(b.n))
     end
-    
+
     # Verify max_dist heuristic
     if length(refbonds) > 10nbonds
         @warn "Found $(length(refbonds)) bonds using max_dist of $max_dist"
@@ -60,14 +60,64 @@ function propagate_reference_bond_for_cell(cryst, b_ref)
     return reduce(vcat, found)
 end
 
+function anisotropy_on_site(sys, i)
+    interactions = isnothing(sys) ? nothing : Sunny.interactions_homog(something(sys.origin, sys))
+    onsite = interactions[i].onsite
+    if onsite isa Sunny.HermitianC64
+        onsite = Sunny.StevensExpansion(Sunny.matrix_to_stevens_coefficients(onsite))
+    end
+    return onsite :: Sunny.StevensExpansion
+end
+
+# Get the quadratic anisotropy as a 3×3 exchange matrix for atom `i` in the
+# chemical cell.
+function quadratic_anisotropy(sys, i)
+    # Get certain Stevens expansion coefficients
+    (; c0, c2) = anisotropy_on_site(sys, i)
 
-# Get the 3×3 exchange matrix for bond `b`
-function exchange_on_bond(interactions, b)
-    isnothing(interactions) && return zero(Sunny.Mat3)
+    # Undo RCS renormalization for quadrupolar anisotropy for spin-s
+    if sys.mode == :dipole
+        s = (sys.Ns[i] - 1) / 2
+        c2 = c2 / Sunny.rcs_factors(s)[2] # Don't mutate c2 in-place!
+    end
+
+    # Stevens quadrupole operators expressed as 3×3 spin 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ˣ
+    ]
+
+    # 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
+
+function coupling_on_bond(interactions, b)
+    isnothing(interactions) && return zero(Mat3)
     pairs = interactions[b.i].pair
     indices = findall(pc -> pc.bond == b, pairs)
-    isempty(indices) && return zero(Sunny.Mat3)
-    return pairs[only(indices)].bilin * Mat3(I)
+    return isempty(indices) ? nothing : pairs[only(indices)]
+end
+
+# Get the 3×3 exchange matrix for bond `b`
+function exchange_on_bond(interactions, b)
+    coupling = coupling_on_bond(interactions, b)
+    return isnothing(coupling) ? zero(Mat3) : coupling.bilin * Mat3(I)
 end
 
 # Get largest exchange interaction scale. For symmetric part, this is the
@@ -100,7 +150,7 @@ function exchange_decomposition(J)
 
     # Now vecs is a pure rotation
     @assert vecs'*vecs ≈ I && det(vecs) ≈ 1
-    
+
     # Quaternion that rotates Cartesian coordinates into principle axes of J.
     axis, angle = Sunny.matrix_to_axis_angle(Mat3(vecs))
     q = iszero(axis) ? Makie.Quaternionf(0,0,0,1) : Makie.qrotation(axis, angle)
@@ -178,7 +228,7 @@ function draw_exchange_geometries(; ax, obs, ionradius, pts, scaled_exchanges)
 end
 
 function draw_bonds(; ax, obs, ionradius, exchange_mag, cryst, interactions, bonds, refbonds, color)
-    
+
     # Map each bond to line segments in global coordinates
     segments = map(bonds) do b
         (; ri, rj) = Sunny.BondPos(cryst, b)
@@ -226,6 +276,10 @@ function draw_bonds(; ax, obs, ionradius, exchange_mag, cryst, interactions, bon
                 dmvecstr = join(Sunny.number_to_simple_string.(dmvec; digits=3), ", ")
                 J_matrix_str *= "\nDM: [$dmvecstr]"
             end
+            c = coupling_on_bond(interactions, b)
+            if !isnothing(c) && (!iszero(c.biquad) || !isempty(c.general.data))
+                J_matrix_str *= "\n  + higher order terms"
+            end
         end
 
         return """
@@ -274,20 +328,33 @@ 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])
         ret = []
 
-        if ismagnetic && is_type_degenerate(cryst, i)
-            c = cryst.classes[i]
-            push!(ret, isempty(typ) ? "Class $c at $rstr" : "'$typ' (class $c) at $rstr")
-        else
-            push!(ret, isempty(typ) ? "Position $rstr" : "'$typ' at $rstr")
-        end
+        (; multiplicity, letter) = Sunny.get_wyckoff(cryst, i)
+        wyckstr = "$multiplicity$letter"
+        typstr = isempty(typ) ? "" : "'$typ', "
+        push!(ret, typstr * "Wyckoff $wyckstr, $rstr")
+
         if ismagnetic
-            # TODO: Show onsite couplings?
+            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
+                aniso = quadratic_anisotropy(sys, i)
+                basis_strs = Sunny.number_to_simple_string.(aniso; digits=3)
+                push!(ret, Sunny.formatted_matrix(basis_strs; prefix="Aniso: "))
+                (; c4, c6) = anisotropy_on_site(sys, i)
+                if !iszero(c4) || !iszero(c6)
+                    push!(ret, "  + higher order terms")
+                end
+            end
         end
         join(ret, "\n")
     end
@@ -336,10 +403,10 @@ 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
-            
+
             # Show dipoles. Mostly consistent with code in plot_spins.
             if !isnothing(sys) && ismagnetic
                 sites = Sunny.position_to_site.(Ref(sys), rs)

src/Operators/Rotation.jl

@@ -143,8 +143,10 @@ end
 Rotates the local quantum operator `A` according to the ``3×3`` rotation matrix
 `R`.
 """
-function rotate_operator(A::HermitianC64, R)
+function rotate_operator(A::AbstractMatrix{ComplexF64}, R)
     isempty(A) && return A
+    A ≈ A' || error("Non-Hermitian operator")
+    A = hermitianpart(A)
     R = convert(Mat3, R)
     N = size(A, 1)
     U = unitary_irrep_for_rotation(R; N)

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)
@@ -162,15 +166,15 @@ function operator_for_stevens_rotation(k, R)
     return real(V)
 end
 
-# Let c denote coefficients of an operator expansion 𝒜 = c† 𝒪. Under the
+# Let c denote coefficients of an operator expansion 𝒜 = cᵀ 𝒪. Under the
 # rotation R, Stevens operators transform as 𝒪 → V 𝒪. Alternatively, we can
 # treat the Stevens operators as fixed, provided the coefficients transform as
-# c† → c† V, or c → V† c.
-function rotate_stevens_coefficients(c, R::Mat3)
+# cᵀ → cᵀ V, or c → Vᵀ c.
+function rotate_stevens_coefficients(c::AbstractVector{Float64}, R::Mat3)
     N = length(c)
     k = Int((N-1)/2)
     V = operator_for_stevens_rotation(k, R)
-    return V' * c
+    return transpose(V) * c
 end
 
 

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))
@@ -102,7 +102,7 @@ function operator_to_stevens_coefficients(p::DP.AbstractPolynomialLike, S)
     end
 end
 
-function rotate_operator(p, R)
+function rotate_operator(p::DP.AbstractPolynomialLike, R)
     𝒮′ = R * [𝒮ˣ, 𝒮ʸ, 𝒮ᶻ]
     DP.subs(p, 𝒮ˣ => 𝒮′[1], 𝒮ʸ => 𝒮′[2], 𝒮ᶻ => 𝒮′[3])
 end

src/Symmetry/AllowedAnisotropy.jl

@@ -16,30 +16,38 @@ function transform_spherical_to_stevens_coefficients(k, c)
     return transpose(stevens_αinv[k]) * c
 end
 
+# The bare Stevens expansion coefficients b::Vector{Float64} are not understood
+# by rotate_operator(). Embed them in a large StevensExpansion object that can
+# be passed to is_anisotropy_valid().
+function is_stevens_expansion_valid(cryst, i, b)
+    N = length(b)
+    c0 = (N == 1) ? b : zeros(1)
+    c2 = (N == 5) ? b : zeros(5)
+    c4 = (N == 9) ? b : zeros(9)
+    c6 = (N == 13) ? b : zeros(13)
+    stvexp = StevensExpansion(6, c0, c2, c4, c6)
+    return is_anisotropy_valid(cryst, i, stvexp)
+end
 
-function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, R::Mat3)
-    # The symmetry operations for the point group at atom i. Each one encodes a
-    # rotation/reflection.
+function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int, R_global::Mat3, atol::Float64)
+    # The symmetry operations for the point group at atom i.
     symops = symmetries_for_pointgroup_of_atom(cryst, i)
 
     # The Wigner D matrices for each symop
     Ds = map(symops) do s
         # R is an orthogonal matrix that transforms positions, x → x′ = R x. It
-        # might or might not include a reflection, i.e., det R = ±1.
-        sR = cryst.latvecs * s.R * inv(cryst.latvecs)
-
-        # Unlike position x, spin S = [Sx, Sy, Sz] is a _pseudo_ vector, which
-        # means that, under reflection, the output gains an additional minus
-        # sign. That is, the orthogonal transformation R applied to spin has the
-        # action, S → S′ = ± R S, where the minus sign corresponds to the case
-        # det(R) = -1. More simply, we may write S′ = Q S, where Q = det(R) R.
-        Q = det(sR) * sR
-
-        # The Wigner D matrix, whose action on a spherical tensor corresponds to
-        # the 3x3 rotation Q (see more below).
+        # might or might not include an inversion, i.e., det R = ±1.
+        R = cryst.latvecs * s.R * inv(cryst.latvecs)
+
+        # Unlike position, spin angular momentum is a pseudo-vector, which means
+        # it is invariant under global inversion. The transformation of spin is
+        # always a pure rotation, S → S′ = Q S, where Q = R det R.
+        Q = det(R) * R
+
+        # The Wigner D matrix corresponding to Q (see more below).
         return unitary_irrep_for_rotation(Q; N=2k+1)
     end
-    
+
     # A general operator in the spin-k representation can be decomposed in the
     # basis of spherical tensors, 𝒜 = ∑_q c_q T_kq, for some coefficients c_q.
     # Spherical tensors transform as T_kq → D^{*}_qq′ T_kq′. Alternatively, we
@@ -53,31 +61,56 @@ function basis_for_symmetry_allowed_anisotropies(cryst::Crystal, i::Int; k::Int,
     # operator 𝒜.
     C = sum(D' for D in Ds)
 
-    # Transform coefficients c to c′ in rotated Stevens operators, T′ = D* T,
-    # where the Wigner D matrix is associated with the rotation R. That is, find
-    # c′ satisfying c′ᵀ T′ = c T. Recall c′ᵀ T′ = (c′ᵀ D*) T = (D† c′)ᵀ T. The
-    # constraint becomes D† c′ = c. Since D is unitary, we have c′ = D c. We
-    # apply this transformation to each column c of C.
-    D = unitary_irrep_for_rotation(R; N=2k+1)
-    C = D * C
-
     # Transform columns c of C to columns b of B such that 𝒜 = cᵀ T = bᵀ 𝒪.
     # Effectively, this reexpresses the symmetry-allowed operators 𝒜 in the
     # basis of Stevens operators 𝒪.
     B = mapslices(C; dims=1) do c
         transform_spherical_to_stevens_coefficients(k, c)
     end
 
+    # Apply a global rotation to the Cartesian coordinate system. Stevens
+    # operators rotate as 𝒪′ = V 𝒪. Coefficients satisfying b′ᵀ 𝒪′ = bᵀ 𝒪
+    # must therefore transform as b′ = V⁻ᵀ b.
+    V = operator_for_stevens_rotation(k, R_global)
+    B = transpose(V) \ B
+
     # If 𝒜 is symmetry-allowed, then its Hermitian and anti-Hermitian parts are
     # independently symmetry-allowed. These are associated with the real and
     # imaginary parts of B. Use the real and imaginary parts of B to construct
     # an over-complete set of symmetry-allowed operators that are guaranteed to
     # be Hermitian. Create a widened real matrix from these two parts, and
     # eliminate linearly-dependent vectors from the column space.
-    B = colspace(hcat(real(B), imag(B)); atol=1e-12)
+    B = colspace(hcat(real(B), imag(B)); atol)
 
     # Find linear combination of columns that sparsifies B
-    return sparsify_columns(B; atol=1e-12)
+    B = sparsify_columns(B; atol)
+
+    # Scale each column to make the final expression prettier
+    return map(eachcol(B)) do b
+        b /= argmax(abs, b)
+        if any(x -> atol < abs(x) < sqrt(atol), b)
+            @info """Found a very small but nonzero expansion coefficient.
+                        This may indicate a slightly misaligned reference frame."""
+        end
+
+        # Magnify by up to 60× if it makes all coefficients integer
+        denoms = denominator.(rationalize.(b; tol=atol))
+        if all(<=(60), denoms)
+            factor = lcm(denominator.(rationalize.(b; tol=atol)))
+            if factor <= 60
+                b .*= factor
+            end
+        end
+
+        # Check that the operator bᵀ 𝒪 satisifies point group symmetries in the
+        # original coordinate system. Multiply by Vᵀ to effectively undo the
+        # global rotation of Stevens operators 𝒪.
+        @assert is_stevens_expansion_valid(cryst, i, transpose(V) * b)
+
+        b
+    end
+
+    return B
 end