From 1aa7e6ae7c83f2d72dfa460697af983ada4fe537 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 09:14:24 -0600
Subject: [PATCH 01/13] Docs fine tuning

---
 examples/01_LSWT_CoRh2O4.jl                   | 10 +--
 examples/03_LSWT_SU3_FeI2.jl                  | 45 +++++-----
 examples/04_GSD_FeI2.jl                       | 56 ++++++------
 examples/06_CP2_Skyrmions.jl                  | 85 +++++++++----------
 examples/07_Dipole_Dipole.jl                  | 27 +++---
 examples/08_Momentum_Conventions.jl           | 19 +++--
 examples/extra/Plotting/plotting2d.jl         |  8 +-
 .../spinw_tutorials/SW19_Different_Ions.jl    |  9 ++
 8 files changed, 133 insertions(+), 126 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index cb388b6ba..2f09b04d9 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -75,8 +75,8 @@ sys = System(cryst, [1 => Moment(s=3/2, g=2)], :dipole)
 # Previous work demonstrated that inelastic neutron scattering data for CoRh₂O₄
 # is well described with a single antiferromagnetic nearest neighbor exchange,
 # `J = 0.63` meV. Use [`set_exchange!`](@ref) with the bond that connects atom 1
-# to atom 3, and has zero displacement between chemical cells. Applying the
-# symmetries of spacegroup 227, Sunny will propagate this interaction to the
+# to atom 3, and has zero displacement between chemical cells. Consistent with
+# the symmetries of spacegroup 227, this interaction will be propagated to all
 # other nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now
 # shows the antiferromagnetic Heisenberg interactions as blue polkadot spheres.
 
@@ -180,8 +180,8 @@ plot_intensities(res; units, saturation=1.0)
 
 # ### What's next?
 #
-# * For more spin wave calculations of this traditional type, one can browse the
-#   [SpinW tutorials ported to Sunny](@ref "SW01 - FM Heisenberg chain").
+# * For more spin wave calculations of this type, browse the [SpinW tutorials
+#   ported to Sunny](@ref "SW01 - FM Heisenberg chain").
 # * Spin wave theory neglects thermal fluctuations of the magnetic order. Our
 #   [next tutorial](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*")
 #   demonstrates how to sample spins in thermal equilibrium, and measure
@@ -189,5 +189,5 @@ plot_intensities(res; units, saturation=1.0)
 # * Sunny also offers features that go beyond the dipole approximation of a
 #   quantum spin via the theory of SU(_N_) coherent states. This can be
 #   especially useful for systems with strong single-ion anisotropy, as
-#   demonstrated in our [tutorial on FeI₂](@ref "3. Multi-flavor spin wave
+#   demonstrated in the [FeI₂ tutorial](@ref "3. Multi-flavor spin wave
 #   simulations of FeI₂").
diff --git a/examples/03_LSWT_SU3_FeI2.jl b/examples/03_LSWT_SU3_FeI2.jl
index 77afc49ea..431de64d7 100644
--- a/examples/03_LSWT_SU3_FeI2.jl
+++ b/examples/03_LSWT_SU3_FeI2.jl
@@ -1,8 +1,7 @@
 # # 3. Multi-flavor spin wave simulations of FeI₂
 # 
-# This tutorial illustrates various powerful features in Sunny, including
-# symmetry analysis, energy minimization, and spin wave theory with multi-flavor
-# bosons.
+# This tutorial illustrates various advanced features such as symmetry analysis,
+# energy minimization, and spin wave theory with multi-flavor bosons.
 #
 # Our context will be FeI₂, an effective spin-1 material with strong single-ion
 # anisotropy. Quadrupolar fluctuations give rise to a single-ion bound state
@@ -42,10 +41,10 @@ positions = [[0, 0, 0], [1/3, 2/3, 1/4], [2/3, 1/3, 3/4]]
 types = ["Fe", "I", "I"]
 cryst = Crystal(latvecs, positions; types)
 
-# Observe that Sunny inferred the space group, 'P -3 m 1' (164) and labeled the
-# atoms according to their point group symmetries. Only the Fe atoms are
-# magnetic, so we focus on them with [`subcrystal`](@ref). Importantly, this
-# operation preserves the spacegroup symmetries.
+# Observe that the space group 'P -3 m 1' (164) has been inferred, as well as
+# point group symmetries. Only the Fe atoms are magnetic, so we focus on them
+# with [`subcrystal`](@ref). Importantly, this operation preserves the
+# spacegroup symmetries.
 
 cryst = subcrystal(cryst, "Fe")
 view_crystal(cryst)
@@ -143,14 +142,15 @@ set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)
 
 # ### Finding the ground state
 
-# This model has been carefully designed so that energy minimization produces
-# the physically correct magnetic ordering. Using [`set_dipole!`](@ref), this
-# magnetic structure can be entered manually. Sunny also provides tools to
-# search for an unknown magnetic order, as we will now demonstrate.
+# These model parameters have already been fitted so that energy minimization
+# yields the physically correct ground state. Knowing this, we could set the
+# magnetic configuration manually by calling [`set_dipole!`](@ref) on each site
+# in the system. Another approach, as we now demonstrate, is to use the built-in
+# optimization tools to search for the ground-state in an automated way.
 #
 # To minimize bias in the search, use [`resize_supercell`](@ref) to create a
-# relatively large system of 4×4×4 chemical cells. Randomize all spins (as SU(3)
-# coherent states) and minimize the energy.
+# relatively large system of 4×4×4 chemical cells. Randomize all spins
+# (represented as SU(3) coherent states) and minimize the energy.
 
 sys = resize_supercell(sys, (4, 4, 4))
 randomize_spins!(sys)
@@ -164,11 +164,12 @@ minimize_energy!(sys)
 plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
 # To better understand the spin configuration, we could inspect the static
-# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. For
-# this, Sunny provides [`SampledCorrelationsStatic`](@ref). Here, however, we
-# will use [`print_wrapped_intensities`](@ref), which gives static intensities
-# averaged over the individual Bravais sublattices (in effect, all ``𝐪``
-# intensities are periodically wrapped to the first Brillouin zone).
+# structure factor ``\mathcal{S}(𝐪)`` in the 3D space of momenta ``𝐪``. The
+# general tool for this analysis is [`SampledCorrelationsStatic`](@ref). For the
+# present purposes, however, it is more convenient to use
+# [`print_wrapped_intensities`](@ref), which reports ``\mathcal{S}(𝐪)`` with
+# periodical wrapping of all commensurate ``𝐪`` wavevectors into the first
+# Brillouin zone.
 
 print_wrapped_intensities(sys)
 
@@ -204,14 +205,14 @@ plot_spins(sys_min; color=[S[3] for S in sys_min.dipoles], ghost_radius=12)
 # SU(3) coherent states, this calculation will require a multi-flavor boson
 # generalization of the usual spin wave theory.
 
+swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min))
+
 # Calculate and plot the spectrum along a momentum-space path that connects a
-# sequence of high-symmetry ``𝐪``-points. These Sunny commands are identical to
-# those described in [`our previous CoRh₂O₄ tutorial`](@ref "1. Spin wave
-# simulations of CoRh₂O₄").
+# sequence of high-symmetry ``𝐪``-points. This interface abstracts over the
+# underlying calculator type.
 
 qs = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]
 path = q_space_path(cryst, qs, 500)
-swt = SpinWaveTheory(sys_min; measure=ssf_perp(sys_min))
 res = intensities_bands(swt, path)
 plot_intensities(res; units)
 
diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl
index 8910d9dc0..0f64c1f52 100644
--- a/examples/04_GSD_FeI2.jl
+++ b/examples/04_GSD_FeI2.jl
@@ -2,7 +2,7 @@
 
 # The [previous FeI₂ tutorial](@ref "3. Multi-flavor spin wave simulations of
 # FeI₂") used multi-flavor spin wave theory to calculate the dynamical spin
