diff --git a/.gitignore b/.gitignore
index cce9fc8a7..867b24065 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
docs/build/
+docs/contributed-tmp/
docs/src/examples/
Manifest.toml
profile.pb.gz
@@ -6,4 +7,4 @@ profile.pb.gz
.vscode
*.ipynb
*.lyx~
-*.lyx#
\ No newline at end of file
+*.lyx#
diff --git a/README.md b/README.md
index 822be18e6..80fd4c398 100644
--- a/README.md
+++ b/README.md
@@ -13,7 +13,7 @@
Sunny is a Julia package for modeling atomic-scale magnetism using classical spin dynamics with quantum corrections. It provides powerful tools to estimate dynamical structure factor intensities, $\mathcal{S}(𝐪,ω)$, enabling quantitative comparison with experimental scattering data, e.g., neutrons or x-rays.
-A unique feature of Sunny is its treatment of spins as [SU(_N_) coherent states](https://doi.org/10.48550/arXiv.2106.14125). Each quantum spin-_S_ state is a superposition of $N=2S+1$ levels, and evolves under unitary, SU(_N_) transformations. Through neglect of entanglement, the formalism allows to generalize the Landau-Lifshitz dynamics of spin dipoles to a dynamics of spin multipoles. The theory becomes especially useful for modeling materials with strong single-ion anisotropy effects (see our [FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/fei2_tutorial.html)). In the future, the theory could also be used to model explicit spin-orbit coupling, or 'units' of locally entangled spins.
+A unique feature of Sunny is its treatment of spins as [SU(_N_) coherent states](https://doi.org/10.48550/arXiv.2106.14125). Each quantum spin-_S_ state is a superposition of $N=2S+1$ levels, and evolves under unitary, SU(_N_) transformations. Through neglect of entanglement, the formalism allows to generalize the Landau-Lifshitz dynamics of spin dipoles to a dynamics of spin multipoles. The theory becomes especially useful for modeling materials with strong single-ion anisotropy effects (see our [FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/01_LSWT_SU3_FeI2.html)). In the future, the theory could also be used to model explicit spin-orbit coupling, or 'units' of locally entangled spins.
At low-temperatures, Sunny supports the usual linear spin wave theory for spin dipoles, and its ['multi-boson' generalization](https://doi.org/10.48550/arXiv.1307.7731). At finite temperatures, the full classical dynamics (with quantum correction factors) may be preferable to capture strongly nonlinear effects. The [coupling of SU(_N_) spins to a thermal bath](https://doi.org/10.48550/arXiv.2209.01265) also makes possible the study of various non-equilibrium dynamics, e.g., thermal transport, pump-probe experiments, and spin-glass relaxation.
@@ -22,7 +22,7 @@ Sunny provides a number of tools to facilitate the specification and solution of
## Try it out!
-To see Sunny in action, a good starting point is our **[FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/fei2_tutorial.html)**. This compound includes effective spin-1 moments with strong easy-axis anisotropy, and exemplifies the power of simulating SU(3) coherent states.
+A good starting point is our **[FeI₂ tutorial](https://sunnysuite.github.io/Sunny.jl/dev/examples/01_LSWT_SU3_FeI2.html)**. This compound includes effective spin-1 moments with strong easy-axis anisotropy, and exemplifies the power of simulating SU(3) coherent states.
@@ -38,7 +38,7 @@ Sunny is inspired by SpinW, especially regarding symmetry analysis, model specif
| -- | -- | -- | -- |
| Symmetry-guided modeling | ❌ | ✅ | ✅ |
| Interactive graphics | ❌ | ✅ | ✅ |
-| General single-ion anisotropy | ✅ | ❌ | ✅ |
+| General spin-multipole interactions | ✅ | ❌ | ✅ |
| [Interaction renormalization for dipoles](https://arxiv.org/abs/2304.03874) | ❌ | ❌ | ✅ |
| [Multi-flavor spin wave theory](https://arxiv.org/abs/1307.7731) | ✅ | ❌ | ✅ |
| [Classical SU(_N_) spin dynamics](https://arxiv.org/abs/2209.01265) | ❌ | ❌ | ✅ |
diff --git a/docs/make.jl b/docs/make.jl
index 5df3ffc42..3ba988f48 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -3,7 +3,7 @@
import Literate, Documenter, Git
using Sunny, GLMakie, WriteVTK # Load packages to enable Documenter references
-draft = false # set `true` to disable cell evaluation
+isdraft = false # set `true` to disable cell evaluation
# Remove existing Documenter `build` directory
build_path = joinpath(@__DIR__, "build")
@@ -13,10 +13,16 @@ notebooks_path = joinpath(build_path, "assets", "notebooks")
scripts_path = joinpath(build_path, "assets", "scripts")
mkpath.([notebooks_path, scripts_path])
+ismarkdown(name) = splitext(name)[2] == ".md"
+isjulia(name) = splitext(name)[2] == ".jl"
+
function build_examples(example_sources, destdir)
assetsdir = joinpath(fill("..", length(splitpath(destdir)))..., "assets")
+ destpath = joinpath(@__DIR__, "src", destdir)
+ isdir(destpath) && rm(destpath; recursive=true)
+
# Transform each Literate source file to Markdown for subsequent processing by
# Documenter.
for source in example_sources
@@ -33,8 +39,7 @@ function build_examples(example_sources, destdir)
""" * str
end
# Write to `src/$destpath/$name.md`
- dest = joinpath(@__DIR__, "src", destdir)
- Literate.markdown(source, dest; preprocess, credit=false)
+ Literate.markdown(source, destpath; preprocess, credit=false)
end
# Create Jupyter notebooks and Julia script for each Literate example. These
@@ -59,18 +64,12 @@ function build_examples(example_sources, destdir)
end
end
-
function prepare_contributed()
- function is_markdown(name)
- if split(name, ".")[end] == "md"
- return true
- end
- false
- end
-
- # Perform a sparse checkout of the `build` directory from SunnyContributed.
- # This directory contains the markdown files and images generated with
- # Literate on the SunnyContributed repo.
+ # Perform a sparse checkout of the `build` directory from SunnyContributed.
+ # This directory contains the markdown files and images generated with
+ # Literate on the SunnyContributed repo. TODO: If directory exists, should
+ # we just `git pull` instead for speed?
+ isdir("contributed-tmp") && rm("contributed-tmp"; recursive=true, force=true)
mkpath("contributed-tmp")
cd("contributed-tmp")
run(`$(Git.git()) init`)
@@ -89,21 +88,19 @@ function prepare_contributed()
# Generate the base names for each contributed markdown file and prepare
# paths for Documenter (relative to `src/`)
- contrib_names = filter(is_markdown, contrib_files)
+ contrib_names = filter(ismarkdown, contrib_files)
return [joinpath("examples", "contributed", name) for name in contrib_names]
end
-
-
-example_names = ["fei2_tutorial", "out_of_equilibrium", "powder_averaging", "fei2_classical", "ising2d"]
-example_sources = [pkgdir(Sunny, "examples", "$name.jl") for name in example_names]
+example_sources = filter(isjulia, readdir(abspath(pkgdir(Sunny, "examples")), join=true))
+example_names = [splitext(basename(src))[1] for src in example_sources]
example_mds = build_examples(example_sources, "examples")
-spinw_names = ["08_Kagome_AFM", "15_Ba3NbFe3Si2O14"]
-spinw_sources = [pkgdir(Sunny, "examples", "spinw_ports", "$name.jl") for name in spinw_names]
+spinw_sources = filter(isjulia, readdir(abspath(pkgdir(Sunny, "examples", "spinw_tutorials")), join=true))
+spinw_names = [splitext(basename(src))[1] for src in spinw_sources]
spinw_mds = build_examples(spinw_sources, joinpath("examples", "spinw"))
-contributed_mds = prepare_contributed()
+contributed_mds = isdraft ? [] : prepare_contributed()
# Build docs as HTML, including the `examples/name.md` markdown built above
@@ -114,7 +111,7 @@ Documenter.makedocs(;
"index.md",
"Examples" => [
example_mds...,
- "SpinW ports" => spinw_mds,
+ "SpinW tutorials" => spinw_mds,
"Contributed" => contributed_mds,
"Advanced" => [
"parallelism.md",
@@ -137,13 +134,9 @@ Documenter.makedocs(;
size_threshold_warn = 200*1024, # 200KB -- library.html gets quite large
size_threshold = 300*2024, # 300KB
),
- draft
+ draft = isdraft
)
-# Remove the (sparsely checked out) SunnyContributed git repo. This is only
-# necessary to enable repeated local builds. It has no effect on the CI.
-rm("contributed-tmp"; force=true, recursive=true)
-
# Attempt to push to gh-pages branch for deployment
Documenter.deploydocs(
repo = "github.com/SunnySuite/Sunny.jl.git",
diff --git a/docs/src/parallelism.md b/docs/src/parallelism.md
index 1fb4fbe8e..a186bdb7b 100644
--- a/docs/src/parallelism.md
+++ b/docs/src/parallelism.md
@@ -13,13 +13,13 @@ copied and pasted into your preferred Julia development environment.
## Review of the serial workflow
-The serial approach to calculating a structure factor, covered in [FeI₂ at
-Finite Temperature](@ref), involves thermalizing a spin `System` and then
-calling [`add_sample!`](@ref). `add_sample!` uses the state of the `System` as
-an initial condition for the calculation of a dynamical trajectory. The
-correlations of the trajectory are calculated and accumulated into a running
-average of the ``𝒮(𝐪,ω)``. This sequence is repeated to generate additional
-samples.
+The serial approach to calculating a structure factor, covered in the [FeI₂
+tutorial](@ref "4. Generalized spin dynamics of FeI₂ at finite *T*"), involves
+thermalizing a spin `System` and then calling [`add_sample!`](@ref).
+`add_sample!` uses the state of the `System` as an initial condition for the
+calculation of a dynamical trajectory. The correlations of the trajectory are
+calculated and accumulated into a running average of the ``𝒮(𝐪,ω)``. This
+sequence is repeated to generate additional samples.
To illustrate, we'll set up a a simple model: a spin-1 antiferromagnet on an FCC
crystal.
diff --git a/docs/src/versions.md b/docs/src/versions.md
index 27338a7e4..af3504ecf 100644
--- a/docs/src/versions.md
+++ b/docs/src/versions.md
@@ -229,7 +229,7 @@ An external field can be applied to a single site with
### Structure factor rewrite
The calculation of structure factors has been completely rewritten. For the new
-interface, see the [FeI₂ at Finite Temperature](@ref) page.
+interface see the documentation tutorials.
### Various
diff --git a/examples/fei2_tutorial.jl b/examples/01_LSWT_SU3_FeI2.jl
similarity index 95%
rename from examples/fei2_tutorial.jl
rename to examples/01_LSWT_SU3_FeI2.jl
index 2134b2345..87d8d4910 100644
--- a/examples/fei2_tutorial.jl
+++ b/examples/01_LSWT_SU3_FeI2.jl
@@ -1,4 +1,4 @@
-# # Case Study: FeI₂
+# # 1. Multi-flavor spin wave simulations of FeI₂ (Showcase)
#
# FeI₂ is an effective spin-1 material with strong single-ion anisotropy.
# Quadrupolar fluctuations give rise to a single-ion bound state that cannot be
@@ -36,8 +36,8 @@
# guide. Sunny requires Julia 1.9 or later.
#
# From the Julia prompt, load `Sunny`. For plotting, one can choose either
-# `GLMakie` (a pop-up window) or `WGLMakie` (inline plots for a Jupyter notebook
-# or VSCode).
+# `GLMakie` (a pop-up window) or `WGLMakie` (inline plots for a Jupyter
+# notebook).
using Sunny, GLMakie
@@ -69,7 +69,7 @@ FeI2 = Crystal(latvecs, positions; types)
cryst = subcrystal(FeI2, "Fe")
-# Observe that `cryst` retains the spacegroup symmetry of the full FeI₂ crystal.
+# Importantly, `cryst` retains the spacegroup symmetry of the full FeI₂ crystal.
# This information will be used, for example, to propagate exchange interactions
# between symmetry-equivalent bonds.
#
@@ -363,18 +363,12 @@ disp, is = dssf(swt, path);
# The full SU(_N_) coherent state dynamics, with appropriate quantum correction
# factors, can be useful to model finite temperature scattering data. In
# particular, it captures certain anharmonic effects due to thermal
-# fluctuations. This is the subject of our [FeI₂ at Finite Temperature](@ref)
-# tutorial.
+# fluctuations. See our [generalized spin dynamics tutorial](@ref "3.
+# Generalized spin dynamics of FeI₂ at finite *T*").
#
# The classical dynamics is also a good starting point to study non-equilibrium
# phenomena. Empirical noise and damping terms can be used to model [coupling to
# a thermal bath](https://arxiv.org/abs/2209.01265). This yields a Langevin
-# dynamics of SU(_N_) coherent states. Our [CP² Skyrmion Quench](@ref)
-# tutorial shows how this dynamics gives rise to the formation of novel
-# topological defects in a temperature quench.
-#
-# Relative to LSWT calculations, it can take much more time to estimate
-# $\mathcal{S}(𝐪,ω)$ intensities using classical dynamics simulation. See the
-# [SunnyTutorials
-# notebooks](https://nbviewer.org/github/SunnySuite/SunnyTutorials/tree/main/Tutorials/)
-# for examples of "production-scale" simulations.
+# dynamics of SU(_N_) coherent states. Our [dynamical SU(_N_) quench](@ref "5.
+# Dynamical quench into CP² skyrmion liquid") tutorial illustrates how a
+# temperature quench can give rise to novel liquid phase of CP² skyrmions.
diff --git a/examples/powder_averaging.jl b/examples/02_LSWT_CoRh2O4.jl
similarity index 83%
rename from examples/powder_averaging.jl
rename to examples/02_LSWT_CoRh2O4.jl
index 130b8e2b2..ffc77231b 100644
--- a/examples/powder_averaging.jl
+++ b/examples/02_LSWT_CoRh2O4.jl
@@ -1,9 +1,9 @@
-# # Powder Averaged CoRh₂O₄
+# # 2. Spin wave simulations of CoRh₂O₄
#
-# This tutorial illustrates the calculation of the powder-averaged structure
-# factor by performing an orientational average. We consider a simple model of
-# the diamond-cubic crystal CoRh₂O₄, with parameters extracted from [Ge et al.,
-# Phys. Rev. B 96, 064413](https://doi.org/10.1103/PhysRevB.96.064413).
+# This tutorial illustrates the conventional spin wave theory of dipoles. We
+# consider a simple model of the diamond-cubic crystal CoRh₂O₄, with parameters
+# extracted from [Ge et al., Phys. Rev. B 96,
+# 064413](https://doi.org/10.1103/PhysRevB.96.064413).
using Sunny, GLMakie
@@ -23,15 +23,14 @@ cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1")
view_crystal(cryst, 8.0)
-# Construct a [`System`](@ref) with an antiferromagnetic nearest neighbor
-# interaction `J`. Because the diamond crystal is bipartite, the ground state
-# will have unfrustrated Néel order. Selecting `latsize=(1,1,1)` is sufficient
-# because the ground state is periodic over each cubic unit cell. By passing an
-# explicit `seed`, the system's random number generator will give repeatable
-# results.
+# Construct a [`System`](@ref) with quantum spin ``S=3/2`` constrained to the
+# space of dipoles. Including an antiferromagnetic nearest neighbor interaction
+# `J` will favor Néel order. To optimize this magnetic structure, it is
+# sufficient to employ a magnetic lattice consisting of a single crystal unit
+# cell, `latsize=(1,1,1)`. Passing an explicit random number `seed` will ensure
+# repeatable results.
-latsize = (2, 2, 2)
-seed = 0
+latsize = (1, 1, 1)
S = 3/2
J = 7.5413*meV_per_K # (~ 0.65 meV)
sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole; seed=0)
@@ -83,7 +82,7 @@ formula = intensity_formula(swt, :perp; kernel, formfactors)
# [`intensities_broadened`](@ref) function collects intensities along this path
# for the given set of energy values.
-qpoints = [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.0], [0.0, 0.0, 0.0]]
+qpoints = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]
path, xticks = reciprocal_space_path(cryst, qpoints, 50)
energies = collect(0:0.01:6)
is = intensities_broadened(swt, path, energies, formula)
diff --git a/examples/03_LLD_CoRh2O4.jl b/examples/03_LLD_CoRh2O4.jl
new file mode 100644
index 000000000..74e026c49
--- /dev/null
+++ b/examples/03_LLD_CoRh2O4.jl
@@ -0,0 +1,191 @@
+# # 3. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*
+
+using Sunny, GLMakie, Statistics
+
+# ### System construction
+
+# Construct the system as in the previous [CoRh₂O₄ tutorial](@ref "2. Spin wave
+# simulations of CoRh₂O₄"). After optimization, the system will be in an
+# unfrustrated antiferromagnetic ground state.
+
+a = 8.5031 # (Å)
+latvecs = lattice_vectors(a, a, a, 90, 90, 90)
+cryst = Crystal(latvecs, [[0,0,0]], 227, setting="1")
+latsize = (2, 2, 2)
+S = 3/2
+J = 0.63 # (meV)
+sys = System(cryst, latsize, [SpinInfo(1; S, g=2)], :dipole; seed=0)
+set_exchange!(sys, J, Bond(1, 3, [0,0,0]))
+randomize_spins!(sys)
+minimize_energy!(sys)
+
+# Use [`resize_supercell`](@ref) to build a new system with a lattice of
+# 10×10×10 unit cells. The desired Néel order is retained.
+
+sys = resize_supercell(sys, (10, 10, 10))
+@assert energy_per_site(sys) ≈ -2J*S^2
+
+# Use the stochastic Landau-Lifshitz dynamics to thermalize system into
+# equilibrium at finite temperature. This is a [`Langevin`](@ref) equation,
+# which includes damping and noise terms. The strength of the noise term is
+# automatically fixed according to the damping time scale `λ` and the target
+# temperature, according to a fluctuation-dissipation theorem.
+
+Δt = 0.05/abs(J*S) # Time step
+λ = 0.1 # Dimensionless damping time-scale
+kT = 16 * meV_per_K # 16K, a temperature slightly below ordering
+langevin = Langevin(Δt; λ, kT);
+
+# Because the magnetic order has been initialized correctly, relatively few
+# additional Langevin time-steps are required to reach thermal equilibrium.
+
+energies = [energy_per_site(sys)]
+for _ in 1:1000
+ step!(sys, langevin)
+ push!(energies, energy_per_site(sys))
+end
+plot(energies, color=:blue, figure=(resolution=(600,300),), axis=(xlabel="Time steps", ylabel="Energy (meV)"))
+
+# Thermal fluctuations are apparent in the spin configuration.
+
+S_ref = sys.dipoles[1,1,1,1]
+plot_spins(sys; color=[s'*S_ref for s in sys.dipoles])
+
+# ### Instantaneous structure factor
+
+# To visualize the instantaneous (equal-time) structure factor, create an object
+# [`instant_correlations`](@ref) and use [`add_sample!`](@ref) to accumulated
+# data for the equilibrated spin configuration.
+
+sc = instant_correlations(sys)
+add_sample!(sc, sys) # Accumulate the newly sampled structure factor into `sf`
+
+# Collect 20 additional decorrelated samples. For each sample, about 100
+# Langevin time-steps is sufficient to collect approximately uncorrelated
+# statistics.
+
+for _ in 1:20
+ for _ in 1:100
+ step!(sys, langevin)
+ end
+ add_sample!(sc, sys)
+end
+
+# Define a slice of momentum space. Wavevectors are specified in reciprocal
+# lattice units (RLU). The notation `q1s in -10:0.1:10` indicates that the first
+# ``q``-component ranges from -10 to 10 in intervals of 0.1. That is, ``q``
+# spans over 20 Brillouin zones. To convert to absolute momentum units, each
+# component of ``q`` would need to be scaled by a reciprocal lattice vector.
+
+q1s = -10:0.1:10
+q2s = -10:0.1:10
+qs = [[q1, q2, 0.0] for q1 in q1s, q2 in q2s];
+
+# Plot the instantaneous structure factor for the given ``q``-slice. We employ
+# the appropriate [`FormFactor`](@ref) for Co2⁺. An [`intensity_formula`](@ref)
+# defines how dynamical correlations correspond to the observable structure
+# factor. The function [`instant_intensities_interpolated`](@ref) calculates
+# intensities at the target `qs` by interpolating over the data available at
+# discrete reciprocal-space lattice points.
+
+formfactors = [FormFactor("Co2")]
+instant_formula = intensity_formula(sc, :perp; formfactors)
+iq = instant_intensities_interpolated(sc, qs, instant_formula);
+
+# Plot the resulting intensity data ``I(𝐪)``. The color scale is clipped to 50%
+# of the maximum intensity.
+
+heatmap(q1s, q2s, iq;
+ colorrange = (0, maximum(iq)/2),
+ axis = (
+ xlabel="Momentum Transfer Qx (r.l.u)", xlabelsize=16,
+ ylabel="Momentum Transfer Qy (r.l.u)", ylabelsize=16,
+ aspect=true,
+ )
+)
+
+
+# ### Dynamical structure factor
+
+# To collect statistics for the dynamical structure factor intensities
+# ``I(𝐪,ω)`` at finite temperature, use [`dynamical_correlations`](@ref). Now,
+# each call to `add_sample!` will run a classical spin dynamics trajectory.
+# Longer-time trajectories will be required to achieve greater energy
+# resolution, as controlled by `nω`. Here, we pick a moderate number of
+# energies, `nω = 50`, which will make the simulation run quickly.
+
+ωmax = 6.0 # Maximum energy to resolve (meV)
+nω = 50 # Number of energies to resolve
+sc = dynamical_correlations(sys; Δt, nω, ωmax, process_trajectory=:symmetrize)
+
+# Each sample requires running a full dynamical trajectory to measure
+# correlations, so we here restrict to 5 samples.
+
+for _ in 1:5
+ for _ in 1:100
+ step!(sys, langevin)
+ end
+ add_sample!(sc, sys)
+end
+
+# Define points that define a piecewise-linear path through reciprocal space,
+# and a sampling density.
+
+points = [[3/4, 3/4, 0],
+ [ 0, 0, 0],
+ [ 0, 1/2, 1/2],
+ [1/2, 1, 0],
+ [ 0, 1, 0],
+ [1/4, 1, 1/4],
+ [ 0, 1, 0],
+ [ 0, -4, 0]]
+density = 50 # (Å)
+path, xticks = reciprocal_space_path(cryst, points, density);
+
+# Calculate ``I(𝐪, ω)`` intensities along this path with Lorentzian broadening
+# on the scale of 0.1 meV.
+
+formula = intensity_formula(sc, :perp; formfactors, kT=langevin.kT)
+η = 0.1
+iqw = intensities_interpolated(sc, path, formula)
+iqwc = broaden_energy(sc, iqw, (ω, ω₀) -> lorentzian(ω-ω₀, η));
+
+# Plot the intensity data on a clipped color scale
+
+ωs = available_energies(sc)
+heatmap(1:size(iqwc, 1), ωs, iqwc;
+ colorrange = (0, maximum(iqwc)/50),
+ axis = (;
+ xlabel="Momentum Transfer (r.l.u)",
+ ylabel="Energy Transfer (meV)",
+ xticks,
+ xticklabelrotation=π/5,
+ aspect = 1.4,
+ )
+)
+
+# ### Powder averaged intensity
+
+# Define spherical shells in reciprocal space via their radii, in absolute units
+# of 1/Å. For each shell, calculate and average the intensities at 100
+# ``𝐪``-points, sampled approximately uniformly.
+
+radii = 0:0.05:3.5 # (1/Å)
+output = zeros(Float64, length(radii), length(ωs))
+for (i, radius) in enumerate(radii)
+ pts = reciprocal_space_shell(sc.crystal, radius, 100)
+ is = intensities_interpolated(sc, pts, formula)
+ is = broaden_energy(sc, is, (ω,ω₀)->lorentzian(ω-ω₀, η))
+ output[i, :] = mean(is , dims=1)[1,:]
+end
+
+# Plot resulting powder-averaged structure factor
+
+heatmap(radii, ωs, output;
+ axis = (
+ xlabel="|Q| (Å⁻¹)",
+ ylabel="Energy Transfer (meV)",
+ aspect = 1.4,
+ ),
+ colorrange = (0, 20.0)
+)
diff --git a/examples/fei2_classical.jl b/examples/04_GSD_FeI2.jl
similarity index 97%
rename from examples/fei2_classical.jl
rename to examples/04_GSD_FeI2.jl
index 785ed9d22..d887d2580 100644
--- a/examples/fei2_classical.jl
+++ b/examples/04_GSD_FeI2.jl
@@ -1,11 +1,12 @@
-# # FeI₂ at Finite Temperature
+# # 4. Generalized spin dynamics of FeI₂ at finite *T*
using Sunny, LinearAlgebra, GLMakie
-# In our previous [Case Study: FeI₂](@ref), we used linear spin wave theory
-# (LSWT) to calculate the dynamical structure factor. Here, we perform a similar
-# calculation using classical spin dynamics. Because we are interested in the
-# coupled dynamics of spin dipoles and quadrupoles, we employ a [classical
+# In the [previous FeI₂ tutorial](@ref "1. Multi-flavor spin wave simulations of
+# FeI₂ (Showcase)"), we used multi-flavor spin wave theory to calculate the
+# dynamical structure factor. Here, we perform a similar calculation using
+# classical spin dynamics at finite temperature. Because we are interested in
+# the coupled dynamics of spin dipoles and quadrupoles, we employ a [classical
# dynamics of SU(3) coherent states](https://arxiv.org/abs/2209.01265) that
# generalizes the Landau-Lifshitz equation.
#
diff --git a/examples/ising2d.jl b/examples/05_MC_Ising.jl
similarity index 98%
rename from examples/ising2d.jl
rename to examples/05_MC_Ising.jl
index ee3d8068d..0dac4ecce 100644
--- a/examples/ising2d.jl
+++ b/examples/05_MC_Ising.jl
@@ -1,4 +1,4 @@
-# # Classical Ising model
+# # 5. Monte Carlo sampling of the Ising model
#
# This tutorial illustrates simulation of the classical 2D Ising model.
diff --git a/examples/out_of_equilibrium.jl b/examples/06_CP2_Skyrmions.jl
similarity index 97%
rename from examples/out_of_equilibrium.jl
rename to examples/06_CP2_Skyrmions.jl
index 94679480e..848f4f79b 100644
--- a/examples/out_of_equilibrium.jl
+++ b/examples/06_CP2_Skyrmions.jl
@@ -1,4 +1,4 @@
-# # CP² Skyrmion Quench
+# # 6. Dynamical quench into CP² skyrmion liquid
#
# This example demonstrates Sunny's ability to simulate the out-of-equilibrium
# dynamics of generalized spin systems. We will implement the model Hamiltonian
@@ -108,9 +108,9 @@ end
# function `plot_triangular_plaquettes` is not part of the core Sunny package,
# but rather something you could define yourself. We are using the definition in
# `plotting2d.jl` from the Sunny [`examples/extra`
-# directory](https://github.com/SunnySuite/Sunny.jl/tree/main/examples/extra).
+# directory](https://github.com/SunnySuite/Sunny.jl/tree/main/examples/extra/Plotting).
-include(pkgdir(Sunny, "examples", "extra", "plotting2d.jl"))
+include(pkgdir(Sunny, "examples", "extra", "Plotting", "plotting2d.jl"))
function sun_berry_curvature(z₁, z₂, z₃)
z₁, z₂, z₃ = normalize.((z₁, z₂, z₃))
diff --git a/examples/extra/PT_WHAM_ising2d.jl b/examples/extra/Advanced_MC/PT_WHAM_ising2d.jl
similarity index 100%
rename from examples/extra/PT_WHAM_ising2d.jl
rename to examples/extra/Advanced_MC/PT_WHAM_ising2d.jl
diff --git a/examples/extra/REWL_ising2d.jl b/examples/extra/Advanced_MC/REWL_ising2d.jl
similarity index 100%
rename from examples/extra/REWL_ising2d.jl
rename to examples/extra/Advanced_MC/REWL_ising2d.jl
diff --git a/examples/extra/WL_ising2d.jl b/examples/extra/Advanced_MC/WL_ising2d.jl
similarity index 100%
rename from examples/extra/WL_ising2d.jl
rename to examples/extra/Advanced_MC/WL_ising2d.jl
diff --git a/examples/extra/plotting2d.jl b/examples/extra/Plotting/plotting2d.jl
similarity index 100%
rename from examples/extra/plotting2d.jl
rename to examples/extra/Plotting/plotting2d.jl
diff --git a/examples/extra/binning_tutorial.jl b/examples/extra/binning_tutorial.jl
deleted file mode 100644
index af10d90dc..000000000
--- a/examples/extra/binning_tutorial.jl
+++ /dev/null
@@ -1,187 +0,0 @@
-# # Histogram Binning Tutorial
-
-using Sunny, GLMakie
-
-# This Tutorial demonstrates how to use Sunny's histogram binning
-# capabilities (via [`intensities_binned`](@ref)). This functionality allows
-# the simulation data produced by Sunny to be compared to experimental data produced by
-# Inelastic Neutron Scattering (INS) in an apples-to-apples fashion.
-# Experimental data can be loaded from a `MDHistoWorkspace` stored in a `.nxs` file
-# by the [Mantid software](https://www.mantidproject.org/) using [`load_nxs`](@ref).
-
-# For this example, we will examine the CTFD compound, which is crystallographically approximately
-# a square lattice.
-# We specify the crystal lattice structure of CTFD using the lattice parameters specified by
-#
-# W Wan et al 2020 J. Phys.: Condens. Matter 32 374007 DOI 10.1088/1361-648X/ab757a
-latvecs = lattice_vectors(8.113,8.119,12.45,90,100,90)
-positions = [[0,0,0]]
-types = ["Cu"]
-formfactors = [FormFactor("Cu2")]
-xtal = Crystal(latvecs,positions;types);
-
-# We will use a somewhat small periodic lattice size of 6x6x4 in order
-# to showcase the effect of a finite lattice size.
-latsize = (6,6,4);
-
-# In this system, the magnetic lattice is the same as the chemical lattice,
-# and there is a spin-1/2 dipole on each site.
-magxtal = xtal;
-valS = 1/2;
-sys = System(magxtal, latsize, [SpinInfo(1, S=valS, g=2)], :dipole; seed=1);
-
-# Quoted value of J = +6.19(2) meV (antiferromagnetic) between
-# nearest neighbors on the square lattice
-J = 6.19 # meV
-characteristic_energy_scale = abs(J * valS)
-set_exchange!(sys,J,Bond(1,1,[1,0,0]))
-set_exchange!(sys,J,Bond(1,1,[0,1,0]))
-
-# Thermalize the system using the Langevin intergator.
-# The timestep and the temperature are roughly based off the characteristic
-# energy scale of the problem.
-Δt = 0.05 / characteristic_energy_scale
-kT = 0.01 * characteristic_energy_scale
-langevin = Langevin(Δt; λ=0.1, kT=kT)
-randomize_spins!(sys);
-for i in 1:10_000 # Long enough to reach equilibrium
- step!(sys, langevin)
-end
-
-# The neutron spectrometer used in the experiment had an incident neutron energy of 36.25 meV.
-# Since this is the most amount of energy that can be deposited by the neutron into the sample, we
-# don't need to consider energies higher than this.
-ωmax = 36.25;
-
-# We choose the resolution in energy (specified by the number of nω modes resolved)
-# to be ≈20× better than the experimental resolution in order to demonstrate
-# the effect of over-resolving in energy.
-nω = 480;
-sc = dynamical_correlations(sys; Δt=Δt, nω=nω, ωmax=ωmax, process_trajectory=:symmetrize)
-add_sample!(sc, sys)
-
-# We re-sample from the thermal equilibrium distribution several times to increase our sample size
-nsamples = 3
-for _ in 1:nsamples
- for _ in 1:8000
- step!(sys, langevin)
- end
- add_sample!(sc, sys)
-end
-
-# Since the SU(N)NY crystal has only finitely many lattice sites, there are finitely
-# many ways for a neutron to scatter off of the sample.
-# We can visualize this discreteness by plotting each possible Qx and Qz, for example:
-
-params = unit_resolution_binning_parameters(sc)#hide
-Sunny.bin_rlu_as_absolute_units!(params,sc)#hide
-lower_aabb_q, upper_aabb_q = Sunny.binning_parameters_aabb(params)#hide
-lower_aabb_cell = floor.(Int64,lower_aabb_q .* latsize .+ 1)#hide
-upper_aabb_cell = ceil.(Int64,upper_aabb_q .* latsize .+ 1)#hide
-
-Qx = zeros(Float64,0)#hide
-Qz = zeros(Float64,0)#hide
-for cell in CartesianIndices(Tuple(((:).(lower_aabb_cell,upper_aabb_cell))))#hide
- q = (cell.I .- 1) ./ latsize # q is in R.L.U.#hide
- push!(Qx,q[1])#hide
- push!(Qz,q[3])#hide
-end#hide
-fig = Figure()#hide
-ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Qz [r.l.u.]")#hide
-## Compute some scattering vectors at and around the first BZ...
-scatter!(ax,Qx,Qz)
-fig#hide
-
-# One way to display the structure factor is to create a histogram with
-# one bin centered at each discrete scattering possibility using [`unit_resolution_binning_parameters`](@ref) to create a set of [`BinningParameters`](@ref).
-params = unit_resolution_binning_parameters(sc)
-
-# Since this is a 4D histogram,
-# it further has to be integrated over two of those directions in order to be displayed.
-# Here, we integrate over Qy and Energy using [`integrate_axes!`](@ref):
-integrate_axes!(params;axes = [2,4]) # Integrate over Qy (2) and E (4)
-
-# Now that we have parameterized the histogram, we can bin our data.
-# In addition to the [`BinningParameters`](@ref), an [`intensity_formula`](@ref) needs to be
-# provided to specify which dipole, temperature, and atomic form factor
-# corrections should be applied during the intensity calculation.
-formula = intensity_formula(sc, :perp; kT, formfactors)
-intensity,counts = intensities_binned(sc, params; formula)
-normalized_intensity = intensity ./ counts;
-
-# With the data binned, we can now plot it. The axes labels give the bin centers of each bin, as given by [`axes_bincenters`](@ref).
-function plot_data(params) #hide
-intensity,counts = intensities_binned(sc, params; formula)#hide
-normalized_intensity = intensity ./ counts;#hide
-bin_centers = axes_bincenters(params);
-
-fig = Figure()#hide
-ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Qz [r.l.u.]")#hide
-heatmap!(ax,bin_centers[1],bin_centers[3],normalized_intensity[:,1,:,1])
-scatter!(ax,Qx,Qz)
-xlims!(ax,params.binstart[1],params.binend[1])
-ylims!(ax,params.binstart[3],params.binend[3])
-return fig#hide
-
-end#hide
-
-plot_data(params)#hide
-
-# Note that some bins have no scattering vectors at all when the bin size is made too small:
-params = unit_resolution_binning_parameters(sc)#hide
-integrate_axes!(params;axes = [2,4])#hide
-params.binwidth[1] /= 1.2
-params.binwidth[3] /= 2.5
-plot_data(params)#hide
-
-# Conversely, making the bins bigger doesn't break anything, but loses resolution:
-params = unit_resolution_binning_parameters(sc)#hide
-integrate_axes!(params;axes = [2,4])#hide
-params.binwidth[1] *= 2
-params.binwidth[3] *= 2
-plot_data(params)#hide
-
-# Recall that while we under-resolved in Q by choosing a small lattice, we
-# over-resolved in energy:
-x = zeros(Float64,0)#hide
-y = zeros(Float64,0)#hide
-for omega = ωs(sc), qx = unique(Qx)#hide
- push!(x,qx)#hide
- push!(y,omega)#hide
-end#hide
-ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Energy [meV]")#hide
-scatter!(ax,x,y)
-
-# Let's make a new histogram which includes the energy axis.
-# The x-axis of the histogram will be a 1D cut from `Q = [0,0,0]` to `Q = [1,1,0]`.
-# See [`slice_2D_binning_parameters`](@ref).
-x_axis_bin_count = 10
-cut_width = 0.3
-params = slice_2D_binning_parameters(sc,[0,0,0],[1,1,0],x_axis_bin_count,cut_width)
-
-# There are no longer any scattering vectors exactly in the plane of the cut. Instead,
-# as described in the `BinningParameters` output above, the transverse Q
-# directions are integrated over, so slightly out of plane points are included.
-#
-# We plot the intensity on a log-scale to improve visibility.
-intensity,counts = intensities_binned(sc, params; formula)
-log_intensity = log10.(intensity ./ counts);
-bin_centers = axes_bincenters(params);#hide
-fig = Figure()#hide
-ax = Axis(fig[1,1];xlabel="Progress along cut [r.l.u]",ylabel="Energy [meV]")#hide
-heatmap!(ax,bin_centers[1],bin_centers[4],log_intensity[:,1,1,:])
-xlims!(ax,params.binstart[1],params.binend[1])#hide
-ylims!(ax,params.binstart[4],params.binend[4])#hide
-fig#hide
-
-# By reducing the number of energy bins to be closer to the number of bins on the x-axis, we can make the dispersion curve look nicer:
-params.binwidth[4] *= 20
-intensity,counts = intensities_binned(sc, params; formula)#hide
-log_intensity = log10.(intensity ./ counts);#hide
-bin_centers = axes_bincenters(params);#hide
-fig = Figure()#hide
-ax = Axis(fig[1,1];xlabel="Progress along cut [r.l.u]",ylabel="Energy [meV]")#hide
-heatmap!(ax,bin_centers[1],bin_centers[4],log_intensity[:,1,1,:])
-xlims!(ax,params.binstart[1],params.binend[1])#hide
-ylims!(ax,params.binstart[4],params.binend[4])#hide
-fig#hide
diff --git a/examples/extra/one_dim_chain.jl b/examples/extra/one_dim_chain.jl
deleted file mode 100644
index 8401df9b5..000000000
--- a/examples/extra/one_dim_chain.jl
+++ /dev/null
@@ -1,358 +0,0 @@
-# # Fitting model parameters in a 1D spin-1 ferromagnetic chain
-
-using Sunny, LinearAlgebra, GLMakie, Optim
-
-# In this Example, we consider a 1D chain of spin-1 sites.
-# The sites along the chain interact via a ferromagnetic nearest-neighbor
-# interaction ``J\sum_{\langle i,j\rangle} \mathbf{S}_i \cdot \mathbf{S}_j``, with ``J < 0``.
-# By default, the ground state would be ferromagnetic and highly degenerate, since the spins
-# can align in any direction.
-# An on-site interaction, ``D\sum_i (S^z_i)^2`` breaks this isotropy by making it easier for
-# the spins to align in the ``\pm z`` direction than in any other orientation.
-# Thus, the entire Hamiltonian is:
-#
-# ``\mathcal{H} = \overbrace{J\sum_{\langle i,j\rangle} \mathbf{S}_i \cdot \mathbf{S}_j}^{\text{Ferromagnetic}}\;\;\;\;\;\; \overbrace{-D\sum_i (S^z_i)^2}^{\text{Easy-axis single-ion anisotropy}}``
-#
-# The goal of this Example is to illustrate how to determine the parameters ``J`` and ``D``
-# from "experiment" data by fitting using Sunny's implementation of Linear Spin Wave Theory.
-# In our case, the "experiment" data will actually be simulation data produced
-# using Landau-Lifschitz dynamics.
-
-# ## Creating simulated "experiment" data using Landau-Lifschitz dynamics
-
-# Our simulated data will use ground truth values ``J_0 = -1\,\text{meV}`` and ``D_0 = 10\,\text{meV}`` with a lattice spacing ``a = 10`` angstrom.
-
-# We begin with a 1D chain of spin-1 sites along the ``x`` direction.
-
-## Establish geometry of the unit cell.
-## "P1" is required due to the rotational symmetry about the
-## x-axis being broken.
-chain_spacing = 10. # Angstrom
-latvecs = chain_spacing * I(3)
-one_dimensional_chain = Crystal(latvecs,[[0,0,0]],"P1")
-
-## Establish geometry of the whole chain.
-chain_length = 16 # Number of atoms
-latsize = (chain_length,1,1) # 1D chain is Nx1x1 lattice
-spin_one_chain = System(one_dimensional_chain, latsize, [SpinInfo(1,S=1,g=2)], :SUN)
-
-# Configure the nearest-neighbor interaction:
-
-## Scalar J indicates J*(Sᵢ⋅Sⱼ)
-J_groundtruth = -1.
-
-## Interaction is with the left and right next neighbor along the chain (x-direction)
-nearest_neighbor_right = Bond(1,1,(1,0,0))
-nearest_neighbor_left = Bond(1,1,(-1,0,0))
-
-set_exchange!(spin_one_chain,J_groundtruth,nearest_neighbor_right)
-set_exchange!(spin_one_chain,J_groundtruth,nearest_neighbor_left)
-
-# Configure the symmetry-breaking easy-axis term:
-D_groundtruth = 10.
-Sz = spin_matrices(1)[3]
-set_onsite_coupling!(spin_one_chain, -D_groundtruth*Sz^2, 1)
-
-# With the ground-truth hamiltonian in place, we use Sunny's classical dynamics
-# to generate ficticious experiment data at temperature `kT = 0.1`.
-
-Δt = 0.05/D_groundtruth
-λ = 0.1
-kT = 0.1
-langevin = Langevin(Δt; kT, λ);
-
-function viz_chain(sys;kwargs...)#hide
- ##ups = map(x -> abs2(x[1]), sys.coherents)[:];#hide
- ##zs = map(x -> abs2(x[2]), sys.coherents)[:];#hide
- ##downs = map(x -> abs2(x[3]), sys.coherents)[:];#hide
-###hide
- ##f = Figure()#hide
- ##ax = LScene(f[1,1];show_axis = false)#hide
- ##_ = Makie.cam3d!(ax.scene, projectiontype=Makie.Orthographic)#hide
-###hide
- ##linewidth = 5.#hide
- ##arrowsize = 10.#hide
- ##lengthscale = 15.#hide
- ##pts = [Point3f(Sunny.global_position(sys,site)) for site in eachsite(sys)][:]#hide
-###hide
- #### Ups#hide
- ##vecs = [Vec3f([0,0,1]) for site in eachsite(sys)][:]#hide
- ##cols = map(x -> (:blue,x), ups)#hide
- ##Makie.arrows!(ax, pts .+ 0.5 .* vecs, vecs;#hide
- ##linecolor = cols, arrowcolor = cols,#hide
- ##lengthscale, arrowsize, linewidth, kwargs...)#hide
-###hide
- #### Downs#hide
- ##vecs = [Vec3f([0,0,-1]) for site in eachsite(sys)][:]#hide
- ##cols = map(x -> (:red,x), downs)#hide
- ##Makie.arrows!(ax, pts .+ 0.5 .* vecs, vecs;#hide
- ##linecolor = cols, arrowcolor = cols,#hide
- ##lengthscale, arrowsize, linewidth, kwargs...)#hide
-###hide
- ##cols = map(x -> (:green,x), zs)#hide
- ##meshscatter!(ax,pts, markersize = 7., color = cols)#hide
- ##f#hide
- Sunny.Plotting.plot_coherents(sys;quantization_axis = [0,0,1],kwargs...)
-end#hide
-randomize_spins!(spin_one_chain)
-viz_chain(spin_one_chain)
-
-# In this plot, the z-axis has been used as the quantization axis for each site,
-# with the up/down arrows and circle representing the ``\pm \hbar`` and ``0\hbar`` spin projections
-# onto the z-axis respectively. The opacity of each object represents the probability (absolute value
-# squared), and the color represents the phase.
-# Since we are using classical dynamics to simulate the data, the phase will be mostly random.
-#
-# First, we thermalize the chain, and then take several samples
-# in order get reasonably good "experiment" data.
-
-nStep = 50_000#hide
-for _ in 1:nStep#hide
- step!(spin_one_chain, langevin)#hide
-end#hide
-## ... thermalize ...
-viz_chain(spin_one_chain)
-
-#
-
-sc = dynamical_correlations(spin_one_chain; Δt, nω = 80, ωmax = 20.);#hide
-
-for _ in 1:10_000#hide
- step!(spin_one_chain, langevin)#hide
-end#hide
-add_sample!(sc, spin_one_chain)#hide
-## ... some time later ...
-viz_chain(spin_one_chain)
-
-#
-
-for _ in 1:10_000#hide
- step!(spin_one_chain, langevin)#hide
-end#hide
-add_sample!(sc, spin_one_chain)#hide
-## ... some time later ...
-viz_chain(spin_one_chain)
-
-#
-
-for _ in 1:20#hide
- for _ in 1:10_000#hide
- step!(spin_one_chain, langevin)#hide
- end#hide
- add_sample!(sc, spin_one_chain)#hide
-end#hide
-## ... some time later ...
-viz_chain(spin_one_chain)
-
-# Now that we have collected several samples,
-
-sc
-
-# we are ready to generate the intensity data.
-# Since this is supposed to represent an experiment, the intensity data will go in a histogram:
-
-SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS = unit_resolution_binning_parameters(sc)
-
-# Here's what the experiment data looks like:
-
-formula = intensity_formula(sc,:perp;kT)
-is, counts = intensities_binned(sc,SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS,formula)
-
-SIMULATED_EXPERIMENT_DATA = (is ./ counts)[:,1,1,:]
-
-bcs = axes_bincenters(SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS)
-f = Figure()#hide
-ax = Axis(f[1,1])#hide
-heatmap!(ax,bcs[1],bcs[4],log10.(SIMULATED_EXPERIMENT_DATA))
-f#hide
-
-# ## Fitting to the experiment data
-
-# To fit this data, we first model the known aspects of the system in Sunny.
-# The first steps are the same whether we are simulating a known system or modelling an
-# unknown system:
-
-## Same as before
-chain_spacing = 10. # Angstrom
-latvecs = chain_spacing * I(3)
-one_dimensional_chain = Crystal(latvecs,[[0,0,0]],"P1")
-chain_length = 16 # Number of atoms
-latsize = (chain_length,1,1) # 1D chain is Nx1x1 lattice
-spin_one_chain = System(one_dimensional_chain, latsize, [SpinInfo(1,S=1,g=2)], :SUN)
-
-
-
-# Originally, the next step would have been to configure the hamiltonian by
-# specifying the ``J`` and ``D`` values. However, since these are unknowns,
-# we will avoid using them as long as possible, and instead proceed to set up the
-# bonds, spin operators, and `Langevin` integrator--none of which require the values:
-
-Δt = 0.05
-λ = 0.1
-kT = 0. # LSWT uses zero temperature
-langevin = Langevin(Δt; kT, λ);
-
-nearest_neighbor_right = Bond(1,1,(1,0,0))
-nearest_neighbor_left = Bond(1,1,(-1,0,0))
-
-# After this setup work is done *once*, we create a function `forward_problem(J_trial,D_trial)`
-# which will compute the Linear Spin Wave Theoretic spectrum at the trial values of the
-# ``J`` and ``D`` fitting parameters.
-# In other words, the part of the original calculation which depends on the fitting
-# parameters gets wrapped into a function:
-
-function forward_problem(J_trial, D_trial)
-
- ## Ensure there is no phase transition (or else LSWT will throw errors)
- J_trial = min(J_trial,0)
- D_trial = max(D_trial,0)
-
- ## Uses J_trial
- set_exchange!(spin_one_chain,J_trial,nearest_neighbor_right)
- set_exchange!(spin_one_chain,J_trial,nearest_neighbor_left)
-
- ## Uses D_trial
- set_onsite_coupling!(spin_one_chain, -D_trial*Sz^2, 1)
-
- ## Perform spin wave calculation, continued below...
-
-#+
-# Note that `forward_problem` refers to variables defined outside
-# of the scope of the function. This allows us to reuse those variables in each call to
-# `forward_problem`, without reconstructing them each time. In general, the more that is known
-# about the system you are modelling, the later in the code `function forward_problem(...)` can be
-# inserted, and the more setup work can be re-used.
-#
-# !!! tip "`forward_problem` is a closure"
-# In computer progrogramming parlance, `forward_problem` is said to 'capture' variables such
-# as `spin_one_chain` from the enviroment. Since the result of calling `forward_problem` depends
-# not only on `J_trial` and `D_trial`, but also on `spin_one_chain`, it's no longer a function
-# of only its arguments.
-#
-# Since `forward_problem` is not a closed system, but
-# `forward_problem + (captured variables)` _is_ a closed system, the latter is called
-# the 'closure' of the former.
-
-# ## Spin wave calculation
-
-# We can leverage our knowledge that the ground state should be ferromagnetic
-# to simplify the spin wave calculation. Since the ferrommagnetic unit cell is just one site,
-# the simplified system is extremely simple:
-
- ## ... perform spin wave calculation, continued from above.
- one_site_system = reshape_supercell(spin_one_chain,[1 0 0; 0 1 0; 0 0 1])
-
-#+
-# After restricting to a single site, it's best to re-thermalize the system
-# at zero temperature to ensure a good classical ground state for LSWT:
-
- langevin.kT = 0.
- nStep = 1_000
- for _ in 1:nStep
- step!(one_site_system, langevin)
- end
-
-#+
-# The spin wave intensity data must be placed in a histogram with the same parameters
-# as the experiment data, in order to ensure a good comparision.
-#
-# The kernel and `intensities_bin_centers` used here are temporary, until a better
-# binning method is written.
-
- swt = SpinWaveTheory(one_site_system)
- formula = intensity_formula(swt,:perp; kernel = lorentzian(0.5))
- params = SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS
- is_swt = Sunny.intensities_bin_centers(swt, params, formula)
-
- return is_swt[:,1,1,:]
-end # end of forward_problem
-
-# We can see the different possible results from LSWT by plotting the dispersion:
-
-function plot_forward(J,D)
- is_swt = forward_problem(J,D)
- bcs = axes_bincenters(SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS)
- heatmap(bcs[1],bcs[4],log10.(is_swt))
-end
-
-plot_forward(-1,10)
-
-#
-
-plot_forward(-6,2)
-
-#
-
-plot_forward(-0.01,15)
-
-# Now, we can easily define a least-squares loss function comparing the "experiment" data to the LSWT result:
-
-function get_loss(parameters)
- J,D = parameters
- is_swt = forward_problem(J,D)
- sqrt(sum(abs2.(SIMULATED_EXPERIMENT_DATA .- is_swt)))
-end
-
-# Sweeping the parameters over a range containing the true value reveals
-# that the loss is minimized near the true parameters (dot).
-# The minimum loss is not exactly at the ground truth parameters in this case.
-# Gradient descent (finite-differenced) can be used to find the actual minimizer:
-
-nJ = 30
-nD = 35
-loss_landscape = zeros(Float64,nJ,nD)
-Js = range(-2,0,length=nJ)
-Ds = range(8,12,length=nD)
-for (ij,J) in enumerate(Js)
- for (id,D) in enumerate(Ds)
- loss_landscape[ij,id] = get_loss([J,D])
- end
-end
-
-fig = Figure()
-ax = Axis(fig[1,1],xlabel = "J [meV]", ylabel = "D [meV]")
-contourf!(ax,Js,Ds,loss_landscape)
-
-x0 = [-2,9.5]
-opt_result = optimize(get_loss,x0,method=GradientDescent(alphaguess=1e-3),store_trace=true,extended_trace = true,time_limit=10.)
-lines!(ax,Point2f.(Optim.x_trace(opt_result)))
-scatter!(ax,-1,10)
-fig
-
-# The fit can be verified by plotting the LSWT band structure over top of the experiment data:
-
-bcs = axes_bincenters(SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS)
-f = Figure()#hide
-ax = Axis(f[1,1]; xlabel="Q [R.L.U.]", ylabel="Energy (meV)")#hide
-heatmap!(ax,bcs[1],bcs[4],log10.(SIMULATED_EXPERIMENT_DATA), colormap = :deepsea)
-f#hide
-
-
-J_trial, D_trial = opt_result.minimizer
-set_exchange!(spin_one_chain,J_trial,nearest_neighbor_right)#hide
-set_exchange!(spin_one_chain,J_trial,nearest_neighbor_left)#hide
-
-set_onsite_coupling!(spin_one_chain, -D_trial*Sz^2, 1)#hide
-one_site_system = reshape_supercell(spin_one_chain,[1 0 0; 0 1 0; 0 0 1])#hide
-
-langevin.kT = 0.#hide
-nStep = 1_000#hide
-for _ in 1:nStep#hide
- step!(one_site_system, langevin)#hide
-end#hide
-
-swt = SpinWaveTheory(one_site_system)#hide
-params = SIMULATED_EXPERIMENT_HISTOGRAM_PARAMS
-
-path = [[q,0,0] for q in bcs[1]]
-disp, intensity = intensities_bands(swt, path, intensity_formula(swt,:perp, kernel = delta_function_kernel))
-
-for i in axes(disp, 2)
- lines!(ax, bcs[1], disp[:,i]; color=intensity[:,i], colormap = :turbo,linewidth = 5,colorrange = (0.,1.))
-end
-Colorbar(f[1,2],colormap = :turbo, limits = (0.,1.))
-Colorbar(f[1,3],colormap = :deepsea, limits = (0.,1.))
-f
-
-
-
diff --git a/examples/longer_examples/CoRh2O4-tutorial.jl b/examples/longer_examples/CoRh2O4-tutorial.jl
deleted file mode 100644
index b0e58aec7..000000000
--- a/examples/longer_examples/CoRh2O4-tutorial.jl
+++ /dev/null
@@ -1,281 +0,0 @@
-# > ![](https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/assets/sunny_logo.jpg)
-# _This is a [tutorial](https://github.com/SunnySuite/SunnyTutorials/tree/main/tutorials)
-# for the [Sunny](https://github.com/SunnySuite/Sunny.jl/) package,
-# which enables dynamical simulations of ordered and thermally disordered spins with dipole
-# and higher order moments._
-#
-# ## Welcome to a Sunny Tutorial on the Diamond Lattice System CoRh₂O₄
-# **Script**: Diamond Lattice Finite Temperature Calculation
-# **Inspired by**: [Ge et al., Phys. Rev. B 96, 064413 (2017)](https://doi.org/10.1103/PhysRevB.96.064413)
-# **Authors**: Martin Mourigal, David Dahlbom
-# **Date**: August 21, 2023 (Sunny 0.5.0)
-# **Goal**: This script is to calculate the temperature dependence of the magnon excitations in the
-# spin-3/2 Heisenberg Diamond Antiferromagnet and compare to powder-averaged results obtained for
-# the compound CoRh₂O₄
-
-# ---
-# #### Loading Packages
-using Sunny, GLMakie, ProgressMeter, Statistics, Random, Brillouin
-Sunny.offline_viewers()
-cif_path = pkgdir(Sunny, "examples", "longer_examples", "CoRh2O4_#109301.cif");
-
-# #### Defining Custom functions
-
-# The function `quench!` randomizes the spins of a given `System`, fixes a
-# target temperature, and lets the system relax at this temperature for `nrelax`
-# integration steps.
-function quench!(sys, integrator; kTtarget, nrelax)
- randomize_spins!(sys);
- integrator.kT = kTtarget;
- prog = Progress(nrelax; dt=10.0, desc="Quenched and now relaxing: ", color=:green);
- for _ in 1:nrelax
- step!(sys, integrator)
- next!(prog)
- end
-end
-
-# `dwell!` takes a `System`, sets a target temperature, and has the system
-# dwell at this temperature for `ndwell` integration steps.
-function dwell!(sys, integrator; kTtarget, ndwell)
- integrator.kT = kTtarget;
- prog = Progress(ndwell; dt=10.0, desc="Dwelling: ", color=:green);
- for _ in 1:ndwell
- step!(sys, integrator)
- next!(prog)
- end
-end
-
-# `anneal!` takes a temperature schedule and cools the `System` through it,
-# with `ndwell` steps of the integrator at each temperature in the schedule.
-# Returns the energy at the end of the dwell for each scheduled temperature.
-function anneal!(sys, integrator; kTschedule, ndwell)
- nspins = prod(size(sys.dipoles));
- ensys = zeros(length(kTschedule))
- prog = Progress(ndwell*length(kTschedule); dt=10.0, desc="Annealing: ", color=:red);
- for (i, kT) in enumerate(kTschedule)
- integrator.kT = kT
- for _ in 1:ndwell
- step!(sys, integrator)
- next!(prog)
- end
- ensys[i] = energy(sys)
- end
- return ensys/nspins
-end
-
-# `sample_sf!` samples a structure factor, which may be either an instant or
-# dynamical structure factor. The integrator is run `ndecorr` times before each
-# one of the samples is taken.
-function sample_sf!(sf, sys, integrator; nsamples, ndecorr)
- prog = Progress(nsamples*ndecorr; dt=10.0, desc="Sampling SF: ", color=:red);
- for _ in 1:nsamples
- for _ in 1:ndecorr
- step!(sys, integrator)
- next!(prog)
- end
- add_sample!(sf, sys) # Accumulate the newly sampled structure factor into `sf`
- end
-end
-
-# `powder_average` powder averages a structure factor. Works for both instant
-# and dynamical structure factors. To prevent smearing, removes Bragg peaks
-# before introducing energy broadening. Bragg peaks are added back at ω=0 after
-# broadening.
-function powder_average(sc, rs, npts, formula; η=0.1)
- prog = Progress(length(rs); dt=10., desc="Powder Averaging: ", color=:blue)
- ωs = available_energies(sc)
- output = zeros(Float64, length(rs), length(ωs))
- for (i, r) in enumerate(rs)
- qs = reciprocal_space_shell(sc.crystal, r, npts)
- vals = intensities_interpolated(sc, qs, formula)
- bragg_idxs = findall(x -> x > maximum(vals)*0.9, vals)
- bragg_vals = vals[bragg_idxs]
- vals[bragg_idxs] .= 0
- vals = broaden_energy(sc, vals, (ω,ω₀)->lorentzian(ω-ω₀, η))
- vals[bragg_idxs] .= bragg_vals
- output[i,:] .= mean(vals, dims=1)[1,:]
- next!(prog)
- end
- return output
-end
-
-# ---
-# ### System Definition for CoRh₂O₄
-
-# Define the crystal structure of CoRh₂O₄ in the conventional cell
-xtal = Crystal(cif_path; symprec=1e-4)
-magxtal = subcrystal(xtal,"Co1")
-view_crystal(magxtal,6.0)
-print_symmetry_table(magxtal, 4.0)
-
-# Assign local Hilbert space
-S = 3/2
-lhs = [SpinInfo(1, S=S, g=2)]
-formfactors = [FormFactor("Co2")];
-
-# Create `System` and randomize it
-sunmode = :dipole
-latsize = (10,10,10)
-sys = System(magxtal, latsize, lhs, sunmode; seed=1)
-randomize_spins!(sys)
-plot_spins(sys)
-
-# Define Exchange Interactions
-scaleJ = 0.63
-valJ1 = 1.00*scaleJ
-set_exchange!(sys, valJ1, Bond(1, 3, [0, 0, 0]));
-
-# ---
-# ### System thermalization to an ordered, yet finite temperature, state
-
-# Define Langevin Integrator and Initialize it
-Δt0 = 0.05/abs(scaleJ*S); ## Time steps in Langevin
-λ0 = 0.1; ## Langevin damping, usually 0.05 or 0.1 is good.
-kT0 = 0.01*abs(scaleJ*S); ## Initialize at some temperature
-integrator = Langevin(Δt0; λ=λ0, kT=kT0);
-
-# Thermalization
-# Option 1: Quench the system from infinite temperature to a target temperature.
-# Note: this may lead to a poorly thermalized sample
-quench!(sys, integrator; kTtarget=kT0, nrelax=10000);
-
-# Option 2: Anneal (according to a temperature schedule) than dwell once reach base
-# Note: starting from very high temperature here
-
-## kTs = [abs(scaleJ)*valS*100 * 0.9^k for k in 0:100]
-## anneal!(sys,integrator;kTschedule=kTs,ndwell=500)
-## dwell!(sys,integrator;kTtarget=kTs[end],ndwell=2000)
-
-# Plot the resulting spin system to check ordering in real space
-plot_spins(sys)
-
-# ---
-# ### Calculation of Neutron Scattering Responses
-
-
-# #### Fourier transformed instantaneous two-point correlation functions
-
-# Calculate the Instantaneous/Equal-time Structure Factor
-eqsf = instant_correlations(sys)
-
-# If desired, add additional samples by decorrelating and then re-calculating the eqsf
-nsamples = 1
-ndecorr = 1000
-@time sample_sf!(eqsf, sys, integrator; nsamples=nsamples, ndecorr=ndecorr);
-
-# Project onto a constant Q-Slice in momentum space
-nQpts = 200
-Qxpts = range(-10.0, 10.0, length=nQpts)
-Qypts = range(-10.0, 10.0, length=nQpts)
-qz = 1.0
-Qpts = [[qx, qy, qz] for qx in Qxpts, qy in Qypts]
-instant_formula = intensity_formula(eqsf, :perp; formfactors)
-iq = instant_intensities_interpolated(eqsf, Qpts, instant_formula)
-
-# Plot the resulting I(Q)
-heatmap(Qxpts, Qypts, iq;
- colorrange = (0, maximum(iq)/20),
- axis = (
- xlabel="Momentum Transfer Qx (r.l.u)", xlabelsize=16,
- ylabel="Momentum Transfer Qy (r.l.u)", ylabelsize=16,
- aspect=true,
- )
-)
-
-
-# #### Dynamical and energy-integrated two-point correlation functions
-
-# Calculate the Time Traces and Fourier Transform: Dynamical Structure Factor (first sample)
-ωmax = 6.0 # Maximum energy to resolve
-nω = 100 # Number of energies to resolve
-sc = dynamical_correlations(sys; Δt=Δt0, nω=nω, ωmax=ωmax, process_trajectory=:symmetrize)
-@time add_sample!(sc, sys) # Add a sample trajectory
-
-# If desired, add additional decorrelated samples.
-nsamples = 10
-ndecorr = 1000
-@time sample_sf!(sc, sys, integrator; nsamples=nsamples, ndecorr=ndecorr);
-
-# Can use the Brillouin package for help on determining high symmetry points
-kp = irrfbz_path(227,[[1,0,0], [0,1,0], [0,0,1]])
-kpc = cartesianize(kp)
-
-# Project onto a constant QE-Slice in momentum-energy space.
-densQpts = 50
-symQpts = [[0.75, 0.75, 0.00], # List of wave vectors that define a path
- [0.00, 0.00, 0.00],
- [0.50, 0.50, 0.50],
- [0.50, 1.00, 0.00],
- [0.00, 1.00, 0.00],
- [0.25, 1.00, 0.25],
- [0.00, 1.00, 0.00],
- [0.00,-4.00, 0.00]]
-(Qpts, xticks) = reciprocal_space_path(magxtal, symQpts, densQpts)
-formula = intensity_formula(sc, :perp; formfactors, kT=integrator.kT)
-iqw = intensities_interpolated(sc, Qpts, formula);
-
-# If desired, broaden the sc in energy
-η = 0.1 ## Lorentzian energy broadening parameter
-iqwc = broaden_energy(sc, iqw, (ω, ω₀) -> lorentzian(ω-ω₀, η));
-
-# If desired, calculated the energy-integrated structure factor
-iqt = instant_intensities_interpolated(sc, Qpts, formula);
-
-# Plot the resulting I(Q,W)
-ωs = available_energies(sc)
-heatmap(1:size(iqwc, 1), ωs, iqwc;
- colorrange = (0, maximum(iqwc)/20000.0),
- axis = (;
- xlabel="Momentum Transfer (r.l.u)",
- ylabel="Energy Transfer (meV)",
- xticks,
- xticklabelrotation=π/5,
- aspect = 1.4,
- )
-)
-
-# Projection into a powder-averaged neutron scattering intensity
-Qmax = 3.5
-nQpts = 100
-Qpow = range(0, Qmax, nQpts)
-npoints = 100
-pqw = powder_average(sc, Qpow, npoints, formula; η);
-
-# Plot resulting Ipow(Q,W)
-heatmap(Qpow, ωs, pqw;
- axis = (
- xlabel="|Q| (Å⁻¹)",
- ylabel="Energy Transfer (meV)",
- aspect = 1.4,
- ),
- colorrange = (0, 40.0)
-)
-
-# ---
-# ### Calculation of temperature-dependent powder average spectrum
-
-# Define a temperature schedule
-kTs = [60 40 25 20 15 12 10 4] * Sunny.meV_per_K
-pqw_res = []
-for kT in kTs
- dwell!(sys, integrator; kTtarget=kT, ndwell=1000);
- sc_loc = dynamical_correlations(sys; Δt=Δt0, nω, ωmax, process_trajectory=:symmetrize);
- add_sample!(sc_loc, sys)
- formula = intensity_formula(sc, :perp; formfactors, kT)
- push!(pqw_res, powder_average(sc_loc, Qpow, npoints, formula; η))
-end
-
-# Plot the resulting Ipow(Q,W) as a function of temperature,
-# to compare with Fig.6 of https://arxiv.org/abs/1706.05881
-fig = Figure(; resolution=(1200,600))
-for i in 1:8
- r, c = fldmod1(i, 4)
- ax = Axis(fig[r, c];
- title = "kT = "*string(round(kTs[9-i], digits=3))*" (meV)",
- xlabel = r == 2 ? "|Q| (Å⁻¹)" : "",
- ylabel = c == 1 ? "Energy Transfer (meV)" : "",
- aspect = 1.4,
- )
- heatmap!(ax, Qpow, ωs, pqw_res[9-i]; colorrange = (0, 20.0))
-end
-fig
\ No newline at end of file
diff --git a/examples/longer_examples/CoRh2O4_#109301.cif b/examples/longer_examples/CoRh2O4_#109301.cif
deleted file mode 100644
index a29e48e83..000000000
--- a/examples/longer_examples/CoRh2O4_#109301.cif
+++ /dev/null
@@ -1,252 +0,0 @@
-
-#(C) 2022 by FIZ Karlsruhe - Leibniz Institute for Information Infrastructure. All rights reserved.
-data_109301-ICSD
-_database_code_ICSD 109301
-_audit_creation_date 2007-04-01
-_chemical_name_common 'Cobalt dirhodium(III) oxide'
-_chemical_formula_structural 'Co Rh2 O4'
-_chemical_formula_sum 'Co1 O4 Rh2'
-_chemical_name_structure_type Spinel#MgAl2O4
-_exptl_crystal_density_diffrn 7.04
-_citation_title 'Rhodites spinelles'
-loop_
-_citation_id
-_citation_journal_full
-_citation_year
-_citation_journal_volume
-_citation_page_first
-_citation_page_last
-_citation_journal_id_ASTM
-primary 'Comptes Rendus Hebdomadaires des Seances de l'Academie des Sciences'
-1959 249 726 728 COREAF
-loop_
-_citation_author_citation_id
-_citation_author_name
-primary 'Bertaut, F.'
-primary 'Forrat, F.'
-primary 'Dulac, J.F.'
-_cell_length_a 8.495
-_cell_length_b 8.495
-_cell_length_c 8.495
-_cell_angle_alpha 90.
-_cell_angle_beta 90.
-_cell_angle_gamma 90.
-_cell_volume 613.04
-_cell_formula_units_Z 8
-_space_group_name_H-M_alt 'F d -3 m S'
-_space_group_IT_number 227
-loop_
-_space_group_symop_id
-_space_group_symop_operation_xyz
-1 'z+1/4, y+1/4, -x+1/4'
-2 'y+1/4, x+1/4, -z+1/4'
-3 'x+1/4, z+1/4, -y+1/4'
-4 'z+1/4, x+1/4, -y+1/4'
-5 'y+1/4, z+1/4, -x+1/4'
-6 'x+1/4, y+1/4, -z+1/4'
-7 'z+1/4, -y+1/4, x+1/4'
-8 'y+1/4, -x+1/4, z+1/4'
-9 'x+1/4, -z+1/4, y+1/4'
-10 'z+1/4, -x+1/4, y+1/4'
-11 'y+1/4, -z+1/4, x+1/4'
-12 'x+1/4, -y+1/4, z+1/4'
-13 '-z+1/4, y+1/4, x+1/4'
-14 '-y+1/4, x+1/4, z+1/4'
-15 '-x+1/4, z+1/4, y+1/4'
-16 '-z+1/4, x+1/4, y+1/4'
-17 '-y+1/4, z+1/4, x+1/4'
-18 '-x+1/4, y+1/4, z+1/4'
-19 '-z+1/4, -y+1/4, -x+1/4'
-20 '-y+1/4, -x+1/4, -z+1/4'
-21 '-x+1/4, -z+1/4, -y+1/4'
-22 '-z+1/4, -x+1/4, -y+1/4'
-23 '-y+1/4, -z+1/4, -x+1/4'
-24 '-x+1/4, -y+1/4, -z+1/4'
-25 '-z, -y, x'
-26 '-y, -x, z'
-27 '-x, -z, y'
-28 '-z, -x, y'
-29 '-y, -z, x'
-30 '-x, -y, z'
-31 '-z, y, -x'
-32 '-y, x, -z'
-33 '-x, z, -y'
-34 '-z, x, -y'
-35 '-y, z, -x'
-36 '-x, y, -z'
-37 'z, -y, -x'
-38 'y, -x, -z'
-39 'x, -z, -y'
-40 'z, -x, -y'
-41 'y, -z, -x'
-42 'x, -y, -z'
-43 'z, y, x'
-44 'y, x, z'
-45 'x, z, y'
-46 'z, x, y'
-47 'y, z, x'
-48 'x, y, z'
-49 'z+1/4, y+3/4, -x+3/4'
-50 'y+1/4, x+3/4, -z+3/4'
-51 'x+1/4, z+3/4, -y+3/4'
-52 'z+1/4, x+3/4, -y+3/4'
-53 'y+1/4, z+3/4, -x+3/4'
-54 'x+1/4, y+3/4, -z+3/4'
-55 'z+1/4, -y+3/4, x+3/4'
-56 'y+1/4, -x+3/4, z+3/4'
-57 'x+1/4, -z+3/4, y+3/4'
-58 'z+1/4, -x+3/4, y+3/4'
-59 'y+1/4, -z+3/4, x+3/4'
-60 'x+1/4, -y+3/4, z+3/4'
-61 '-z+1/4, y+3/4, x+3/4'
-62 '-y+1/4, x+3/4, z+3/4'
-63 '-x+1/4, z+3/4, y+3/4'
-64 '-z+1/4, x+3/4, y+3/4'
-65 '-y+1/4, z+3/4, x+3/4'
-66 '-x+1/4, y+3/4, z+3/4'
-67 '-z+1/4, -y+3/4, -x+3/4'
-68 '-y+1/4, -x+3/4, -z+3/4'
-69 '-x+1/4, -z+3/4, -y+3/4'
-70 '-z+1/4, -x+3/4, -y+3/4'
-71 '-y+1/4, -z+3/4, -x+3/4'
-72 '-x+1/4, -y+3/4, -z+3/4'
-73 '-z, -y+1/2, x+1/2'
-74 '-y, -x+1/2, z+1/2'
-75 '-x, -z+1/2, y+1/2'
-76 '-z, -x+1/2, y+1/2'
-77 '-y, -z+1/2, x+1/2'
-78 '-x, -y+1/2, z+1/2'
-79 '-z, y+1/2, -x+1/2'
-80 '-y, x+1/2, -z+1/2'
-81 '-x, z+1/2, -y+1/2'
-82 '-z, x+1/2, -y+1/2'
-83 '-y, z+1/2, -x+1/2'
-84 '-x, y+1/2, -z+1/2'
-85 'z, -y+1/2, -x+1/2'
-86 'y, -x+1/2, -z+1/2'
-87 'x, -z+1/2, -y+1/2'
-88 'z, -x+1/2, -y+1/2'
-89 'y, -z+1/2, -x+1/2'
-90 'x, -y+1/2, -z+1/2'
-91 'z, y+1/2, x+1/2'
-92 'y, x+1/2, z+1/2'
-93 'x, z+1/2, y+1/2'
-94 'z, x+1/2, y+1/2'
-95 'y, z+1/2, x+1/2'
-96 'x, y+1/2, z+1/2'
-97 'z+3/4, y+1/4, -x+3/4'
-98 'y+3/4, x+1/4, -z+3/4'
-99 'x+3/4, z+1/4, -y+3/4'
-100 'z+3/4, x+1/4, -y+3/4'
-101 'y+3/4, z+1/4, -x+3/4'
-102 'x+3/4, y+1/4, -z+3/4'
-103 'z+3/4, -y+1/4, x+3/4'
-104 'y+3/4, -x+1/4, z+3/4'
-105 'x+3/4, -z+1/4, y+3/4'
-106 'z+3/4, -x+1/4, y+3/4'
-107 'y+3/4, -z+1/4, x+3/4'
-108 'x+3/4, -y+1/4, z+3/4'
-109 '-z+3/4, y+1/4, x+3/4'
-110 '-y+3/4, x+1/4, z+3/4'
-111 '-x+3/4, z+1/4, y+3/4'
-112 '-z+3/4, x+1/4, y+3/4'
-113 '-y+3/4, z+1/4, x+3/4'
-114 '-x+3/4, y+1/4, z+3/4'
-115 '-z+3/4, -y+1/4, -x+3/4'
-116 '-y+3/4, -x+1/4, -z+3/4'
-117 '-x+3/4, -z+1/4, -y+3/4'
-118 '-z+3/4, -x+1/4, -y+3/4'
-119 '-y+3/4, -z+1/4, -x+3/4'
-120 '-x+3/4, -y+1/4, -z+3/4'
-121 '-z+1/2, -y, x+1/2'
-122 '-y+1/2, -x, z+1/2'
-123 '-x+1/2, -z, y+1/2'
-124 '-z+1/2, -x, y+1/2'
-125 '-y+1/2, -z, x+1/2'
-126 '-x+1/2, -y, z+1/2'
-127 '-z+1/2, y, -x+1/2'
-128 '-y+1/2, x, -z+1/2'
-129 '-x+1/2, z, -y+1/2'
-130 '-z+1/2, x, -y+1/2'
-131 '-y+1/2, z, -x+1/2'
-132 '-x+1/2, y, -z+1/2'
-133 'z+1/2, -y, -x+1/2'
-134 'y+1/2, -x, -z+1/2'
-135 'x+1/2, -z, -y+1/2'
-136 'z+1/2, -x, -y+1/2'
-137 'y+1/2, -z, -x+1/2'
-138 'x+1/2, -y, -z+1/2'
-139 'z+1/2, y, x+1/2'
-140 'y+1/2, x, z+1/2'
-141 'x+1/2, z, y+1/2'
-142 'z+1/2, x, y+1/2'
-143 'y+1/2, z, x+1/2'
-144 'x+1/2, y, z+1/2'
-145 'z+3/4, y+3/4, -x+1/4'
-146 'y+3/4, x+3/4, -z+1/4'
-147 'x+3/4, z+3/4, -y+1/4'
-148 'z+3/4, x+3/4, -y+1/4'
-149 'y+3/4, z+3/4, -x+1/4'
-150 'x+3/4, y+3/4, -z+1/4'
-151 'z+3/4, -y+3/4, x+1/4'
-152 'y+3/4, -x+3/4, z+1/4'
-153 'x+3/4, -z+3/4, y+1/4'
-154 'z+3/4, -x+3/4, y+1/4'
-155 'y+3/4, -z+3/4, x+1/4'
-156 'x+3/4, -y+3/4, z+1/4'
-157 '-z+3/4, y+3/4, x+1/4'
-158 '-y+3/4, x+3/4, z+1/4'
-159 '-x+3/4, z+3/4, y+1/4'
-160 '-z+3/4, x+3/4, y+1/4'
-161 '-y+3/4, z+3/4, x+1/4'
-162 '-x+3/4, y+3/4, z+1/4'
-163 '-z+3/4, -y+3/4, -x+1/4'
-164 '-y+3/4, -x+3/4, -z+1/4'
-165 '-x+3/4, -z+3/4, -y+1/4'
-166 '-z+3/4, -x+3/4, -y+1/4'
-167 '-y+3/4, -z+3/4, -x+1/4'
-168 '-x+3/4, -y+3/4, -z+1/4'
-169 '-z+1/2, -y+1/2, x'
-170 '-y+1/2, -x+1/2, z'
-171 '-x+1/2, -z+1/2, y'
-172 '-z+1/2, -x+1/2, y'
-173 '-y+1/2, -z+1/2, x'
-174 '-x+1/2, -y+1/2, z'
-175 '-z+1/2, y+1/2, -x'
-176 '-y+1/2, x+1/2, -z'
-177 '-x+1/2, z+1/2, -y'
-178 '-z+1/2, x+1/2, -y'
-179 '-y+1/2, z+1/2, -x'
-180 '-x+1/2, y+1/2, -z'
-181 'z+1/2, -y+1/2, -x'
-182 'y+1/2, -x+1/2, -z'
-183 'x+1/2, -z+1/2, -y'
-184 'z+1/2, -x+1/2, -y'
-185 'y+1/2, -z+1/2, -x'
-186 'x+1/2, -y+1/2, -z'
-187 'z+1/2, y+1/2, x'
-188 'y+1/2, x+1/2, z'
-189 'x+1/2, z+1/2, y'
-190 'z+1/2, x+1/2, y'
-191 'y+1/2, z+1/2, x'
-192 'x+1/2, y+1/2, z'
-loop_
-_atom_type_symbol
-_atom_type_oxidation_number
-Co2+ 2
-Rh3+ 3
-O2- -2
-loop_
-_atom_site_label
-_atom_site_type_symbol
-_atom_site_symmetry_multiplicity
-_atom_site_Wyckoff_symbol
-_atom_site_fract_x
-_atom_site_fract_y
-_atom_site_fract_z
-_atom_site_B_iso_or_equiv
-_atom_site_occupancy
-Co1 Co2+ 8 a 0 0 0 . 1.
-Rh1 Rh3+ 16 d 0.625 0.625 0.625 . 1.
-O1 O2- 32 e 0.389 0.389 0.389 . 1.
-#End of TTdata_109301-ICSD
\ No newline at end of file
diff --git a/examples/longer_examples/MgCr2O4-tutorial.jl b/examples/longer_examples/MgCr2O4-tutorial.jl
deleted file mode 100644
index 86490238d..000000000
--- a/examples/longer_examples/MgCr2O4-tutorial.jl
+++ /dev/null
@@ -1,349 +0,0 @@
-# > ![](https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/assets/sunny_logo.jpg)
-# _This is a
-# [tutorial](https://github.com/SunnySuite/SunnyTutorials/)
-# for the [Sunny](https://github.com/SunnySuite/Sunny.jl/) package, which
-# enables dynamical simulations of ordered and thermally disordered spins with
-# dipole and higher order moments._
-#
-# # Spin Dynamics of the Heisenberg pyrochlore antiferromagnet and applications to MgCr2O4
-# **Author**: Martin Mourigal
-# **Date**: September 9, 2022 (Updated by David Dahlbom on August 21, 2023 using Sunny 0.5.0)
-#
-#
-# In this tutorial, we will walk through a first example in Sunny and calculate
-# the spin dynamical properties of the Heisenberg pyrochlore antiferromagnet and
-# apply this knowledge to MgCr2O4 and ZnCr2O4, which are known to approximate
-# this model. Relevant publications include:
-#
-# [1] P. H. Conlon and J. T. Chalker, Phys. Rev. Lett. 102, 237206 (2009)
-# (https://doi.org/10.1103/PhysRevLett.102.237206)
-#
-# [2] P. H. Conlon and J. T. Chalker, Phys. Rev. B 81, 224413 (2010)
-# (https://doi.org/10.1103/PhysRevB.81.224413)
-#
-# [3] X. Bai, J. A. M. Paddison, _et al._ Phys. Rev. Lett. 122, 097201 (2019)
-# (https://doi.org/10.1103/PhysRevLett.122.097201)
-#
-# ## Setting up Julia
-# To run the examples in the tutorial, you will need a working installation of
-# the Julia programming language and the Sunny package. Some useful references
-# for getting started are:
-#
-# https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia-and-Sunny
-#
-# https://sunnysuite.github.io/Sunny.jl/dev/
-#
-# We will begin by loading the relevant packages.
-
-using Sunny # The main package
-using GLMakie # Plotting package
-
-# ## Setting up the crystal structure
-#
-# Before specifying the interactions of our system, we first must set up the
-# crystal. We will begin by specifying the pyrochlore lattice (illustrated
-# below) in the manner that is typical of theorists.
-
-# ![](https://www.oist.jp/sites/default/files/photos/20200110-diagram-of-pyrochlore-lattice.jpg)
-#
-# _Picture Credits: Theory of Quantum Matter Unit, OIST_
-
-# ### "Theorist" Method
-
-# In this approach, we directly define the lattice vectors and specify the position
-# of each atom in fractional coordinates.
-latvecs = lattice_vectors(8.3342, 8.3342, 8.3342, 90, 90, 90)
-positions = [[0.875, 0.625, 0.375],
- [0.625, 0.125, 0.625],
- [0.875, 0.875, 0.125],
- [0.625, 0.875, 0.375],
- [0.875, 0.125, 0.875],
- [0.625, 0.625, 0.125],
- [0.875, 0.375, 0.625],
- [0.625, 0.375, 0.875],
- [0.375, 0.625, 0.875],
- [0.125, 0.125, 0.125],
- [0.375, 0.875, 0.625],
- [0.125, 0.875, 0.875],
- [0.375, 0.125, 0.375],
- [0.125, 0.625, 0.625],
- [0.375, 0.375, 0.125],
- [0.125, 0.375, 0.375]];
-types = ["B" for _ in positions]
-xtal_pyro = Crystal(latvecs, positions; types) # We will call this crystal the Theoretical Pyrochlore
-
-# To examine the result interactively, we can call `view_crystal`.
-
-view_crystal(xtal_pyro, 3.2)
-
-# ### "Experimentalist" Method #1 (Incorrect)
-# A real crystal is more complicated than this, however, and we will now
-# construct the system using the actual CIF file of MgCr2O4 from ICSD. This can
-# be done by copying over the data from a CIF file by hand, but this can be
-# dangerous, as shown below.
-
-latvecs = lattice_vectors(8.3342, 8.3342, 8.3342, 90, 90, 90)
-positions = [[0.1250, 0.1250, 0.1250],
- [0.5000, 0.5000, 0.5000],
- [0.2607, 0.2607, 0.2607]]
-types = ["Mg","Cr","O"]
-xtal_mgcro_1 = Crystal(latvecs, positions; types)
-
-# Sunny returned a valid crystal, but it did get right space group for MgCr2O4.
-# This can be fixed by modifying the input to include the space group and the
-# setting.
-#
-# ### "Experimentalist" Method #2 (Correct)
-
-# As above, we will define the crystal structure of MgCr2O4 by copying the info
-# from a CIF file, but we will also specify the space group and setting
-# explicitly.
-
-latvecs = lattice_vectors(8.3342, 8.3342, 8.3342, 90, 90, 90)
-positions = [[0.1250, 0.1250, 0.1250],
- [0.5000, 0.5000, 0.5000],
- [0.2607, 0.2607, 0.2607]]
-types = ["Mg","Cr","O"]
-spacegroup = 227 # Space Group Number
-setting = "2" # Space Group setting
-xtal_mgcro_2 = Crystal(latvecs, positions, spacegroup; types, setting)
-
-# This result is correct, but at this point we might as well import the CIF file
-# directly, which we now proceed to do.
-
-# ### "Experimentalist" Method #3 (Correct -- if your CIF file is)
-
-# To import a CIF file, simply give the path to `Crystal`. One may optionally
-# specify a precision parameter to apply to the symmetry analysis.
-
-cif = pkgdir(Sunny, "examples", "longer_examples", "MgCr2O4_160953_2009.cif")
-xtal_mgcro_3 = Crystal(cif; symprec=0.001)
-
-# Finally, we wish to restrict attention to the magnetic atoms in the unit cell
-# while maintaining symmetry information for the full crystal, which is required
-# to determine the correct exchange and g-factor anisotropies. This can be
-# achieved with the `subcrystal` function.
-
-xtal_mgcro = subcrystal(xtal_mgcro_2,"Cr")
-
-# ## Making a `System` and assigning interactions
-# ### Making a `System`
-# Before assigning any interactions, we first have to set up a `System` using
-# our `Crystal`.
-
-dims = (6, 6, 6) # Supercell dimensions
-spininfos = [SpinInfo(1, S=3/2, g=2)] # Specify spin information, note that all sites are symmetry equivalent
-sys_pyro = System(xtal_pyro, dims, spininfos, :dipole) # Make a system in dipole (Landau-Lifshitz) mode on pyrochlore lattice
-sys_mgcro = System(xtal_mgcro, dims, spininfos, :dipole); # Same on MgCr2O4 crystal
-
-# To understand what interactions we may assign to this system, we have to
-# understand the the symmetry properties of the crystal, which we turn to next.
-#
-# ### Symmetry analysis for exchange and single-ion anisotropies
-#
-# `print_symmetry_table` reports all the exchange interactions, single-site
-# anisotropies, and g-factors allowed on our crystal. It takes a `Cyrstal` and a
-# distance. We'll look at both the "theorist's" pyrochlore lattice,
-
-print_symmetry_table(xtal_pyro, 5.9)
-
-# and for the the MgCrO4 crystal,
-
-print_symmetry_table(xtal_mgcro, 6.0)
-
-# Note that the exchange anisotropies allowed on the the pyrochlore lattice are
-# slightly different from those on the MgCr2O4 cyrstal due to the role of Mg
-# and O in the bonds. Also note that Sunny has correctly identified the two
-# inequivalent bonds 3a and 3b having the same length. A question may arises to
-# know which bond is J3a and which is J3b, let's plot the structure.
-
-view_crystal(xtal_mgcro, 5.9)
-
-# The crystal viewer shows that the second interaction -- cyan color with
-# distance of 5.89Å -- is in fact the one hopping through a chromium site,
-# meaning it is J3a! We will need to be careful with that later.
-#
-# ### Building the exchange interactions for our system
-#
-# We begin by setting the scale of our exchange interactions on each bond.
-
-J1 = 3.27 # value of J1 in meV from Bai's PRL paper
-J_pyro = [1.00,0.000,0.000,0.000]*J1 # pure nearest neighbor pyrochlore
-J_mgcro = [1.00,0.0815,0.1050,0.085]*J1; # further neighbor pyrochlore relevant for MgCr2O4
-## val_J_mgcro = [1.00,0.000,0.025,0.025]*val_J1; # this is a funny setting from Conlon-Chalker
-
-# These values are then assigned to their corresponding bonds with `set_exchange!`.
-
-## === Assign exchange interactions to pyrochlore system ===
-set_exchange!(sys_pyro, J_pyro[1], Bond(1, 3, [0,0,0])) # J1
-set_exchange!(sys_pyro, J_pyro[2], Bond(1, 2, [0,0,0])) # J2
-set_exchange!(sys_pyro, J_pyro[3], Bond(2, 6, [0,0,0])) # J3a
-set_exchange!(sys_pyro, J_pyro[4], Bond(1, 5, [0,0,0])) # J3b
-
-## === Assign exchange interactions to MgCr2O4 system ===
-set_exchange!(sys_mgcro, J_mgcro[1], Bond(1, 2, [0,0,0])) # J1
-set_exchange!(sys_mgcro, J_mgcro[2], Bond(1, 7, [0,0,0])) # J2
-set_exchange!(sys_mgcro, J_mgcro[3], Bond(1, 3, [0,0,0])) # J3a -- Careful here!
-set_exchange!(sys_mgcro, J_mgcro[4], Bond(1, 3, [1,0,0])); # J3b -- And here!
-
-# We will not be assigning any single-ion anisotropies, so we have finished
-# specifying our models. For good measure, we will set both systems to a random
-# (infinite temperature) initial condition.
-
-randomize_spins!(sys_pyro)
-randomize_spins!(sys_mgcro);
-
-# ## Cooling our `System` amd calculating the instantaneous and dynamic structure factors at the final temperature ##
-#
-# We begin by thermalizing our system at a particular temperature. We will
-# accomplish this by running Langevin dynamics. To do this, we must set up a
-# Langevin integrator.
-
-Δt = 0.05 # Integration time step in meV^-1
-λ = 0.1 # Phenomenological damping parameter
-kT = 1.8 # Desired temperature in meV
-langevin = Langevin(Δt; λ, kT); # Construct integrator
-
-# We can now thermalize our systems by running the integrator.
-for _ in 1:2000
- step!(sys_pyro, langevin)
- step!(sys_mgcro, langevin)
-end
-
-# As a sanity check, we'll plot the real-space spin configurations of both
-# systems after themalization. First the pyrochlore,
-
-plot_spins(sys_pyro)
-
-# and then the MgCr2O4,
-
-plot_spins(sys_mgcro)
-
-# ## Instantaneous Structure Factor
-# Next we can examine the instantaneous structure factor.
-
-isf_pyro = instant_correlations(sys_pyro)
-isf_mgcro = instant_correlations(sys_mgcro);
-
-# These are currently empty. Let's add correlation data from 10 trajectories.
-
-for _ in 1:10
- ## Run dynamics to decorrelate
- for _ in 1:500
- step!(sys_pyro, langevin)
- step!(sys_mgcro, langevin)
- end
- ## Add samples
- add_sample!(isf_pyro, sys_pyro)
- add_sample!(isf_mgcro, sys_mgcro)
-end
-
-# To retrieve the intensities, we set up a formula and call
-# `intensities_interpolated` on an array of wave vectors.
-
-qvals = -4.0:0.025:4.0
-qs = [(qa, qb, 0) for qa in qvals, qb in qvals] # Wave vectors to query
-
-formula_pyro = intensity_formula(isf_pyro, :perp)
-formula_mgcro = intensity_formula(isf_mgcro, :perp)
-Sq_pyro = instant_intensities_interpolated(isf_pyro, qs, formula_pyro)
-Sq_mgcro = instant_intensities_interpolated(isf_mgcro, qs, formula_mgcro);
-
-# Finally, we can plot the results.
-
-fig = Figure(; resolution=(1200,500))
-axparams = (aspect = true, xticks=-4:4, yticks=-4:4, titlesize=20,
- xlabel = "𝐪a (2π a⁻¹)", ylabel = "𝐪b (2π b⁻¹)", xlabelsize = 18, ylabelsize=18,)
-ax_pyro = Axis(fig[1,1]; title="Pyrochlore", axparams...)
-ax_mgcro = Axis(fig[1,3]; title="MgCr2O4", axparams...)
-hm = heatmap!(ax_pyro, qvals, qvals, Sq_pyro)
-Colorbar(fig[1,2], hm)
-hm = heatmap!(ax_mgcro, qvals, qvals, Sq_mgcro)
-Colorbar(fig[1,4], hm)
-fig
-
-# ## Dynamical Structure Factor
-# We can also estimate the dynamical structure factor.
-
-sc_pyro = dynamical_correlations(sys_pyro; Δt, ωmax = 10.0, nω = 100)
-sc_mgcro = dynamical_correlations(sys_mgcro; Δt, ωmax = 10.0, nω = 100);
-
-# Next we add some sample trajectories.
-
-for _ in 1:3
- ## Run dynamics to decorrelate
- for _ in 1:500
- step!(sys_pyro, langevin)
- step!(sys_mgcro, langevin)
- end
- ## Add samples
- add_sample!(sc_pyro, sys_pyro)
- add_sample!(sc_mgcro, sys_mgcro)
-end
-
-# We can now examine the structure factor intensities along a path in momentum space.
-
-fig = Figure(; resolution=(1200,900))
-axsqw = (xticks=-4:4, yticks=0:2:10, ylabel="E (meV)", ylabelsize=18, xlabelsize=18, )
-
-qbs = 0.0:0.5:1.5 # Determine q_b for each slice
-for (i, qb) in enumerate(qbs)
- path, _ = reciprocal_space_path(xtal_pyro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
- ## vectors through the BZ
- formula = intensity_formula(sc_pyro, :perp; kT) # Temperature keyword enables intensity rescaling
- Sqω_pyro = intensities_interpolated(sc_pyro, path, formula)
-
- ax = Axis(fig[fldmod1(i,2)...]; xlabel = "q = (x, $qb, 0)", axsqw...)
- ωs = available_energies(sc_pyro)
- heatmap!(ax, [p[1] for p in path], ωs, Sqω_pyro; colorrange=(0.0, 4.0))
-end
-fig
-
-# And let's take a look at the same slices for MgCr2O4.
-
-fig = Figure(; resolution=(1200,900))
-
-qbs = 0.0:0.5:1.5
-for (i, qb) in enumerate(qbs)
- path, _ = reciprocal_space_path(xtal_mgcro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40) # Generate a path of wave
- ## vectors through the BZ
- formula = intensity_formula(sc_mgcro, :perp; kT) # Temperature keyword enables intensity rescaling
- Sqω_mgcro = intensities_interpolated(sc_mgcro, path, formula)
-
- ax = Axis(fig[fldmod1(i,2)...]; xlabel = "q = (x, $qb, 0)", axsqw...)
- ωs = available_energies(sc_mgcro)
- heatmap!(ax, [p[1] for p in path], ωs, Sqω_mgcro; colorrange=(0.0, 4.0))
-end
-fig
-
-# ### Instantaneous structure factor from a dynamical structure factor
-
-# Finally, we note that the instant structure factor can be calculated from the
-# dynamical structure factor. We simply call `instant_intensities` rather than
-# `intensities`. This will calculate the instantaneous structure factor from
-# from ``𝒮(𝐪,ω)`` by integrating out ``ω`` . An advantage of doing this (as
-# opposed to using `instant_correlations`) is that Sunny is able to apply a
-# temperature- and energy-dependent intensity rescaling before integrating out
-# the dynamical information. The results of this approach are more suitable for
-# comparison with experimental data.
-
-qvals = -4.0:0.05:4.0
-qs = [(qa, qb, 0) for qa in qvals, qb in qvals] # Wave vectors to query
-
-formula_pyro = intensity_formula(sc_pyro, :perp; kT) # Temperature keyword enables intensity rescaling
-formula_mgcro = intensity_formula(sc_mgcro, :perp; kT) # Temperature keyword enables intensity rescaling
-
-Sq_pyro = instant_intensities_interpolated(sc_pyro, qs, formula_pyro)
-Sq_mgcro = instant_intensities_interpolated(sc_mgcro, qs, formula_mgcro);
-
-# We can plot the results below. It is useful to compare these to the plot above
-# generated with an `instant_correlations`.
-
-fig = Figure(; resolution=(1200,500))
-ax_pyro = Axis(fig[1,1]; title="Pyrochlore", axparams...)
-ax_mgcro = Axis(fig[1,3]; title="MgCr2O4", axparams...)
-hm = heatmap!(ax_pyro, qvals, qvals, Sq_pyro)
-Colorbar(fig[1,2], hm)
-hm = heatmap!(ax_mgcro, qvals, qvals, Sq_mgcro)
-Colorbar(fig[1,4], hm)
-fig
\ No newline at end of file
diff --git a/examples/longer_examples/MgCr2O4_160953_2009.cif b/examples/longer_examples/MgCr2O4_160953_2009.cif
deleted file mode 100644
index a2210c746..000000000
--- a/examples/longer_examples/MgCr2O4_160953_2009.cif
+++ /dev/null
@@ -1,254 +0,0 @@
-
-#(C) 2022 by FIZ Karlsruhe - Leibniz Institute for Information Infrastructure. All rights reserved.
-data_160953-ICSD
-_database_code_ICSD 160953
-_audit_creation_date 2009-02-01
-_chemical_name_common 'Magnesium dichromate(III)'
-_chemical_formula_structural 'Mg (Cr2 O4)'
-_chemical_formula_sum 'Cr2 Mg1 O4'
-_chemical_name_structure_type Spinel#MgAl2O4
-_exptl_crystal_density_diffrn 4.42
-_diffrn_ambient_temperature 19.
-_citation_title 'Low temperature neutron diffraction study of Mg Cr2 O4 spinel'
-loop_
-_citation_id
-_citation_journal_full
-_citation_year
-_citation_journal_volume
-_citation_page_first
-_citation_page_last
-_citation_journal_id_ASTM
-primary 'Journal of Physics: Condensed Matter' 2008 20 1 5 JCOMEL
-loop_
-_citation_author_citation_id
-_citation_author_name
-primary 'Ortega-San-Martin, L.'
-primary 'Williams, A.J.'
-primary 'Gordon, C.D.'
-primary 'Klemme, S.'
-primary 'Attfield, J.P.'
-_cell_length_a 8.3329(1)
-_cell_length_b 8.3329(1)
-_cell_length_c 8.3329(1)
-_cell_angle_alpha 90.
-_cell_angle_beta 90.
-_cell_angle_gamma 90.
-_cell_volume 578.61
-_cell_formula_units_Z 8
-_space_group_name_H-M_alt 'F d -3 m Z'
-_space_group_IT_number 227
-loop_
-_space_group_symop_id
-_space_group_symop_operation_xyz
-1 '-z, y+3/4, x+3/4'
-2 'z+3/4, -y, x+3/4'
-3 'z+3/4, y+3/4, -x'
-4 '-z, -y, -x'
-5 'y+3/4, x+3/4, -z'
-6 '-y, x+3/4, z+3/4'
-7 'y+3/4, -x, z+3/4'
-8 '-y, -x, -z'
-9 'x+3/4, -z, y+3/4'
-10 'x+3/4, z+3/4, -y'
-11 '-x, z+3/4, y+3/4'
-12 '-x, -z, -y'
-13 '-z, x+3/4, y+3/4'
-14 'z+3/4, x+3/4, -y'
-15 'z+3/4, -x, y+3/4'
-16 '-z, -x, -y'
-17 'y+3/4, -z, x+3/4'
-18 '-y, z+3/4, x+3/4'
-19 'y+3/4, z+3/4, -x'
-20 '-y, -z, -x'
-21 'x+3/4, y+3/4, -z'
-22 'x+3/4, -y, z+3/4'
-23 '-x, y+3/4, z+3/4'
-24 '-x, -y, -z'
-25 'z, -y+1/4, -x+1/4'
-26 '-z+1/4, y, -x+1/4'
-27 '-z+1/4, -y+1/4, x'
-28 'z, y, x'
-29 '-y+1/4, -x+1/4, z'
-30 'y, -x+1/4, -z+1/4'
-31 '-y+1/4, x, -z+1/4'
-32 'y, x, z'
-33 '-x+1/4, z, -y+1/4'
-34 '-x+1/4, -z+1/4, y'
-35 'x, -z+1/4, -y+1/4'
-36 'x, z, y'
-37 'z, -x+1/4, -y+1/4'
-38 '-z+1/4, -x+1/4, y'
-39 '-z+1/4, x, -y+1/4'
-40 'z, x, y'
-41 '-y+1/4, z, -x+1/4'
-42 'y, -z+1/4, -x+1/4'
-43 '-y+1/4, -z+1/4, x'
-44 'y, z, x'
-45 '-x+1/4, -y+1/4, z'
-46 '-x+1/4, y, -z+1/4'
-47 'x, -y+1/4, -z+1/4'
-48 'x, y, z'
-49 '-z, y+1/4, x+1/4'
-50 'z+3/4, -y+1/2, x+1/4'
-51 'z+3/4, y+1/4, -x+1/2'
-52 '-z, -y+1/2, -x+1/2'
-53 'y+3/4, x+1/4, -z+1/2'
-54 '-y, x+1/4, z+1/4'
-55 'y+3/4, -x+1/2, z+1/4'
-56 '-y, -x+1/2, -z+1/2'
-57 'x+3/4, -z+1/2, y+1/4'
-58 'x+3/4, z+1/4, -y+1/2'
-59 '-x, z+1/4, y+1/4'
-60 '-x, -z+1/2, -y+1/2'
-61 '-z, x+1/4, y+1/4'
-62 'z+3/4, x+1/4, -y+1/2'
-63 'z+3/4, -x+1/2, y+1/4'
-64 '-z, -x+1/2, -y+1/2'
-65 'y+3/4, -z+1/2, x+1/4'
-66 '-y, z+1/4, x+1/4'
-67 'y+3/4, z+1/4, -x+1/2'
-68 '-y, -z+1/2, -x+1/2'
-69 'x+3/4, y+1/4, -z+1/2'
-70 'x+3/4, -y+1/2, z+1/4'
-71 '-x, y+1/4, z+1/4'
-72 '-x, -y+1/2, -z+1/2'
-73 'z, -y+3/4, -x+3/4'
-74 '-z+1/4, y+1/2, -x+3/4'
-75 '-z+1/4, -y+3/4, x+1/2'
-76 'z, y+1/2, x+1/2'
-77 '-y+1/4, -x+3/4, z+1/2'
-78 'y, -x+3/4, -z+3/4'
-79 '-y+1/4, x+1/2, -z+3/4'
-80 'y, x+1/2, z+1/2'
-81 '-x+1/4, z+1/2, -y+3/4'
-82 '-x+1/4, -z+3/4, y+1/2'
-83 'x, -z+3/4, -y+3/4'
-84 'x, z+1/2, y+1/2'
-85 'z, -x+3/4, -y+3/4'
-86 '-z+1/4, -x+3/4, y+1/2'
-87 '-z+1/4, x+1/2, -y+3/4'
-88 'z, x+1/2, y+1/2'
-89 '-y+1/4, z+1/2, -x+3/4'
-90 'y, -z+3/4, -x+3/4'
-91 '-y+1/4, -z+3/4, x+1/2'
-92 'y, z+1/2, x+1/2'
-93 '-x+1/4, -y+3/4, z+1/2'
-94 '-x+1/4, y+1/2, -z+3/4'
-95 'x, -y+3/4, -z+3/4'
-96 'x, y+1/2, z+1/2'
-97 '-z+1/2, y+3/4, x+1/4'
-98 'z+1/4, -y, x+1/4'
-99 'z+1/4, y+3/4, -x+1/2'
-100 '-z+1/2, -y, -x+1/2'
-101 'y+1/4, x+3/4, -z+1/2'
-102 '-y+1/2, x+3/4, z+1/4'
-103 'y+1/4, -x, z+1/4'
-104 '-y+1/2, -x, -z+1/2'
-105 'x+1/4, -z, y+1/4'
-106 'x+1/4, z+3/4, -y+1/2'
-107 '-x+1/2, z+3/4, y+1/4'
-108 '-x+1/2, -z, -y+1/2'
-109 '-z+1/2, x+3/4, y+1/4'
-110 'z+1/4, x+3/4, -y+1/2'
-111 'z+1/4, -x, y+1/4'
-112 '-z+1/2, -x, -y+1/2'
-113 'y+1/4, -z, x+1/4'
-114 '-y+1/2, z+3/4, x+1/4'
-115 'y+1/4, z+3/4, -x+1/2'
-116 '-y+1/2, -z, -x+1/2'
-117 'x+1/4, y+3/4, -z+1/2'
-118 'x+1/4, -y, z+1/4'
-119 '-x+1/2, y+3/4, z+1/4'
-120 '-x+1/2, -y, -z+1/2'
-121 'z+1/2, -y+1/4, -x+3/4'
-122 '-z+3/4, y, -x+3/4'
-123 '-z+3/4, -y+1/4, x+1/2'
-124 'z+1/2, y, x+1/2'
-125 '-y+3/4, -x+1/4, z+1/2'
-126 'y+1/2, -x+1/4, -z+3/4'
-127 '-y+3/4, x, -z+3/4'
-128 'y+1/2, x, z+1/2'
-129 '-x+3/4, z, -y+3/4'
-130 '-x+3/4, -z+1/4, y+1/2'
-131 'x+1/2, -z+1/4, -y+3/4'
-132 'x+1/2, z, y+1/2'
-133 'z+1/2, -x+1/4, -y+3/4'
-134 '-z+3/4, -x+1/4, y+1/2'
-135 '-z+3/4, x, -y+3/4'
-136 'z+1/2, x, y+1/2'
-137 '-y+3/4, z, -x+3/4'
-138 'y+1/2, -z+1/4, -x+3/4'
-139 '-y+3/4, -z+1/4, x+1/2'
-140 'y+1/2, z, x+1/2'
-141 '-x+3/4, -y+1/4, z+1/2'
-142 '-x+3/4, y, -z+3/4'
-143 'x+1/2, -y+1/4, -z+3/4'
-144 'x+1/2, y, z+1/2'
-145 '-z+1/2, y+1/4, x+3/4'
-146 'z+1/4, -y+1/2, x+3/4'
-147 'z+1/4, y+1/4, -x'
-148 '-z+1/2, -y+1/2, -x'
-149 'y+1/4, x+1/4, -z'
-150 '-y+1/2, x+1/4, z+3/4'
-151 'y+1/4, -x+1/2, z+3/4'
-152 '-y+1/2, -x+1/2, -z'
-153 'x+1/4, -z+1/2, y+3/4'
-154 'x+1/4, z+1/4, -y'
-155 '-x+1/2, z+1/4, y+3/4'
-156 '-x+1/2, -z+1/2, -y'
-157 '-z+1/2, x+1/4, y+3/4'
-158 'z+1/4, x+1/4, -y'
-159 'z+1/4, -x+1/2, y+3/4'
-160 '-z+1/2, -x+1/2, -y'
-161 'y+1/4, -z+1/2, x+3/4'
-162 '-y+1/2, z+1/4, x+3/4'
-163 'y+1/4, z+1/4, -x'
-164 '-y+1/2, -z+1/2, -x'
-165 'x+1/4, y+1/4, -z'
-166 'x+1/4, -y+1/2, z+3/4'
-167 '-x+1/2, y+1/4, z+3/4'
-168 '-x+1/2, -y+1/2, -z'
-169 'z+1/2, -y+3/4, -x+1/4'
-170 '-z+3/4, y+1/2, -x+1/4'
-171 '-z+3/4, -y+3/4, x'
-172 'z+1/2, y+1/2, x'
-173 '-y+3/4, -x+3/4, z'
-174 'y+1/2, -x+3/4, -z+1/4'
-175 '-y+3/4, x+1/2, -z+1/4'
-176 'y+1/2, x+1/2, z'
-177 '-x+3/4, z+1/2, -y+1/4'
-178 '-x+3/4, -z+3/4, y'
-179 'x+1/2, -z+3/4, -y+1/4'
-180 'x+1/2, z+1/2, y'
-181 'z+1/2, -x+3/4, -y+1/4'
-182 '-z+3/4, -x+3/4, y'
-183 '-z+3/4, x+1/2, -y+1/4'
-184 'z+1/2, x+1/2, y'
-185 '-y+3/4, z+1/2, -x+1/4'
-186 'y+1/2, -z+3/4, -x+1/4'
-187 '-y+3/4, -z+3/4, x'
-188 'y+1/2, z+1/2, x'
-189 '-x+3/4, -y+3/4, z'
-190 '-x+3/4, y+1/2, -z+1/4'
-191 'x+1/2, -y+3/4, -z+1/4'
-192 'x+1/2, y+1/2, z'
-loop_
-_atom_type_symbol
-_atom_type_oxidation_number
-Mg2+ 2
-Cr3+ 3
-O2- -2
-loop_
-_atom_site_label
-_atom_site_type_symbol
-_atom_site_symmetry_multiplicity
-_atom_site_Wyckoff_symbol
-_atom_site_fract_x
-_atom_site_fract_y
-_atom_site_fract_z
-_atom_site_U_iso_or_equiv
-_atom_site_occupancy
-Mg1 Mg2+ 8 a 0.125 0.125 0.125 0.002 1.
-Cr1 Cr3+ 16 d 0.5 0.5 0.5 0.002 1.
-O1 O2- 32 e 0.2612(1) 0.2612(1) 0.2612(1) 0.002 1.
-#End of TTdata_160953-ICSD
\ No newline at end of file
diff --git a/examples/spinw_ports/08_Kagome_AFM.jl b/examples/spinw_tutorials/SW08_Kagome_AFM.jl
similarity index 98%
rename from examples/spinw_ports/08_Kagome_AFM.jl
rename to examples/spinw_tutorials/SW08_Kagome_AFM.jl
index 4de716581..65ce90507 100644
--- a/examples/spinw_ports/08_Kagome_AFM.jl
+++ b/examples/spinw_tutorials/SW08_Kagome_AFM.jl
@@ -1,4 +1,4 @@
-# # Tutorial 8 - Kagome Antiferromagnet
+# # SW8 - Kagome Antiferromagnet
#
# This is a Sunny port of [SpinW Tutorial
# 8](https://spinw.org/tutorials/08tutorial), originally authored by Bjorn Fak
diff --git a/examples/spinw_ports/15_Ba3NbFe3Si2O14.jl b/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
similarity index 98%
rename from examples/spinw_ports/15_Ba3NbFe3Si2O14.jl
rename to examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
index ac0d8c043..4d547252c 100644
--- a/examples/spinw_ports/15_Ba3NbFe3Si2O14.jl
+++ b/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
@@ -1,4 +1,4 @@
-# # Tutorial 15 - Ba₃NbFe₃Si₂O₁₄
+# # SW15 - Ba₃NbFe₃Si₂O₁₄
#
# This is a Sunny port of [SpinW Tutorial
# 15](https://spinw.org/tutorials/15tutorial), originally authored by Sandor