Wavevectors should be specified in absolute units.

docs/src/structure-factor.md

@@ -170,9 +170,9 @@ e.g., `bzsize=(2,2,2)`. The resulting list of wave vectors may then be passed to
 Alternatively, [`intensities_binned`](@ref) can be used to place the exact data
 into histogram bins for comparison with experiment.
 
-The convenience function [`connected_path`](@ref) returns a list of wavevectors
-sampled along a path that connects specified $𝐪$ points. This list can be used
-as an input to `intensities`. Another convenience method,
+The convenience function [`connected_path_from_rlu`](@ref) returns a list of
+wavevectors sampled along a path that connects specified $𝐪$ points. This list
+can be used as an input to `intensities`. Another convenience method,
 [`spherical_shell`](@ref) will provide a list of wave vectors on a sphere of a
 specified radius. This is useful for powder averaging. 
 

docs/src/versions.md

@@ -38,6 +38,12 @@ Remove `intensities` function. Instead, use one of
 [`intensities_interpolated`](@ref) or [`intensities_binned`](@ref). These will
 require an [`intensity_formula`](@ref), which defines a calculator (e.g., LSWT).
 
+Sunny now expects all wavevectors in units of inverse Angstrom (1/Å). This
+facilitates orientational averaging. Replace `connected_path` with
+[`connected_path_from_rlu`](@ref), which returns wavevectors in 1/Å. Now
+[`spherical_shell`](@ref) needs one fewer argument, and returns wavevectors in
+1/Å.
+
 Rename `polarize_spin!` to [`set_dipole!`](@ref) for consistency with
 [`set_coherent!`](@ref). The behavior of the former function is unchanged: the
 spin at a given site will still be polarized along the provided direction.

examples/fei2_tutorial.jl

@@ -245,27 +245,19 @@ plot_spins(sys_min; arrowlength=2.5, linewidth=0.75, arrowsize=1.5)
 
 swt = SpinWaveTheory(sys_min);
 
-# Next, select a sequence of wavevectors that will define a piecewise linear
-# path in reciprocal space. The coefficients of each $q$-vector are given in
-# reciprocal lattice units (RLU), i.e., are multiples of the recripocal lattice
-# vectors. The function [`connected_path`](@ref) will linearly sample between
-# the provided $q$-points with some `density`.
-points = [[0,   0, 0],
-          [1,   0, 0],
-          [0,   1, 0],
-          [1/2, 0, 0],
-          [0,   1, 0],
-          [0,   0, 0]] 
-labels = string.(points)
-density = 600
-path, markers = connected_path(swt, points, density);
-
-# The dispersion relation can be calculated by providing `dispersion` with a
-# `SpinWaveTheory` and a list of wave vectors. For each wave vector,
-# `dispersion` will return a list of energies, one for each band. Frequently one
-# wants dispersion information along lines that connect special wave vectors.
-# The function [`connected_path`](@ref) linearly samples between provided
-# $𝐪$-points, with a given sample density.
+# Select a sequence of wavevectors that will define a piecewise linear
+# interpolation in reciprocal lattice units (RLU).
+
+points_rlu = [[0,0,0], [1,0,0], [0,1,0], [1/2,0,0], [0,1,0], [0,0,0]]
+
+# The function [`connected_path_from_rlu`](@ref) will linearly sample between
+# the provided $q$-points with a given `density`. The `path` return value is a
+# list of wavevectors in absolute units (inverse Å). The `xticks` return value
+# keeps track of the locations of the special $𝐪$-points, and provides
+# human-readable labels for use in plotting.
+
+density = 50
+path, xticks = connected_path_from_rlu(cryst, points_rlu, density);
 
 # The [`dispersion`](@ref) function defines the quasiparticle excitation
 # energies $ω_i(𝐪)$ for each point $𝐪$ along the reciprocal space path.
@@ -288,67 +280,55 @@ disp, intensity = intensities_bands(swt, path; formula);
 # These can be plotted in GLMakie.
 
 fig = Figure()