-# structure factor. Here we perform an analogous calculation at finite
+# structure factor. This tutorial performs an analogous calculation at finite
 # temperature using the [classical dynamics of SU(_N_) coherent
 # states](https://doi.org/10.1103/PhysRevB.106.054423).
 #
@@ -13,14 +13,14 @@
 # enables the study of [highly non-equilibrium
 # processes](https://doi.org/10.1103/PhysRevB.106.235154).
 #
-# The structure of this tutorial largely follows our [previous study of CoRh₂O₄
+# The structure of this tutorial largely follows the [previous study of CoRh₂O₄
 # at finite *T*](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*").
 # The main difference is that CoRh₂O₄ can be well described with `:dipole` mode,
-# whereas FeI₂ has a strong easy-axis anisotropy that introduces a single-ion
-# bound state and necessitates the use of `:SUN` mode.
+# whereas FeI₂ is best modeled using `:SUN` mode, owing to its strong easy-axis
+# anisotropy.
 #
-# Construct the FeI₂ system as described in the [previous tutorial](@ref "3.
-# Multi-flavor spin wave simulations of FeI₂").
+# Construct the FeI₂ system as [previously described](@ref "3. Multi-flavor spin
+# wave simulations of FeI₂").
 
 using Sunny, GLMakie
 
@@ -64,10 +64,12 @@ set_onsite_coupling!(sys, S -> -D*S[3]^2, 1)
 
 sys = resize_supercell(sys, (16, 16, 4))
 
-# Previously, we used [`minimize_energy!`](@ref) to find a local energy minimum.
-# Here, we will instead use [`Langevin`](@ref) dynamics to relax the system into
-# thermal equilibrium. The temperature 2.3 K ≈ 0.2 meV is within the ordered
-# phase, but large enough so that the dynamics can overcome local energy
+# Direct optimization via [`minimize_energy!`](@ref) is susceptible to trapping
+# in a local minimum. An alternative approach is to simulate the system using
+# [`Langevin`](@ref) spin dynamics. This requires a bit more set-up, but allows
+# sampling from thermal equilibrium at any target temperature. Select the
+# temperature 2.3 K ≈ 0.2 meV. This temperature is small enough to magnetically
+# order, but large enough so that the dynamics can readily overcome local energy
 # barriers and annihilate defects.
 
 langevin = Langevin(; damping=0.2, kT=2.3*units.K)
@@ -89,42 +91,42 @@ for _ in 1:10_000
 end
 plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
-# The two-up, two-down spiral order can be verified by calling
+# Verify the expected two-up, two-down spiral magnetic order by calling
 # [`print_wrapped_intensities`](@ref). A single propagation wavevector ``±𝐤``
-# provides most of the static intensity in ``\mathcal{S}(𝐪)``. A smaller amount
-# of intensity is spread among many other wavevectors due to thermal
-# fluctuations.
+# dominates the static intensity in ``\mathcal{S}(𝐪)``, indicating the expected
+# 2 up, 2 down magnetic spiral order. A smaller amount of intensity is spread
+# among many other wavevectors due to thermal fluctuations.
 
 print_wrapped_intensities(sys)
 
-# Calling [`suggest_timestep`](@ref) shows that thermalization has not
-# substantially altered the suggested `dt`.
+# Thermalization has not substantially altered the suggested `dt`.
 
 suggest_timestep(sys, langevin; tol=1e-2)
 
 # ### Structure factor in the paramagnetic phase
 
-# Now we will re-thermalize the system to a temperature of 3.5 K ≈ 0.30 meV,
-# which is in the paramagnetic phase.
+# The remainder of this tutorial will focus on the paramagnetic phase.
+# Re-thermalize the system to the temperature of 3.5 K ≈ 0.30 meV.
 
 langevin.kT = 3.5 * units.K
 for _ in 1:10_000
     step!(sys, langevin)
 end
 
-# At this higher temperature, the suggested timestep has increased slightly.
+# The suggested timestep has increased slightly. Following this suggestion will
+# make the simulations a bit faster.
 
 suggest_timestep(sys, langevin; tol=1e-2)
 langevin.dt = 0.040;
 
-# Now collect dynamical spin structure factor data using a
-# [`SampledCorrelations`](@ref) object. This will involve sampling spin
-# configurations from thermal equilibrium and integrating a [classical spin
-# dynamics for SU(_N_) coherent states](https://arxiv.org/abs/2204.07563).
-# Normal modes appearing in the classical dynamics can be quantized to yield
-# magnetic excitations. The associated structure factor intensities
-# ``S^{αβ}(q,ω)`` can be compared with inelastic neutron scattering data.
-# Incorporate the [`FormFactor`](@ref) appropriate to Fe²⁺. 
+# Collect dynamical spin structure factor data using
+# [`SampledCorrelations`](@ref). This procedure involves sampling spin
+# configurations from thermal equilibrium and using the [spin dynamics of
+# SU(_N_) coherent states](https://arxiv.org/abs/2204.07563) to estimate
+# dynamical correlations. With proper classical-to-quantum corrections, the
+# associated structure factor intensities ``S^{αβ}(q,ω)`` can be compared with
+# finite-temperature inelastic neutron scattering data. Incorporate the
+# [`FormFactor`](@ref) appropriate to Fe²⁺. 
 
 dt = 2*langevin.dt
 energies = range(0, 7.5, 120)
diff --git a/examples/06_CP2_Skyrmions.jl b/examples/06_CP2_Skyrmions.jl
index 403adfa03..ef998dd17 100644
--- a/examples/06_CP2_Skyrmions.jl
+++ b/examples/06_CP2_Skyrmions.jl
@@ -1,44 +1,45 @@
 # # 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
-# of [Zhang et al., Nature Communications **14**, 3626
-# (2023)](https://www.nature.com/articles/s41467-023-39232-8), which supports a
-# novel type of topological defect, a CP² skyrmion, that involves both the
-# dipolar and quadrupolar parts of a quantum spin.
+# This example demonstrates a non-equilibrium study of SU(3) spin dynamics
+# leading to the formation of a CP² skyrmion liquid. As proposed in [Zhang et
+# al., Nature Communications **14**, 3626
+# (2023)](https://www.nature.com/articles/s41467-023-39232-8), CP² skyrmions are
+# topological defects that involve both the dipolar and quadrupolar parts of
+# quantum spin-1, and can be studied using the formalism SU(3) coherent states.
 #
-# Beginning from an initial high-temperature state, a disordered gas of CP²
-# skyrmions can be formed by rapidly quenching to low temperature. To model the
-# coupled dynamics of dipoles and quadrupoles, Sunny uses a recently developed
-# generalization of the Landau-Lifshitz spin dynamics, [Dahlbom et al., Phys.
-# Rev. B **106**, 235154 (2022)](https://doi.org/10.1103/PhysRevB.106.235154).
+# This study uses the SU(N) generalization of Landau-Lifshitz spin dynamics,
+# with Langevin coupling to a thermal bath, as described in [Dahlbom et al.,
+# Phys. Rev. B **106**, 235154
+# (2022)](https://doi.org/10.1103/PhysRevB.106.235154). Beginning from an
+# initial high-temperature state, the dynamics following a rapid quench in
+# temperature gives rise a disordered liquid of CP² skyrmions.
 
 using Sunny, GLMakie
 
