Refactor single ion anisotropies (#93)

* Refactor single ion anisotropies

- Replace `set_anisotropy!` with `set_onsite_coupling!`. The new function
  expects a matrix built from `stevens_operators` or  `spin_operators`. Provide
  `print_stevens_expansion` to decompose polynomials of spin operators
  into a linear combination of Stevens operators.

- Make DynamicPolynomials an optional dependency. This provides symbolic
  functions `print_classical_stevens_expansion` and
  `print_classical_spin_polynomial`.

- Functions to create tensor products of operators.

Project.toml

@@ -7,7 +7,6 @@ version = "0.4.4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CrystalInfoFramework = "6007d9b0-c6b2-11e8-0510-1d10e825f3f1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FilePathsBase = "48062228-2e41-5def-b9a4-89aafe57970f"
 Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
@@ -29,7 +28,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Colors = "0.12.8"
 CrystalInfoFramework = "0.5.0"
 DataStructures = "0.18.13"
-DynamicPolynomials = "0.4.6"
 FFTW = "1.4.5"
 FilePathsBase = "0.9.18"
 Inflate = "0.1.2"

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 Formatting = "59287772-0a20-5a39-b81b-1366585eb4c0"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"

docs/make.jl

@@ -2,6 +2,9 @@
 
 using Literate, Documenter, Sunny
 
+import DynamicPolynomials # get symbolic functions
+import GLMakie # get plotting functions
+
 execute = true # set `false` to disable cell evaluation
 
 example_names = ["fei2_tutorial", "powder_averaging", "ising2d"]

docs/src/quick-start.md

@@ -105,7 +105,7 @@ The next steps are typically the following
    crystal unit cells.
 2. Add interactions to the system using functions like
    [`set_external_field!`](@ref), [`set_exchange!`](@ref), and
-   [`set_anisotropy!`](@ref).
+   [`set_onsite_coupling!`](@ref).
 3. Perform Monte Carlo simulation to equilibrate the spin configuration. Options
    include the continuous [`Langevin`](@ref) dynamics, or single-spin flip
    updates with [`LocalSampler`](@ref). The former can efficiently handle

docs/src/versions.md

@@ -1,3 +1,23 @@
+# Version 0.5.0
+
+This version includes many **breaking changes**.
+
+Replace `set_anisotropy!` with a new function [`set_onsite_coupling!`](@ref)
+(and similarly [`set_onsite_coupling_at!`](@ref)). The latter expects an
+explicit matrix representation for the local Hamiltonian. This can be
+constructed, e.g., as a linear combination of [`stevens_operators`](@ref), or as
+a polynomial of [`spin_operators`](@ref). To understand the mapping between
+these two, the new function [`print_stevens_expansion`](@ref) acts on an
+arbitrary local operator.
+
+Symbolic representations of operators are now hidden unless the package
+`DynamicPolynomials` is explicitly loaded by the user. The functionality of
+`print_anisotropy_as_stevens` has been replaced with
+[`print_classical_stevens_expansion`](@ref), while
+`print_anisotropy_as_classical_spins` has become
+[`print_classical_spin_polynomial`](@ref).
+
+
 # Version 0.4.3
 
 **Experimental** support for linear [`SpinWaveTheory`](@ref), implemented in
@@ -50,7 +70,7 @@ This update includes many breaking changes, and is missing some features of
 
 ### Creating a spin `System`
 
-`SpinSystem` has been renamed [`System`](@ref). Its constructor now has the form,
+Rename `SpinSystem` to [`System`](@ref). Its constructor now has the form,
 
 ```julia
 System(crystal, latsize, infos, mode)
@@ -66,21 +86,21 @@ The parameter `mode` is one of `:SUN` or `:dipole`.
 
 Interactions are now added mutably to an existing `System` using the following
 functions: [`set_external_field!`](@ref), [`set_exchange!`](@ref),
-[`set_anisotropy!`](@ref), [`enable_dipole_dipole!`](@ref).
+[`set_onsite_coupling!`](@ref), [`enable_dipole_dipole!`](@ref).
 
 As a convenience, one can use [`dmvec(D)`](@ref) to convert a DM vector to a
 $3×3$ antisymmetric exchange matrix.
 
 Fully general single-ion anisotropy is now possible. The function
-[`set_anisotropy!`](@ref) expects the single ion anisotropy to be expressed as a
+[`set_onsite_coupling!`](@ref) expects the single ion anisotropy to be expressed as a
 polynomial in symbolic spin operators [`𝒮`](@ref), or as a linear combination
 of symbolic Stevens operators [`𝒪`](@ref). For example, an easy axis anisotropy
 in the direction `n` may be written `D*(𝒮⋅n)^2`.
 
 Stevens operators `𝒪[k,q]` admit polynomial expression in spin operators
 `𝒮[α]`. Conversely, a polynomial of spin operators can be expressed as a linear
 combination of Stevens operators. To see this expansion use
-[`print_anisotropy_as_stevens`](@ref).
+`print_anisotropy_as_stevens`.
 
 
 ### Inhomogeneous field

examples/fei2_tutorial.jl

@@ -25,7 +25,6 @@
 # Begin by importing `Sunny` and `GLMakie`, a plotting package.
 
 using Sunny, GLMakie
-#md Makie.inline!(true); #hide
 #nb Sunny.offline_viewers()  # Inject Javascript code for additional plotting capabilities 
 
 # If you see an error `Package <X> not found in current path`, add the package
@@ -145,18 +144,17 @@ set_exchange!(sys, [J′2apm 0.0    0.0;
                     0.0    J′2apm 0.0;
                     0.0    0.0    J′2azz], Bond(1,1,[1,2,1]))
 
-# The function [`set_anisotropy!`](@ref) assigns a single-ion anisotropy. It
-# takes an abstract operator and an atom index. The operator may be a polynomial
-# of spin operators or a linear combination of Stevens operators. Sunny provides
-# special symbols for their construction: [`𝒮`](@ref) is a vector of the three
-# spin operators and [`𝒪`](@ref) are the symbolic Stevens operators. Here we
-# construct an easy-axis anisotropy.
+# The function [`set_onsite_coupling!`](@ref) assigns a single-ion anisotropy
+# operator. It can be constructed, e.g., from the matrices given by
+# [`spin_operators`](@ref) or [`stevens_operators`](@ref). Here we construct an
+# easy-axis anisotropy along the direction $\hat{z}$.
 
 D = 2.165
-set_anisotropy!(sys, -D*𝒮[3]^2, 1)
+S = spin_operators(sys, 1)
+set_onsite_coupling!(sys, -D*S[3]^2, 1)
 
 # Any anisotropy operator can be converted to a linear combination of Stevens
-# operators with [`print_anisotropy_as_stevens`](@ref).
+# operators with [`print_stevens_expansion`](@ref).
 
 # # Calculating structure factor intensities
 # In the remainder of this tutorial, we will examine Sunny's tools for

examples/powder_averaging.jl

@@ -7,7 +7,6 @@
 
 using Sunny, GLMakie
 using Statistics: mean
-#md Makie.inline!(true); #hide
 
 dims = (8,8,8)               # Lattice dimensions
 seed = 1                     # RNG seed for repeatable behavior 

src/Integrators.jl

@@ -146,7 +146,7 @@ end
 
 # Returns (Λ + B⋅S) Z
 @generated function mul_spin_matrices(Λ, B::Sunny.Vec3, Z::Sunny.CVec{N}) where N
-    S = Sunny.spin_matrices(N)
+    S = spin_matrices(; N)
     out = map(1:N) do i
         out_i = map(1:N) do j
             terms = Any[:(Λ[$i,$j])]

src/Operators/Rotation.jl

@@ -61,7 +61,7 @@ end
 function unitary_for_rotation(R::Mat3; N::Int)
     !(R'*R ≈ I)   && error("Not an orthogonal matrix, R = $R.")
     !(det(R) ≈ 1) && error("Matrix includes a reflection, R = $R.")
-    S = spin_matrices(N)
+    S = spin_matrices(; N)
     n, θ = axis_angle(R)
     return exp(-im*θ*(n'*S))
 end

src/Operators/Spin.jl

@@ -1,5 +1,11 @@
-# Construct spin operators, i.e. generators of su(2), of dimension N
-function spin_matrices(N::Int)
+"""
+    spin_matrices(; N)
+
+Constructs the three spin operators, i.e. the generators of SU(2), in the
+`N`-dimensional irrep. See also [`spin_operators`](@ref), which determines the
+appropriate value of `N` for a given site index.
+"""
+function spin_matrices(; N::Int)
     if N == 0
         return fill(zeros(ComplexF64,0,0), 3)
     end
@@ -17,7 +23,7 @@ end
 
 # Returns ⟨Z|Sᵅ|Z⟩
 @generated function expected_spin(Z::CVec{N}) where N
-    S = spin_matrices(N)
+    S = spin_matrices(; N)
     elems_x = SVector{N-1}(diag(S[1], 1))
     elems_z = SVector{N}(diag(S[3], 0))
     lo_ind = SVector{N-1}(1:N-1)
@@ -39,7 +45,7 @@ end
 # http://www.emis.de/journals/SIGMA/2014/084/
 ket_from_dipole(_::Vec3, ::Val{0}) :: CVec{0} = zero(CVec{0})
 function ket_from_dipole(dip::Vec3, ::Val{N}) :: CVec{N} where N
-    S = spin_matrices(N) 
+    S = spin_matrices(; N)
     λs, vs = eigen(dip' * S)
     return CVec{N}(vs[:, argmax(real.(λs))])
 end

src/Operators/Stevens.jl

@@ -77,7 +77,7 @@ function stevens_matrices(k::Int; N::Int)
     if k >= N
         return fill(zeros(ComplexF64, N, N), 2k+1)
     else
-        return stevens_abstract_polynomials(; J=spin_matrices(N), k)
+        return stevens_abstract_polynomials(; J=spin_matrices(; N), k)
     end
 end
 
@@ -167,3 +167,66 @@ function rotate_stevens_coefficients(c, R::Mat3)
     @assert norm(imag(c′)) < 1e-12
     return real(c′)
 end
+
+
+# Builds matrix representation of Stevens operators on demand
+struct StevensMatrices
+    N::Int
+end
+function Base.getindex(this::StevensMatrices, k::Int, q::Int)
+    k < 0  && error("Stevens operators 𝒪[k,q] require k >= 0.")
+    k > 6  && error("Stevens operators 𝒪[k,q] currently require k <= 6.")
+    !(-k <= q <= k) && error("Stevens operators 𝒪[k,q] require -k <= q <= k.")
+    if k == 0
+        return Matrix{ComplexF64}(I, this.N, this.N)
+    else
+        # Stevens operators are stored in descending order: k, k-1, ... -k.
+        # Therefore q=k should be the first index, q_idx=1.
+        q_idx = k - q + 1
+        # Build all (2k+1) matrices and extract the desired one. This seems a
+        # bit wasteful, but presumably this function won't be called in
+        # performance sensitive contexts.
+        return stevens_matrices(k; this.N)[q_idx]
+    end
+end
+
+
+"""
+    function print_stevens_expansion(op)
+
+Prints a local Hermitian operator as a linear combination of Stevens operators.
+This function works on explicit matrix representations. The analogous function
+in the large-``S`` classical limit is
+[`print_classical_stevens_expansion`](@ref).
+
+# Examples
+
+```julia
+S = spin_matrices(N=5)
+print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
+# Prints: (1/20)𝒪₄₀ + (1/4)𝒪₄₄ + 102/5
+```
+"""
+function print_stevens_expansion(op::Matrix{ComplexF64})
+    op ≈ op' || error("Requires Hermitian operator")
+    terms = String[]
+
+    # Decompose op into Stevens coefficients
+    c = matrix_to_stevens_coefficients(op)
+    for k in 1:6
+        for (c_km, m) in zip(reverse(c[k]), -k:k)
+            abs(c_km) < 1e-12 && continue
+            push!(terms, *(coefficient_to_math_string(c_km), "𝒪", int_to_underscore_string.((k,m))...))
+        end
+    end
+
+    # Handle linear shift specially
+    c_00 = real(tr(op))/size(op, 1)
+    abs(c_00) > 1e-12 && push!(terms, number_to_math_string(c_00))
+
+    # Concatenate with plus signs
+    str = join(terms, " + ")
+    # Remove redundant plus signs and print
+    str = replace(str, "+ -" => "- ")
+    println(str)
+end