-ax = Axis(fig[1,1]; xlabel="𝐪", ylabel="Energy (meV)",
-          xticks=(markers, labels), xticklabelrotation=π/6)
+ax = Axis(fig[1,1]; xlabel="𝐪", ylabel="Energy (meV)", xticks, xticklabelrotation=π/6)
 ylims!(ax, 0.0, 7.5)
 xlims!(ax, 1, size(disp, 1))
 for i in axes(disp)[2]
     lines!(ax, 1:length(disp[:,i]), disp[:,i]; color=intensity[:,i])
 end
 fig
 
-# To make comparisons with inelastic neutron scattering (INS) data, it is helpful to employ
-# an empirical broadening kernel.
+# To make comparisons with inelastic neutron scattering (INS) data, it is
+# helpful to employ an empirical broadening kernel, e.g., a
+# [`lorentzian`](@ref).
 
 γ = 0.15 # width in meV
 broadened_formula = intensity_formula(swt; kernel=lorentzian(γ))
 
-# The broadened intensity can be calculated for using
-# [`intensities_broadened`](@ref). One must space a path in $𝐪$-space and an
-# energy range.
+# The [`intensities_broadened`](@ref) function requires an energy range in
+# addition to the $𝐪$-space path.
 
-energies = collect(0:0.01:7.5)  # 0 < ω < 7.5 (meV).
+energies = collect(0:0.01:10)  # 0 < ω < 10 (meV).
 is1 = intensities_broadened(swt, path, energies, broadened_formula);
 
-# In real FeI$_2$, there will be competing magnetic domains associated with the
-# three possible orientations of the symmetry-broken ground state, which derive
-# from the three-fold rotational symmetry of a triangular lattice.
-#
-# The matrix $A$ defined below implements a 120° clockwise rotation on
-# wavevectors $𝐪$ defined in RLU.
+# A real FeI$_2$ sample will exhibit competing magnetic domains associated with
+# spontaneous symmetry breaking of the 6-fold rotational symmetry of the
+# triangular lattice. Note that the wavevectors $𝐪$ and $-𝐪$ are equivalent in
+# the structure factor, which leaves three distinct domain orientations, which
+# are related by 120° rotations. Rather than rotating the spin configuration
+# directly, on can rotate the $𝐪$-space path. Below, we collect and plot
+# intensity data that is averaged over all three possible orientations.
 
-A = [-1 -1 0;
-      1  0 0;
-      0  0 1]
-@assert A^3 ≈ [1 0 0; 0 1 0; 0 0 1] # Rotation by 360° is the identity
-
-# We will average the intensity data over the three symmetry-equivalent
-# rotations of the $𝐪$-space path.
-
-function rotated_intensity(R)
-    rotated_path = [R * q for q in path]
-    intensities_broadened(swt, rotated_path, energies, broadened_formula)
-end
+sθ, cθ = sincos(2π/3)
+R = [cθ -sθ 0; sθ cθ 0; 0  0 1]  # 120° rotation about the ẑ axis
 
-@assert is1 == rotated_intensity(A^0)
-is2 = rotated_intensity(A)
-is3 = rotated_intensity(A^2)
+is2 = intensities_broadened(swt, [R*q for q in path], energies, broadened_formula)
+is3 = intensities_broadened(swt, [R*R*q for q in path], energies, broadened_formula)
 is_averaged = (is1 + is2 + is3) / 3
 
 fig = Figure()
-ax = Axis(fig[1,1], xlabel="(H,0,0)", ylabel="Energy (meV)", xticks=(markers, labels),
-          )
+ax = Axis(fig[1,1]; xlabel="(H,0,0)", ylabel="Energy (meV)", xticks, xticklabelrotation=π/6)
 heatmap!(ax, 1:size(is_averaged, 1), energies, is_averaged)
 fig
 
 # This result can be directly compared to experimental neutron scattering data
-# [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1)
+# from [Bai et al.](https://doi.org/10.1038/s41567-020-01110-1)
 # ```@raw html
 # <img src="https://raw.githubusercontent.com/SunnySuite/Sunny.jl/main/docs/src/assets/FeI2_intensity.jpg">
 # ```
 #