-# The Hamiltonian we will implement, 
-# ```math
-# \mathcal{H} = \sum_{\langle i,j \rangle} J_{ij}( \hat{S}_i^x \hat{S}_j^x + \hat{S}_i^y \hat{S}_j^y + \Delta\hat{S}_i^z \hat{S}_j^z) - h\sum_{i}\hat{S}_i^z + D\sum_{i}(\hat{S}_i^z)^2
-# ```
-# contains competing ferromagnetic nearest-neightbor and antiferromagnetic
-# next-nearest-neighbor exchange terms on a triangular lattice. Both exchanges
-# exhibit anisotropy on the z-term. Additionally, there is an external magnetic
-# field, ``h``, and easy-plane single-ion anisotropy, ``D > 0``. We begin by
-# implementing the [`Crystal`](@ref).
+
+# Begin with a [`Crystal`](@ref) cell for the triangular lattice.
 
 lat_vecs = lattice_vectors(1, 1, 10, 90, 90, 120)
-basis_vecs = [[0,0,0]]
-cryst = Crystal(lat_vecs, basis_vecs)
+positions = [[0, 0, 0]]
+cryst = Crystal(lat_vecs, positions)
 
-# Create a spin [`System`](@ref) containing ``L×L`` cells. Selecting
-# [`Units.theory`](@ref Units) with ``g=-1`` provides a dimensionless Zeeman
-# coupling of the form ``-𝐁⋅𝐬``.
+# Create a spin [`System`](@ref) containing ``L×L`` cells. Following previous
+# worse, select ``g=-1`` so that the Zeeman coupling has the form ``-𝐁⋅𝐬``.
 
 L = 40
 sys = System(cryst, [1 => Moment(s=1, g=-1)], :SUN; dims=(L, L, 1))
 
-# We proceed to implement each term of the Hamiltonian, selecting our parameters
-# so that the system occupies a region of the phase diagram that supports
-# skyrmions. The exchange interactions are set as follows.
+# The Hamiltonian,
+# ```math
+# \mathcal{H} = \sum_{\langle i,j \rangle} J_{ij}( \hat{S}_i^x \hat{S}_j^x + \hat{S}_i^y \hat{S}_j^y + \Delta\hat{S}_i^z \hat{S}_j^z) - h\sum_{i}\hat{S}_i^z + D\sum_{i}(\hat{S}_i^z)^2,
+# ```
+# contains competing ferromagnetic first-neighbor and antiferromagnetic
+# second-neighbor exchange terms on a triangular lattice. Both exchange matrices
+# include anisotropy in the $\hat{z}$ direction. Additionally, there is an
+# external magnetic field ``h`` and an easy-plane single-ion anisotropy ``D``.
+# Select parameters for a point in the [Zhang et
+# al.](https://www.nature.com/articles/s41467-023-39232-8) phase diagram where
+# the CP² skyrmions are stable.
 
 J1 = -1           # Nearest-neighbor ferromagnetic
 J2 = (2.0/(1+√5)) # Tune competing exchange to set skyrmion scale length
@@ -53,14 +54,10 @@ ex2 = J2 * [1 0 0;
 set_exchange!(sys, ex1, Bond(1, 1, [1, 0, 0]))
 set_exchange!(sys, ex2, Bond(1, 1, [1, 2, 0]))
 
-# Next we add the external field,
-
-h = 15.5
+h = 15.5         # External field in energy units
 field = set_field!(sys, [0, 0, h])
 
-# and finally an easy-plane single-ion anisotropy,
-
-D = 19.0
+D = 19.0         # Easy-plane anisotropy
 set_onsite_coupling!(sys, S -> D*S[3]^2, 1)
 
 # Initialize system to an infinite temperature (fully randomized) initial
@@ -88,9 +85,9 @@ suggest_timestep(sys, integrator; tol=0.05)
 
 integrator.dt = 0.01;
 
-# Now run the dynamical quench starting from a randomized configuration. We will
-# record the state of the system at three different times during the quenching
-# process by copying the `coherents` field of the `System`.
+# Now run the dynamical quench starting from a randomized configuration. The
+# field `frames` stores the system spin configuration, as SU(3) coherent states,
+# at three different times during the quenching process.
 
 randomize_spins!(sys)
 τs = [4, 16, 256]   # Times to record snapshots
@@ -104,12 +101,12 @@ for i in eachindex(τs)
     push!(frames, copy(sys.coherents))     # Save a snapshot spin configuration
 end
 
-# To visualize the state of the system contained in each snapshot, we will
-# calculate and plot the skyrmion density on each plaquette of our lattice. The
-# 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/Plotting).
+# Visualize the state of the system contained in each snapshot by plotting the
+# SU(3) Berry phase curvature over triangular plaquettes. This is a measure of
+# CP² skyrmion density. The function `plot_triangular_plaquettes` is not part of
+# the core Sunny package, but rather something you could define yourself using
+# Makie. Find this helper function at
+# [`examples/extra/Plotting/`](https://github.com/SunnySuite/Sunny.jl/tree/main/examples/extra/Plotting).
 
 include(pkgdir(Sunny, "examples", "extra", "Plotting", "plotting2d.jl"))
 
@@ -121,7 +118,7 @@ function sun_berry_curvature(z₁, z₂, z₃)
     return angle(n₁ * n₂ * n₃)
 end
 
-plot_triangular_plaquettes(sun_berry_curvature, frames; size=(600,200),
+plot_triangular_plaquettes(sun_berry_curvature, frames; size=(600, 200),
     offset_spacing=10, texts=["\tt = "*string(τ) for τ in τs], text_offset=(0, 6)
 )
 
diff --git a/examples/07_Dipole_Dipole.jl b/examples/07_Dipole_Dipole.jl
index 34dac2df1..1c1e2a42a 100644
--- a/examples/07_Dipole_Dipole.jl
+++ b/examples/07_Dipole_Dipole.jl
@@ -1,18 +1,21 @@
 # # 7. Long-range dipole interactions
 #
-# This example demonstrates Sunny's ability to incorporate long-range dipole
-# interactions using Ewald summation. The calculation reproduces previous
-# results in  [Del Maestro and Gingras, J. Phys.: Cond. Matter, **16**, 3339
+# This example shows how long-range dipole-dipole interactions can affect a spin
+# wave calculation. These interactions can be included two ways: Ewald summation
+# or real-space space with a distance cutoff. The study follows [Del Maestro and
+# Gingras, J. Phys.: Cond. Matter, **16**, 3339
 # (2004)](https://arxiv.org/abs/cond-mat/0403494).
 
 using Sunny, GLMakie
 
 # Create a Pyrochlore crystal from spacegroup 227.
 
-units = Units(:K)
+units = Units(:K, :angstrom)
 latvecs = lattice_vectors(10.19, 10.19, 10.19, 90, 90, 90)
-positions = [[0, 0, 0]]
-cryst = Crystal(latvecs, positions, 227, setting="2")
+positions = [[1/8, 1/8, 1/8]]
+cryst = Crystal(latvecs, positions, 227, setting="1")
+view_crystal(cryst)
+
 
 sys = System(cryst, [1 => Moment(s=7/2, g=2)], :dipole)
 J1 = 0.304 # (K)
@@ -25,16 +28,16 @@ set_exchange!(sys, J1, Bond(1, 2, [0,0,0]))
 shape = [1/2 1/2 0; 0 1/2 1/2; 1/2 0 1/2]
 sys_prim = reshape_supercell(sys, shape)
 
-set_dipole!(sys_prim, [+1, -1, 0], position_to_site(sys_prim, [0, 0, 0]))
-set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [1/4, 1/4, 0]))
-set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [1/4, 0, 1/4]))
-set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [0, 1/4, 1/4]))
+set_dipole!(sys_prim, [+1, -1, 0], position_to_site(sys_prim, [1/8, 1/8, 1/8]))
+set_dipole!(sys_prim, [-1, +1, 0], position_to_site(sys_prim, [3/8, 3/8, 1/8]))
+set_dipole!(sys_prim, [+1, +1, 0], position_to_site(sys_prim, [3/8, 1/8, 3/8]))
+set_dipole!(sys_prim, [-1, -1, 0], position_to_site(sys_prim, [1/8, 3/8, 3/8]))
 
 plot_spins(sys_prim; ghost_radius=8, color=[:red, :blue, :yellow, :purple])
 
 # Calculate dispersions with and without long-range dipole interactions. The
