Progress

src/System/SingleIonAnisotropy.jl

@@ -245,7 +245,8 @@ function set_anisotropy_at!(sys::System{N}, op::Union{DP.AbstractPolynomialLike,
 end
 
 
-# Evaluate a given linear combination of Stevens operators for classical spin s
+# Evaluate a given linear combination of Stevens operators where each spin
+# operator is replaced by its dipole expectation value.
 function energy_and_gradient_for_classical_anisotropy(s::Vec3, stvexp::StevensExpansion)
     (; kmax, c2, c4, c6) = stvexp
 
@@ -283,13 +284,13 @@ function energy_and_gradient_for_classical_anisotropy(s::Vec3, stvexp::StevensEx
     Jp⁴ = Jp²*Jp²
     Jz⁴ = Jz²*Jz²
 
-    A = (35Jz⁴ - (30X)Jz² + (3X²),
-        7Jz³ - (3X)Jz¹,
-        7Jz² - (X),
+    A = (35Jz⁴ - (30X-25)Jz² + 3X² - 6X,
+        7Jz³ - (3X+1)Jz¹,
+        7Jz² - X - 5,
         Jz¹,
         1)
-    dA_dz = (140Jz³ - (60X)Jz¹,
-            21Jz² - 3X,
+    dA_dz = (140Jz³ - (60X-50)Jz¹,
+            21Jz² - (3X+1),
             14Jz¹,
             1)
     E +=        (c4[1]*real(Jp⁴)+c4[9]*imag(Jp⁴))A[5] +
@@ -316,17 +317,17 @@ function energy_and_gradient_for_classical_anisotropy(s::Vec3, stvexp::StevensEx
     Jp⁶ = Jp³*Jp³
     Jz⁶ = Jz³*Jz³
 
-    A = (231Jz⁶ - (315X)Jz⁴ + (105X²)Jz² - (5X³),
-        33Jz⁵ - (30X)Jz³ + (5X²)Jz¹,
-        33Jz⁴ - (18X)Jz² + (X²),
-        11Jz³ - (3X)Jz¹,
-        11Jz² - (X),
+    A = (231Jz⁶ - (315X-735)Jz⁴ + (105X²-525X+294)Jz² - 5X³ + 40X² - 60X,
+        33Jz⁵ - (30X-15)Jz³ + (5X²-10X+12)Jz¹,
+        33Jz⁴ - (18X+123)Jz² + X² + 10X + 102,
+        11Jz³ - (3X+59)Jz¹,
+        11Jz² - X - 38,
         Jz¹,
         1)
-    dA_dz = (1386Jz⁵ - (1260X)Jz³ + (210X²)Jz¹,
-            165Jz⁴ - (90X)Jz² + 5X²,
-            132Jz³ - (36X)Jz¹,
-            33Jz² - 3X,
+    dA_dz = (1386Jz⁵ - (1260X-2940)Jz³ + (210X²-1050X+588)Jz¹,
+            165Jz⁴ - (90X-45)Jz² + (5X²-10X+12),
+            132Jz³ - (36X+246)Jz¹,
+            33Jz² - (3X+59),
             22Jz¹,
             1)
     E +=        (c6[1]*real(Jp⁶)+c6[13]*imag(Jp⁶))A[7] +

test/test_operators.jl

@@ -209,40 +209,34 @@ end
 
 
 @testitem "Stevens evaluation" begin
-    using LinearAlgebra
+    using LinearAlgebra, FiniteDifferences
 
-    # Test fast evaluation of Stevens functions
-    let
-        import DynamicPolynomials as DP
-
-        s = randn(Sunny.Vec3)
-        p = randn(5)' * Sunny.stevens_operator_symbols[2] + 
-            randn(9)' * Sunny.stevens_operator_symbols[4] +
-            randn(13)' * Sunny.stevens_operator_symbols[6]
-        (_, c2, _, c4, _, c6) = Sunny.operator_to_classical_stevens_coefficients(p, 1.0)
-
-        p_classical = Sunny.operator_to_classical_polynomial(p)
-        grad_p_classical = DP.differentiate(p_classical, Sunny.spin_classical_symbols)
-
-        E_ref = p_classical(Sunny.spin_classical_symbols => s)
-        gradE_ref = [g(Sunny.spin_classical_symbols => s) for g = grad_p_classical]
-
-        stvexp = Sunny.StevensExpansion(c2, c4, c6)
-        E, gradE = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp)
-
-        @test E ≈ E_ref
-
-        # Above, when calculating gradE_ref, the value X = |S|^2 is treated
-        # as varying with S, such that dX/dS = 2S. Conversely, when calculating
-        # gradE, the value X is treated as a constant, such that dX/dS = 0. In
-        # practice, gradE will be used to drive spin dynamics, for which |S| is
-        # constant, and the component of gradE parallel to S will be projected
-        # out anyway. Therefore we only need agreement between the parts of
-        # gradE and gradE_ref that are perpendicular to S.
-        gradE_ref -= (gradE_ref⋅s)*s / (s⋅s) # Orthogonalize to s
-        gradE -= (gradE⋅s)*s / (s⋅s)         # Orthogonalize to s
-        @test gradE_ref ≈ gradE
+    s = randn(Sunny.Vec3)
+    c2 = randn(5)
+    c4 = randn(9)
+    c6 = randn(13)
+
+    function eval_anisotropy(s)
+        c2' * Sunny.stevens_abstract_polynomials(; J=s, k=2) +
+        c4' * Sunny.stevens_abstract_polynomials(; J=s, k=4) +
+        c6' * Sunny.stevens_abstract_polynomials(; J=s, k=6) |> real
     end
+
+    gradE_ref = grad(central_fdm(5, 1), eval_anisotropy, s)[1]
+
+    stvexp = Sunny.StevensExpansion(c2, c4, c6)
+    E, gradE = Sunny.energy_and_gradient_for_classical_anisotropy(s, stvexp)
+
+    # Above, when calculating gradE_ref, the value X = |S|^2 is treated
+    # as varying with S, such that dX/dS = 2S. Conversely, when calculating
+    # gradE, the value X is treated as a constant, such that dX/dS = 0. In
+    # practice, gradE will be used to drive spin dynamics, for which |S| is
+    # constant, and the component of gradE parallel to S will be projected
+    # out anyway. Therefore we only need agreement between the parts of
+    # gradE and gradE_ref that are perpendicular to S.
+    gradE_ref -= (gradE_ref⋅s)*s / (s⋅s) # Orthogonalize to s
+    gradE -= (gradE⋅s)*s / (s⋅s)         # Orthogonalize to s
+    @assert gradE_ref ≈ gradE
 end