-# (Note that the publication figure uses a different coordinate system to label
-# the wave vectors.)
+# (The publication figure accidentally used a non-standard coordinate system to
+# label the wave vectors.)
 # 
 # To get this agreement, the use of SU(3) coherent states is essential (in other
 # words, we needed a theory of multi-flavored bosons). The lower band has large

examples/longer_examples/CoRh2O4-tutorial.jl

@@ -88,7 +88,7 @@ function powder_average(sc, rs, density; η=0.1, mode=:perp, kwargs...)
     prog   = Progress(length(rs); dt=10., desc="Powder Averaging: ", color=:blue);
     output = zeros(Float64, length(rs), length(ωs(sc)))
     for (i, r) in enumerate(rs)
-        qs = spherical_shell(sc, r, density)
+        qs = spherical_shell(r; density)
         if length(qs) == 0
             qs = [[0., 0., 0.]] ## If radius is too small, just look at 0 vector
         end
@@ -108,7 +108,7 @@ end
 # ### System Definition for CoRh<sub>2</sub>O<sub>4</sub>
 
 # Define the crystal structure of CoRh$_2$O$_4$  in the conventional cell
-xtal    = Crystal(cif_path ;symprec=1e-4)
+xtal    = Crystal(cif_path; symprec=1e-4)
 magxtal = subcrystal(xtal,"Co1")
 view_crystal(magxtal,6.0)
 print_symmetry_table(magxtal, 4.0)
@@ -214,7 +214,7 @@ symQpts   = [[0.75, 0.75, 0.00],  # List of wave vectors that define a path
             [0.25, 1.00, 0.25],
             [0.00, 1.00, 0.00],
             [0.00,-4.00, 0.00]];
-(Qpts, symQind) = connected_path(sc, symQpts,densQpts);
+(Qpts, xticks) = connected_path_from_rlu(magxtal, symQpts, densQpts);
 @time  iqw = intensities(sc, Qpts, :perp; interpolation = :none, kT=integrator.kT, formfactors=ffs); 
 
 # If desired, broaden the sc in energy
@@ -227,10 +227,10 @@ iqwc  = broaden_energy(sc, iqw, (ω, ω₀) -> lorentzian(ω-ω₀, η));
 # Plot the resulting I(Q,W)    
 heatmap(1:size(iqwc, 1), ωs(sc), iqwc;
     colorrange = (0, maximum(iqwc)/20000.0),
-    axis = (
+    axis = (;
         xlabel="Momentum Transfer (r.l.u)",
         ylabel="Energy Transfer (meV)", 
-        xticks = (symQind, string.(symQpts)),
+        xticks,
         xticklabelrotation=π/5,
         aspect = 1.4,
     )

examples/longer_examples/MgCr2O4-tutorial.jl

@@ -288,7 +288,7 @@ axsqw = (xticks=-4:4, yticks=0:2:10, ylabel="E (meV)", ylabelsize=18, xlabelsize
 
 qbs = 0.0:0.5:1.5 # Determine q_b for each slice
 for (i, qb) in enumerate(qbs)
-    path, labels = connected_path(sc_pyro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40)  # Generate a path of wave
+    path, _ = connected_path_from_rlu(xtal_pyro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40)  # Generate a path of wave
                                                                          ## vectors through the BZ
     Sqω_pyro  = intensities(sc_pyro, path, :perp; kT)  # Temperature keyword enables intensity rescaling
 
@@ -303,7 +303,7 @@ fig = Figure(; resolution=(1200,900))
 
 qbs = 0.0:0.5:1.5
 for (i, qb) in enumerate(qbs)
-    path, labels = connected_path(sc_mgcro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40)  # Generate a path of wave
+    path, _ = connected_path_from_rlu(xtal_mgcro, [[-4.0, qb, 0.0],[4.0, qb, 0.0]], 40)  # Generate a path of wave
                                                                          ## vectors through the BZ
     Sqω_mgcro  = intensities(sc_mgcro, path, :perp; kT)  # Temperature keyword enables intensity rescaling