-# high-symmetry k-points are specified with respect to the conventional cubic
-# cell.
+# high-symmetry ``𝐪``-points are specified with respect to the conventional
+# cubic cell.
 
 qs = [[0,0,0], [0,1,0], [1,1/2,0], [1/2,1/2,1/2], [3/4,3/4,0], [0,0,0]]
 labels = ["Γ", "X", "W", "L", "K", "Γ"]
diff --git a/examples/08_Momentum_Conventions.jl b/examples/08_Momentum_Conventions.jl
index 46c98e9b0..5586fb8fd 100644
--- a/examples/08_Momentum_Conventions.jl
+++ b/examples/08_Momentum_Conventions.jl
@@ -2,13 +2,14 @@
 #
 # This example illustrates Sunny's conventions for dynamical structure factor
 # intensities, ``\mathcal{S}(𝐪,ω)``, as documented in the page [Structure
-# Factor Conventions](@ref). In the neutron scattering context, the variables
-# ``𝐪`` and ``ω`` describe momentum and energy transfer _to_ the sample.
+# Factor Conventions](@ref). The variables ``𝐪`` and ``ω`` describe momentum
+# and energy transfer _to_ the sample.
 #
 # For systems without inversion-symmetry, the structure factor intensities at
-# ``± 𝐪`` may be inequivalent. To highlight this inequivalence, we will
-# construct a 1D chain with Dzyaloshinskii–Moriya interactions between nearest
-# neighbor bonds, and apply a magnetic field.
+# ``± 𝐪`` may be inequivalent. To highlight this, consider a simple 1D chain
+# that includes only Dzyaloshinskii–Moriya interactions between neighboring
+# sites. Coupling to an external field then breaks time-reversal symmetry,
+# giving rise to an inequivalence of intensities ``\mathcal{S}(±𝐪,ω)``
 
 using Sunny, GLMakie
 
@@ -20,10 +21,10 @@ using Sunny, GLMakie
 latvecs = lattice_vectors(2, 2, 1, 90, 90, 90)
 cryst = Crystal(latvecs, [[0,0,0]], "P1")
 
-# Construct a 1D chain system that extends along ``𝐚₃``. The Hamiltonian
-# includes DM and Zeeman coupling terms, ``ℋ = ∑_j D ẑ ⋅ (𝐒_j × 𝐒_{j+1}) -
-# ∑_j 𝐁 ⋅ μ_j``, where ``μ_j = - μ_B g 𝐒_j`` is the [`magnetic_moment`](@ref)
-# and ``𝐁 ∝ ẑ``.
+# Construct a 1D chain system that extends along the global Cartesian ``ẑ``
+# axis. The Hamiltonian includes DM and Zeeman coupling terms, ``ℋ = ∑_j D ẑ ⋅
+# (𝐒_j × 𝐒_{j+1}) - ∑_j 𝐁 ⋅ μ_j``, where ``μ_j = - g 𝐒_j`` is the
+# [`magnetic_moment`](@ref) and ``𝐁 ∝ ẑ``.
 
 sys = System(cryst, [1 => Moment(s=1, g=2)], :dipole; dims=(1, 1, 25))
 D = 0.1
diff --git a/examples/extra/Plotting/plotting2d.jl b/examples/extra/Plotting/plotting2d.jl
index e7a16e5d9..4d0dea391 100644
--- a/examples/extra/Plotting/plotting2d.jl
+++ b/examples/extra/Plotting/plotting2d.jl
@@ -1,10 +1,7 @@
 using LinearAlgebra
 
 
-#################################################
-# Functions for plotting on triangular plaquettes
-#################################################
-
+## Plotting on triangular plaquettes
 
 function plaquette_idcs(dims::Tuple{Int,Int,Int})
     dx, dy, dz = dims
@@ -47,9 +44,6 @@ function plaquette_map(f::Function, x::Array{T,4}) where T
 end
 
 
-################################################################################
-# Plotting functions
-################################################################################
 function aspect_ratio(x_panel, y_panel, x_offset, y_offset, numrows, numcols;
     adhoc_offset = (0.0, 0.0)
 )
diff --git a/examples/spinw_tutorials/SW19_Different_Ions.jl b/examples/spinw_tutorials/SW19_Different_Ions.jl
index 9a0fdab01..621e02213 100644
--- a/examples/spinw_tutorials/SW19_Different_Ions.jl
+++ b/examples/spinw_tutorials/SW19_Different_Ions.jl
@@ -60,3 +60,12 @@ res = intensities_bands(swt, path)
 plot_intensities!(fig[2, 2], res; units, axisopts=(; title="Fe-Fe correlations"))
 
 fig
+
+# Calculate quantum corrections ``δS`` to spin magnitude, which arise from the
+# zero-point energy of the spin waves. The outputs are ordered following the
+# [`Site`](@ref) indexing scheme of the system: `(cell1, cell2, cell3,
+# sublattice)`, with left-side indices fastest. The two corrections ``δS ≈
+# -0.137`` and ``δS ≈ -0.578`` apply to the Cu and Fe ions, respectively. The
+# larger correction on Fe is due to the relatively weak interchain coupling.
+
+Sunny.magnetization_lswt_correction_dipole(swt; atol=1e-4)

From 5e6737828da1ceff6951eef53a4dd6f2037cc18f Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 09:39:50 -0600
Subject: [PATCH 02/13] Base.show implementations

---
 src/Measurements/Broadening.jl  |  8 ++++++++
 src/Measurements/MeasureSpec.jl | 11 +++++++++++
 2 files changed, 19 insertions(+)

diff --git a/src/Measurements/Broadening.jl b/src/Measurements/Broadening.jl
index 20ea0a9b5..14393e056 100644
--- a/src/Measurements/Broadening.jl
+++ b/src/Measurements/Broadening.jl
@@ -18,6 +18,14 @@ function (b::NonstationaryBroadening)(ϵ, ω)
     b.kernel(ϵ, ω)
 end
 
+function Base.show(io::IO, ::Broadening)
+    print(io, "Broadening")
+end
+
+function Base.show(io::IO, ::NonstationaryBroadening)
+    print(io, "NonstationaryBroadening")
+end
+
 """
     lorentzian(; fwhm)
 
diff --git a/src/Measurements/MeasureSpec.jl b/src/Measurements/MeasureSpec.jl
index 0249e4dae..2da0f9aae 100644
--- a/src/Measurements/MeasureSpec.jl
+++ b/src/Measurements/MeasureSpec.jl
@@ -22,6 +22,17 @@ struct MeasureSpec{Op <: Union{Vec3, HermitianC64}, F, Ret}
     end
 end
 
+function Base.show(io::IO, ::MeasureSpec)
+    print(io, "MeasureSpec")
+end
+
+function Base.show(io::IO, ::MIME"text/plain", m::MeasureSpec)
+    nobs = size(m.observables, 1)
+    ret = eltype(m)
+    println(io, "MeasureSpec [$nobs observables, returns $ret]")
+end
+
+
 Base.eltype(::MeasureSpec{Op, F, Ret}) where {Op, F, Ret} = Ret
 
 # Replicate Sam's observables code for the moment

From a3950fa740e6dafc8d67ae08f07bfb6e104403c6 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 09:40:02 -0600
Subject: [PATCH 03/13] Fix reported energy bins

---
 src/SampledCorrelations/SampledCorrelations.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/SampledCorrelations/SampledCorrelations.jl b/src/SampledCorrelations/SampledCorrelations.jl
index 2a2265ca4..9c8804ab9 100644
--- a/src/SampledCorrelations/SampledCorrelations.jl
+++ b/src/SampledCorrelations/SampledCorrelations.jl
@@ -212,11 +212,13 @@ end
 
 function Base.show(io::IO, ::MIME"text/plain", sc::SampledCorrelations)
     (; crystal, sys_dims, nsamples) = sc
+    nω = Int(size(sc.data, 7)/2+1)
+    Δω = round(sc.Δω, digits=4)
     printstyled(io, "SampledCorrelations"; bold=true, color=:underline)
     println(io," ($(Base.format_bytes(Base.summarysize(sc))))")
     print(io,"[")
     printstyled(io,"S(q,ω)"; bold=true)
-    print(io," | nω = $(round(Int, size(sc.data)[7]/2)), Δω = $(round(sc.Δω, digits=4))")
+    print(io," | nω = $nω, Δω = $Δω")
     println(io," | $nsamples $(nsamples > 1 ? "samples" : "sample")]")
     println(io, supercell_to_str(sys_dims, crystal))
 end

From 72b9788c37c89043298432c62a5e7dd13b02d286 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 13:22:54 -0600
Subject: [PATCH 04/13] Docs tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 27 +++++++++++++++------------
 examples/04_GSD_FeI2.jl     |  9 ++++-----
 2 files changed, 19 insertions(+), 17 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index 2f09b04d9..aa1e2ae17 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -105,12 +105,14 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
 # ### Reshaping the magnetic cell
 
-# The same Néel order can also be described with a magnetic cell that consists
-# of the 2 cobalt atoms in the primitive cell. Columns of the 3×3 `shape` matrix
-# below are the primitive lattice vectors in units of the conventional, cubic
-# lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. Use [`reshape_supercell`](@ref) to
-# construct a system with this shape, and verify that the energy per site is
-# unchanged.
+# The Néel magnetic order can be concisely described with a magnetic cell that
+# coincides with a primitive unit cell. Reshape the system to this primitive
+# cell using the command [`reshape_supercell`](@ref). Columns of the `shape`
+# matrix denote primitive lattice vectors as multiples of the conventional cubic
+# lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. Alternatively, a primitive cell shape
+# can be obtained as `shape = cryst.latvecs \
+# cryst.prim_latvecs`. Verify that the energy per site is unchanged after the
+# reshaping the supercell.
 
 shape = [0 1 1;
          1 0 1;
@@ -158,12 +160,13 @@ energies = range(0, 6, 300)
 res = intensities(swt, path; energies, kernel)
 plot_intensities(res; units)
 
-# To directly compare with the available experimental data, perform a
-# [`powder_average`](@ref) over all possible crystal orientations. Consider 200
-# ``𝐪`` magnitudes ranging from 0 to 3 inverse angstroms. Each magnitude
-# defines spherical shell in reciprocal space, to be sampled with `2000`
-# ``𝐪``-points. The calculation completes in just a couple seconds because the
-# magnetic cell size is small.
+# Sometimes experimental data is only available as a powder average, i.e., as an
+# average over all possible crystal orientations. Use [`powder_average`](@ref)
+# to simulate these intensities. Each ``𝐪``-magnitude defines a spherical shell
+# in reciprocal space. Consider 200 radii from 0 to 3 inverse angstroms, and
+# collect `2000` random samples per spherical shell. As configured, this
+# calculation completes in about two seconds. Had we used the conventional cubic
+# cell, the calculation would be about ``4^3`` times slower.
 
 radii = range(0, 3, 200) # (1/Å)
 res = powder_average(cryst, radii, 2000) do qs
diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl
index 0f64c1f52..e03b4e1db 100644
--- a/examples/04_GSD_FeI2.jl
+++ b/examples/04_GSD_FeI2.jl
@@ -149,11 +149,10 @@ for _ in 1:2
     add_sample!(sc, sys)
 end
 
-# Measure intensities along a path connecting high-symmetry ``𝐪``-points,
-# specified in reciprocal lattice units (RLU). A classical-to-quantum rescaling
-# of normal mode occupations will be performed according to the temperature
-# `kT`. The large statistical noise could be reduced by averaging over more
-# thermal samples.
+# Perform an intensity calculation for two special ``𝐪``-points in reciprocal
+# lattice units (RLU). A classical-to-quantum rescaling of normal mode
+# occupations will be performed according to the temperature `kT`. The large
+# statistical noise could be reduced by averaging over more thermal samples.
 
 res = intensities(sc, [[0, 0, 0], [0.5, 0.5, 0.5]]; energies, langevin.kT)
 fig = lines(res.energies, res.data[:, 1]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)")

From 4ec5fbd0d92fbbc459268172fdec7014077909f3 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 14:00:35 -0600
Subject: [PATCH 05/13] Tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index aa1e2ae17..ae7198ead 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -106,11 +106,11 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])
 # ### Reshaping the magnetic cell
 
 # The Néel magnetic order can be concisely described with a magnetic cell that
-# coincides with a primitive unit cell. Reshape the system to this primitive
-# cell using the command [`reshape_supercell`](@ref). Columns of the `shape`
-# matrix denote primitive lattice vectors as multiples of the conventional cubic
-# lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. Alternatively, a primitive cell shape
-# can be obtained as `shape = cryst.latvecs \
+# coincides with a primitive unit cell. Reshape to a primitive cell with
+# [`reshape_supercell`](@ref). Columns of the `shape` matrix denote primitive
+# lattice vectors as multiples of the conventional cubic lattice vectors
+# ``(𝐚_1, 𝐚_2, 𝐚_3)``. Alternatively, a primitive cell shape can be obtained
+# as `shape = cryst.latvecs \
 # cryst.prim_latvecs`. Verify that the energy per site is unchanged after the
 # reshaping the supercell.
 

From 2cf14d7f8ebfaabeb0e067aa81bb95cc5ac19f6e Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 16:17:22 -0600
Subject: [PATCH 06/13] Tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index ae7198ead..df716b3ad 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -105,12 +105,11 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
 # ### Reshaping the magnetic cell
 
-# The Néel magnetic order can be concisely described with a magnetic cell that
-# coincides with a primitive unit cell. Reshape to a primitive cell with
-# [`reshape_supercell`](@ref). Columns of the `shape` matrix denote primitive
-# lattice vectors as multiples of the conventional cubic lattice vectors
-# ``(𝐚_1, 𝐚_2, 𝐚_3)``. Alternatively, a primitive cell shape can be obtained
-# as `shape = cryst.latvecs \
+# This Néel order can be captured using a primitive unit cell for the magnetic
+# cell. Build such a system with [`reshape_supercell`](@ref), where columns of
+# the `shape` matrix denote primitive lattice vectors as multiples of the
+# conventional cubic lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. One could also use
+# `shape = cryst.latvecs \
 # cryst.prim_latvecs`. Verify that the energy per site is unchanged after the
 # reshaping the supercell.
 

From 2785ad799b2e5a8c3b53633787220fb980b0f122 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Wed, 28 Aug 2024 18:24:59 -0600
Subject: [PATCH 07/13] Formatting

---
 examples/01_LSWT_CoRh2O4.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index df716b3ad..099793bcf 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -109,9 +109,8 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])
 # cell. Build such a system with [`reshape_supercell`](@ref), where columns of
 # the `shape` matrix denote primitive lattice vectors as multiples of the
 # conventional cubic lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. One could also use
-# `shape = cryst.latvecs \
-# cryst.prim_latvecs`. Verify that the energy per site is unchanged after the
-# reshaping the supercell.
+# `shape = cryst.latvecs \ cryst.prim_latvecs`. Verify that the energy per site
+# is unchanged after the reshaping the supercell.
 
 shape = [0 1 1;
          1 0 1;

From 048c992f279bf00b432170ea752f1362dadb530d Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 07:48:54 -0600
Subject: [PATCH 08/13] Tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 42 ++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index 099793bcf..4ca55d6c4 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -1,8 +1,8 @@
 # # 1. Spin wave simulations of CoRh₂O₄
 #
 # This tutorial introduces Sunny through its features for performing
-# conventional spin wave theory calculations. For concreteness, we consider the
-# crystal CoRh₂O₄ and reproduce the calculations of [Ge et al., Phys. Rev. B 96,
+# conventional spin wave theory calculations. We consider the crystal CoRh₂O₄
+# and reproduce the calculations of [Ge et al., Phys. Rev. B 96,
 # 064413](https://doi.org/10.1103/PhysRevB.96.064413).
 
 # ### Get Julia and Sunny
@@ -72,13 +72,13 @@ view_crystal(cryst)
 
 sys = System(cryst, [1 => Moment(s=3/2, g=2)], :dipole)
 
-# Previous work demonstrated that inelastic neutron scattering data for CoRh₂O₄
-# is well described with a single antiferromagnetic nearest neighbor exchange,
-# `J = 0.63` meV. Use [`set_exchange!`](@ref) with the bond that connects atom 1
-# to atom 3, and has zero displacement between chemical cells. Consistent with
-# the symmetries of spacegroup 227, this interaction will be propagated to all
-# other nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now
-# shows the antiferromagnetic Heisenberg interactions as blue polkadot spheres.
+# Ge et al. demonstrated that inelastic neutron scattering data for CoRh₂O₄ is
+# well modeled by antiferromagnetic nearest neighbor exchange, `J = 0.63` meV.
+# Call [`set_exchange!`](@ref) with the bond that connects atom 1 to atom 3, and
+# has zero displacement between chemical cells. Consistent with the symmetries
+# of spacegroup 227, this interaction will be propagated to all other
+# nearest-neighbor bonds. Calling [`view_crystal`](@ref) with `sys` now shows
+# the antiferromagnetic Heisenberg interactions as blue polkadot spheres.
 
 J = +0.63 # (meV)
 set_exchange!(sys, J, Bond(1, 3, [0, 0, 0]))
@@ -105,12 +105,12 @@ plot_spins(sys; color=[S[3] for S in sys.dipoles])
 
 # ### Reshaping the magnetic cell
 
-# This Néel order can be captured using a primitive unit cell for the magnetic
-# cell. Build such a system with [`reshape_supercell`](@ref), where columns of
-# the `shape` matrix denote primitive lattice vectors as multiples of the
-# conventional cubic lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. One could also use
-# `shape = cryst.latvecs \ cryst.prim_latvecs`. Verify that the energy per site
-# is unchanged after the reshaping the supercell.
+# The most compact magnetic cell for this Néel order is a primitive unit cell.
+# Reduce the magnetic cell size using [`reshape_supercell`](@ref), where columns
+# of the `shape` matrix are primitive lattice vectors as multiples of the
+# conventional cubic lattice vectors ``(𝐚_1, 𝐚_2, 𝐚_3)``. One could
+# alternatively use `shape = cryst.latvecs \ cryst.prim_latvecs`. Verify that
+# the energy per site is unchanged after the reshaping the supercell.
 
 shape = [0 1 1;
          1 0 1;
@@ -118,8 +118,8 @@ shape = [0 1 1;
 sys_prim = reshape_supercell(sys, shape)
 @assert energy_per_site(sys_prim) ≈ -2J*(3/2)^2
 
-# Plotting the spins of `sys_prim` shows the primitive cell as a gray wireframe
-# inside the conventional cubic cell.
+# Plotting `sys_prim` shows the two spins within the primitive cell, as well as
+# the larger conventional cubic cell for context.
 
 plot_spins(sys_prim; color=[S[3] for S in sys_prim.dipoles])
 
@@ -183,10 +183,10 @@ plot_intensities(res; units, saturation=1.0)
 #
 # * For more spin wave calculations of this type, browse the [SpinW tutorials
 #   ported to Sunny](@ref "SW01 - FM Heisenberg chain").
-# * Spin wave theory neglects thermal fluctuations of the magnetic order. Our
-#   [next tutorial](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at finite *T*")
-#   demonstrates how to sample spins in thermal equilibrium, and measure
-#   dynamical correlations from the classical spin dynamics.
+# * Spin wave theory neglects thermal fluctuations of the magnetic order. The
+#   [next CoRh₂O₄ tutorial](@ref "2. Landau-Lifshitz dynamics of CoRh₂O₄ at
+#   finite *T*") demonstrates how to sample spins in thermal equilibrium, and
+#   measure dynamical correlations from the classical spin dynamics.
 # * Sunny also offers features that go beyond the dipole approximation of a
 #   quantum spin via the theory of SU(_N_) coherent states. This can be
 #   especially useful for systems with strong single-ion anisotropy, as

From 460859fc642fab5e3db1131e0d85a41bf662bb3f Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 07:50:28 -0600
Subject: [PATCH 09/13] Tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index 4ca55d6c4..cb5b88c09 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -26,9 +26,8 @@ using Sunny, GLMakie
 # ### Units
 
 # The [`Units`](@ref Sunny.Units) object selects reference energy and length
-# scales. This is achieved by providing physical constants. For example,
-# `units.K` would provide one kelvin as 0.086 meV, with the Boltzmann constant
-# implicit.
+# scales, and uses these to provide physical constants. For example, `units.K`
+# would provide one kelvin as 0.086 meV, with the Boltzmann constant implicit.
 
 units = Units(:meV, :angstrom);
 

From f76c014c4ef4d5e9e948c9735d14e0733e28c8df Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 08:31:45 -0600
Subject: [PATCH 10/13] More fixes

---
 examples/01_LSWT_CoRh2O4.jl                   |  2 +-
 .../SW01_FM_Heseinberg_chain.jl               | 10 ++++-----
 .../spinw_tutorials/SW05_Simple_kagome_FM.jl  |  4 ++--
 .../spinw_tutorials/SW08_sqrt3_kagome_AFM.jl  |  6 +++---
 .../spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl    | 18 ++++++----------
 .../spinw_tutorials/SW18_Distorted_kagome.jl  | 16 +++++++-------
 .../spinw_tutorials/SW19_Different_Ions.jl    | 21 ++++++++++++-------
 src/SpinWaveTheory/SpinWaveTheory.jl          |  4 ++--
 8 files changed, 40 insertions(+), 41 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index cb5b88c09..d6e7ebdf3 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -163,7 +163,7 @@ plot_intensities(res; units)
 # in reciprocal space. Consider 200 radii from 0 to 3 inverse angstroms, and
 # collect `2000` random samples per spherical shell. As configured, this
 # calculation completes in about two seconds. Had we used the conventional cubic
-# cell, the calculation would be about ``4^3`` times slower.
+# cell, the calculation would be an order of magnitude slower.
 
 radii = range(0, 3, 200) # (1/Å)
 res = powder_average(cryst, radii, 2000) do qs
diff --git a/examples/spinw_tutorials/SW01_FM_Heseinberg_chain.jl b/examples/spinw_tutorials/SW01_FM_Heseinberg_chain.jl
index 552b720e3..0f430af68 100644
--- a/examples/spinw_tutorials/SW01_FM_Heseinberg_chain.jl
+++ b/examples/spinw_tutorials/SW01_FM_Heseinberg_chain.jl
@@ -36,11 +36,11 @@ view_crystal(cryst; ndims=2, ghost_radius=8)
 
 print_symmetry_table(cryst, 8)
 
-# Use the chemical cell to create a spin [`System`](@ref) with spin S=1 and a
-# magnetic form factor for Cu¹⁺. In this case, it is only necessary to simulate
-# a single chemical cell. The option `:dipole` indicates that, following
-# traditional spin wave theory, we are modeling quantum spins using only their
-# expected dipole moments.
+# Use the chemical cell to create a spin [`System`](@ref) with spin ``s = 1``
+# and a magnetic form factor for Cu¹⁺. In this case, it is only necessary to
+# simulate a single chemical cell. The option `:dipole` indicates that,
+# following traditional spin wave theory, we are modeling quantum spins using
+# only their expected dipole moments.
 
 sys = System(cryst, [1 => Moment(s=1, g=2)], :dipole)
 
diff --git a/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl b/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl
index 232fec5a3..4355beae2 100644
--- a/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl
+++ b/examples/spinw_tutorials/SW05_Simple_kagome_FM.jl
@@ -56,8 +56,8 @@ res = intensities_bands(swt, path)
 plot_intensities(res; units)
 
 # Calculate and plot the powder average with two different magnitudes of
-# Gaussian line-broadening. Pick an explicit colorrange so that the two color
-# scales are consistent.
+# Gaussian line-broadening. Pick an explicit intensity `colorrange` (as a
+# density in meV) so that the two color scales are consistent.
 
 radii = range(0, 2.5, 200)
 energies = range(0, 6.5, 200)
diff --git a/examples/spinw_tutorials/SW08_sqrt3_kagome_AFM.jl b/examples/spinw_tutorials/SW08_sqrt3_kagome_AFM.jl
index 97b743d10..41a7d6b1d 100644
--- a/examples/spinw_tutorials/SW08_sqrt3_kagome_AFM.jl
+++ b/examples/spinw_tutorials/SW08_sqrt3_kagome_AFM.jl
@@ -25,9 +25,9 @@ set_exchange!(sys, J, Bond(2, 3, [0, 0, 0]))
 # Initialize to an energy minimizing magnetic structure, for which
 # nearest-neighbor spins are at 120° angles.
 
-set_dipole!(sys, [cos(0),sin(0),0], (1, 1, 1, 1))
-set_dipole!(sys, [cos(0),sin(0),0], (1, 1, 1, 2))
-set_dipole!(sys, [cos(2π/3),sin(2π/3),0], (1, 1, 1, 3))
+set_dipole!(sys, [cos(0), sin(0), 0], (1, 1, 1, 1))
+set_dipole!(sys, [cos(0), sin(0), 0], (1, 1, 1, 2))
+set_dipole!(sys, [cos(2π/3), sin(2π/3), 0], (1, 1, 1, 3))
 k = [-1/3, -1/3, 0]
 axis = [0, 0, 1]
 sys_enlarged = repeat_periodically_as_spiral(sys, (3, 3, 1); k, axis)
diff --git a/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl b/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
index d85b8a5d9..fc9b3e767 100644
--- a/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
+++ b/examples/spinw_tutorials/SW15_Ba3NbFe3Si2O14.jl
@@ -37,16 +37,10 @@ set_exchange!(sys, J₁, Bond(3, 2, [1,1,0]))
 set_exchange!(sys, J₄, Bond(1, 1, [0,0,1]))
 set_exchange!(sys, J₂, Bond(1, 3, [0,0,0]))
 
-# The final two exchanges define the chirality of the magnetic structure. The
-# crystal chirality, ``\epsilon_T``, the chirality of each triangle, ``ϵ_D`` and
-# the sense of rotation of the spin helices along ``c``, ``ϵ_H``. The three
-# chiralities are related by ``ϵ_T=ϵ_D ϵ_H``. We now assign ``J_3`` and ``J_5``
-# according to the crystal chirality.
-
-ϵD = -1
-ϵH = +1
-ϵT = ϵD * ϵH
+# The final two exchanges are setting according to the desired chirality ``ϵ_T``
+# of the magnetic structure.
 
+ϵT = -1
 if ϵT == -1
     set_exchange!(sys, J₃, Bond(2, 3, [-1,-1,1]))
     set_exchange!(sys, J₅, Bond(3, 2, [1,1,1]))
@@ -54,7 +48,7 @@ elseif ϵT == 1
     set_exchange!(sys, J₅, Bond(2, 3, [-1,-1,1]))
     set_exchange!(sys, J₃, Bond(3, 2, [1,1,1]))
 else
-    throw("Provide a valid chirality")
+    error("Chirality must be ±1")
 end
 
 # This compound is known to have a spiral order with approximate propagation
@@ -88,7 +82,7 @@ qs = [[0, 1, -1], [0, 1, -1+1], [0, 1, -1+2], [0, 1, -1+3]]
 path = q_space_path(cryst, qs, 400)
 energies = range(0, 6, 400)
 res = intensities(swt, path; energies, kernel=gaussian(fwhm=0.25))
-axisopts = (; title=L"$ϵ_T=-1$, $ϵ_Δ=-1$, $ϵ_H=+1$", titlesize=20)
+axisopts = (; title=L"$ϵ_T = %$ϵT$", titlesize=20)
 plot_intensities(res; units, axisopts, saturation=0.7, colormap=:jet)
 
 # Use [`ssf_custom_bm`](@ref) to calculate the imaginary part of
@@ -102,5 +96,5 @@ measure = ssf_custom_bm(sys; u=[0, 1, 0], v=[0, 0, 1]) do q, ssf
 end
 swt = SpinWaveTheorySpiral(sys; measure, k, axis)
 res = intensities(swt, path; energies, kernel=gaussian(fwhm=0.25))
-axisopts = (; title=L"$ϵ_T=-1$, $ϵ_Δ=-1$, $ϵ_H=+1$", titlesize=20)
+axisopts = (; title=L"$ϵ_T = %$ϵT$", titlesize=20)
 plot_intensities(res; units, axisopts, saturation=0.8, allpositive=false)
diff --git a/examples/spinw_tutorials/SW18_Distorted_kagome.jl b/examples/spinw_tutorials/SW18_Distorted_kagome.jl
index 7b7cdaf59..2a2f05690 100644
--- a/examples/spinw_tutorials/SW18_Distorted_kagome.jl
+++ b/examples/spinw_tutorials/SW18_Distorted_kagome.jl
@@ -74,17 +74,17 @@ randomize_spins!(sys2)
 minimize_energy!(sys2)
 energy_per_site(sys2)
 
-# Define a path in q-space
-
-qs = [[0,0,0], [1,0,0]]
-path = q_space_path(cryst, qs, 400)
-
-# Calculate intensities for the incommensurate spiral phase using
-# [`SpinWaveTheorySpiral`](@ref). It is necessary to provide the original `sys`,
-# consisting of a single chemical cell.
+# Return to the original system (with a single chemical cell) and construct
+# [`SpinWaveTheorySpiral`](@ref) for calculations on the incommensurate spiral
+# phase.
 
 measure = ssf_perp(sys; apply_g=false)
 swt = SpinWaveTheorySpiral(sys; measure, k, axis)
+
+# Plot intensities for a path through ``𝐪``-space.
+
+qs = [[0,0,0], [1,0,0]]
+path = q_space_path(cryst, qs, 400)
 res = intensities_bands(swt, path)
 plot_intensities(res; units)
 
diff --git a/examples/spinw_tutorials/SW19_Different_Ions.jl b/examples/spinw_tutorials/SW19_Different_Ions.jl
index 621e02213..c84b451c3 100644
--- a/examples/spinw_tutorials/SW19_Different_Ions.jl
+++ b/examples/spinw_tutorials/SW19_Different_Ions.jl
@@ -7,7 +7,7 @@
 
 using Sunny, GLMakie
 
-# Build a crystal with Cu²⁺ and Fe²⁺ ions
+# Build a crystal with Cu²⁺ and Fe²⁺ ions.
 
 units = Units(:meV, :angstrom)
 a = 3.0
@@ -19,28 +19,33 @@ types = ["Cu2", "Fe2"]
 cryst = Crystal(latvecs, positions, 1; types)
 view_crystal(cryst)
 
-# Set interactions
+# Set exchange interactions.
+
 J_Cu_Cu = 1.0
 J_Fe_Fe = 1.0
 J_Cu_Fe = -0.1
-sys = System(cryst, [1 => Moment(s=1/2, g=2), 2 => Moment(s=2, g=2)], :dipole; dims=(2, 1, 1), seed=0)
+moments = [1 => Moment(s=1/2, g=2), 2 => Moment(s=2, g=2)]
+sys = System(cryst, moments, :dipole; dims=(2, 1, 1), seed=0)
 set_exchange!(sys, J_Cu_Cu, Bond(1, 1, [-1, 0, 0]))
 set_exchange!(sys, J_Fe_Fe, Bond(2, 2, [-1, 0, 0]))
 set_exchange!(sys, J_Cu_Fe, Bond(2, 1, [0, 1, 0]))
 set_exchange!(sys, J_Cu_Fe, Bond(1, 2, [0, 0, 0]))
 
-# Find and plot a minimum energy configuration
+# Find and plot a minimum energy configuration.
 
 randomize_spins!(sys)
 minimize_energy!(sys)
 plot_spins(sys)
 
-# Configure spin wave calculation
+# Define a path through ``𝐪``-space.
 
 qs = [[0,0,0], [1,0,0]]
 path = q_space_path(cryst, qs, 400)
 
-# Plot three types of pair correlation intensities
+# Plot different pair correlation intensities by varying the
+# [`FormFactor`](@ref) on different atom types. Indices 1 and 2 refer to atoms
+# in the original chemical, and are propagated by symmetry. The special "zero"
+# form factor effectively removes the spin moment from the calculation.
 
 fig = Figure(size=(768,600))
 
@@ -63,8 +68,8 @@ fig
 
 # Calculate quantum corrections ``δS`` to spin magnitude, which arise from the
 # zero-point energy of the spin waves. The outputs are ordered following the
-# [`Site`](@ref) indexing scheme of the system: `(cell1, cell2, cell3,
-# sublattice)`, with left-side indices fastest. The two corrections ``δS ≈
+# [`Site`](@ref) indexing scheme for the system `sys`: `(cell1, cell2, cell3,
+# sublattice)`, with left-most indices fastest. The two corrections ``δS ≈
 # -0.137`` and ``δS ≈ -0.578`` apply to the Cu and Fe ions, respectively. The
 # larger correction on Fe is due to the relatively weak interchain coupling.
 
diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl
index 2e6c52679..0dc4becca 100644
--- a/src/SpinWaveTheory/SpinWaveTheory.jl
+++ b/src/SpinWaveTheory/SpinWaveTheory.jl
@@ -293,8 +293,8 @@ end
 
 # i is an "atom" index for the flattened swt.sys. However, `measure` was
 # originally constructed for a system with some crystal containing natoms. Use
-# mod1(i, natoms) to index into measure.formfactors.
+# fld1(i, natoms) to index into measure.formfactors.
 function get_swt_formfactor(measure, μ, i)
     natoms = size(measure.observables, 5)
-    measure.formfactors[μ, mod1(i, natoms)]
+    measure.formfactors[μ, fld1(i, natoms)]
 end

From f9bc23d785e2913ce0b6402e93734a30360fa978 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 09:47:02 -0600
Subject: [PATCH 11/13] Fix better

---
 src/SpinWaveTheory/SpinWaveTheory.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/SpinWaveTheory/SpinWaveTheory.jl b/src/SpinWaveTheory/SpinWaveTheory.jl
index 0dc4becca..8dda1fd84 100644
--- a/src/SpinWaveTheory/SpinWaveTheory.jl
+++ b/src/SpinWaveTheory/SpinWaveTheory.jl
@@ -291,10 +291,11 @@ function swt_data(sys::System{0}, measure)
     return SWTDataDipole(Rs, obs_localized, cs, sqrtS)
 end
 
-# i is an "atom" index for the flattened swt.sys. However, `measure` was
-# originally constructed for a system with some crystal containing natoms. Use
-# fld1(i, natoms) to index into measure.formfactors.
+# i is a site index for the flattened swt.sys. However, `measure` was originally
+# constructed for a system with nontrivial lattice dims. Use fld1(i, prod(dims))
+# to get a true atom index for the unflattened sys. This is needed to index
+# measure.formfactors.
 function get_swt_formfactor(measure, μ, i)
-    natoms = size(measure.observables, 5)
-    measure.formfactors[μ, fld1(i, natoms)]
+    sys_dims = size(measure.observables)[2:4]
+    measure.formfactors[μ, fld1(i, prod(sys_dims))]
 end

From 2c58d36e38a7723d14f810bd3122afdf291cd8b7 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 13:15:45 -0600
Subject: [PATCH 12/13] Minor tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index d6e7ebdf3..f63ce943c 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -13,28 +13,27 @@
 # Started](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia)**
 # guide. Sunny requires Julia 1.10 or later.
 #
-# From the Julia prompt, load both `Sunny` and `GLMakie`. The latter is needed
-# for graphics.
+# From the Julia prompt, load `Sunny` and also `GLMakie` for graphics.
 
 using Sunny, GLMakie
 
-# If these packages are not yet installed, Julia will offer to install them for
-# you. If executing this script gives an error, you may need to `update` Sunny
-# and GLMakie from the [built-in package
+# If these packages are not yet installed, Julia will offer to install them. If
+# executing this tutorial gives an error, you may need to update Sunny and
+# GLMakie from the [built-in package
 # manager](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia#the-built-in-julia-package-manager).
 
 # ### Units
 
 # The [`Units`](@ref Sunny.Units) object selects reference energy and length
 # scales, and uses these to provide physical constants. For example, `units.K`
-# would provide one kelvin as 0.086 meV, with the Boltzmann constant implicit.
+# returns one kelvin as 0.086 meV, where the Boltzmann constant is implicit.
 
 units = Units(:meV, :angstrom);
 
 # ### Crystal cell
 #
-# A crystallographic cell may be loaded from a `.cif` file, or can be specified
-# from atom positions and types.
+# A crystallographic cell may be loaded from a `.cif` file, or specified from
+# atom positions and types.
 #
 # Start by defining the shape of the conventional chemical cell. CoRh₂O₄ has
 # cubic spacegroup 227 (Fd-3m). Its lattice constants are 8.5 Å, and the cell

From d847a1a0768eda9459030508c84baa44a54e26c1 Mon Sep 17 00:00:00 2001
From: Kipton Barros <kbarros@gmail.com>
Date: Thu, 29 Aug 2024 18:02:40 -0600
Subject: [PATCH 13/13] Tweaks

---
 examples/01_LSWT_CoRh2O4.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/examples/01_LSWT_CoRh2O4.jl b/examples/01_LSWT_CoRh2O4.jl
index f63ce943c..e57491d18 100644
--- a/examples/01_LSWT_CoRh2O4.jl
+++ b/examples/01_LSWT_CoRh2O4.jl
@@ -13,7 +13,8 @@
 # Started](https://github.com/SunnySuite/Sunny.jl/wiki/Getting-started-with-Julia)**
 # guide. Sunny requires Julia 1.10 or later.
 #
-# From the Julia prompt, load `Sunny` and also `GLMakie` for graphics.
+# From the Julia prompt, load Sunny and also [GLMakie](https://docs.makie.org/)
+# for graphics.
 
 using Sunny, GLMakie
 
@@ -148,9 +149,9 @@ kernel = lorentzian(fwhm=0.8)
 qs = [[0, 0, 0], [1/2, 0, 0], [1/2, 1/2, 0], [0, 0, 0]]
 path = q_space_path(cryst, qs, 500)
 
-# Calculate the single-crystal scattering [`intensities`](@ref)` along the path,
-# for 300 energy points between 0 and 6 meV. Use [`plot_intensities`](@ref) to
-# visualize the result.
+# Calculate single-crystal scattering [`intensities`](@ref) along this path, for
+# energies between 0 and 6 meV. Use [`plot_intensities`](@ref) to visualize the
+# result.
 
 energies = range(0, 6, 300)
 res = intensities(swt, path; energies, kernel)