Update docs and release notes for 0.8.0 (#463)

Author: kbarros

docs/src/versions.md

@@ -1,7 +1,7 @@
 # Release Notes
 
 ## v0.8.0
-(In development)
+(Nov 30, 2025)
 
 **Breaking changes** in this release.
 
@@ -33,7 +33,7 @@
 * Add missing 3D support for [`q_space_grid`](@ref) ([#457](@ref)).
 * If user-provided lattice vectors do not match crystallographic conventions,
   suggest the use of [`standardize`](@ref) ([#461](@ref)).
-* Fixes to [`SpinWaveTheoryKPM`](@ref) ([#462](@ref)).
+* Improve robustness of [`SpinWaveTheoryKPM`](@ref) ([#462](@ref)).
 
 ## v0.7.8
 (Jul 1, 2025)

ext/PlottingExt/PlotIntensities.jl

@@ -338,9 +338,9 @@ is:
 
     To export the figure to a PDF file, use the CairoMakie backend and call
 
-    ```jl
-    fig = plot_intensities(..., interpolate=true)
-    save("myfile.pdf", fig)
+    ```julia
+        fig = plot_intensities(..., interpolate=true)
+        save("myfile.pdf", fig)
     ```
 
     The choice `interpolate=true` significantly reduces file sizes for plots

src/Binning/Binning.jl

@@ -25,7 +25,7 @@ The convention for the binning scheme is that:
   
 A `value` can be binned by computing its bin index:
 
-```jl
+```julia
     coords = covectors * value
     bin_ix = 1 .+ floor.(Int64, (coords .- binstart) ./ binwidth)
 ```

src/KPM/SpinWaveTheoryKPM.jl

@@ -1,34 +1,43 @@
 """
     SpinWaveTheoryKPM(sys::System; measure, regularization=1e-8, tol=nothing,
-                      niters=nothing, niters_bounds=10, method=:lanczos)
+                      niters=nothing, method=:lanczos)
 
 A variant of [`SpinWaveTheory`](@ref) that estimates [`intensities`](@ref) using
 iterated matrix-vector products. By avoiding direct matrix diagonalization, this
 method reduces computational cost from cubic to linear-scaling in the system
 size ``N``. Large system sizes can arise, e.g., for models of quenched disorder
 or models with nearly incommensurate ordering wavevectors.
 
-Computational cost scales like ``𝒪(N M + M^2)``. The number of iterations ``M``
-may be specified directly with the `niters` parameter. Alternatively, one may
-specify a dimensionless error tolerance `tol` such that `M ≈ -2 log10(tol) Δϵ /
-fwhm` where `Δϵ` is the estimated spectral bandwidth of excitations and `fwhm`
-is the full width at half maximum of the user-supplied broadening `kernel`.
-Common choices for `tol` are `0.05` (more speed) or `0.01` (more accuracy).
-Either `tol` or `niters` is required. The parameter `niters_bounds` selects the
-Krylov subspace dimension to be used for estimating spectral bounds.
-
-!!! warning "Accuracy considerations"  
-    Energy-space resolution scales inversely with the number ``M`` of
-    iterations. Such broadening errors can usually be well controlled via the
-    tolerance parameter `tol`. A more serious problem is intensity loss at bands
-    with small excitation energy, e.g., in the viscinity of Goldstone modes. In
-    certain cases this missing intensity will be due to numerical roundoff error
-    that cannot be fixed by increasing the number of iterations ``M``.
-
-Two `method` options are available, `:lanczos` and `:kpm`.  Lanczos is generally
-preferred, as it achieves near-optimal accuracy for a given number of iterations
-[1, 2]. The Lanczos implementation builds on earlier research using the Kernel
-Polynomial Method [3], which may be of historical interest.
+Energy resolution is controlled by the dimensionless `tol` parameter. Common
+choices are `tol=0.05` (more speed) or `0.01` (more accuracy). This will
+determine the number of iterations as `M ≈ -2 log10(tol) Δϵ / fwhm`, where `Δϵ`
+is the estimated spectral bandwidth of excitations and `fwhm` is the full width
+at half maximum of the user-supplied broadening `kernel`. Computational cost
+scales like ``𝒪(N M + M^2)``. Use `niters` instead of `tol` to directly specify
+``M``.
+
+!!! warning "Intensity loss at low-energy excitations"
+
+    Not all numerical artifacts can be resolved by reducing `tol`. In particular,
+    there may be unavoidable intensity loss at low-energy excitations, e.g., near
+    Goldstone modes. This type of error originates from finite numerical precision
+    and ill-conditioning of the dynamical matrix.
+
+!!! tip "Consider `SampledCorrelations` when calculating powder averages"
+
+    Spin wave theory requires an independent calculation for each ``𝐪`` point of
+    interest. Consequently, it can be very slow to sample a 3D volume of
+    ``𝐪``-space, e.g., as required for a [`powder_average`](@ref). A compelling
+    alternative may be [`SampledCorrelations`](@ref). It uses real-time spin
+    dynamics to calculate structure factor data over the entire 3D grid of
+    commensurate ``𝐪``-vectors in one shot. This may provide a considerable speedup
+    at the cost of: limited ``𝐪``-space resolution and stochastic error due to
+    statistical sampling.
+
+Two choices of `method` are possible. Lanczos is the default because it achieves
+appears near-optimal accuracy at fixed iterations ``M`` [1, 2] and can detect
+energetic instabilities. The alternative, `method=:kpm`, implements the Kernel
+Polynomial Method as described in Ref. [3], which may be of historical interest.
 
 ## References
 

test/test_lswt.jl

@@ -613,7 +613,7 @@ end
     res = intensities(swt, qs; energies, kernel)
     # println(round.(res.data; digits=12))
     data_ref = [0.042551644188 0.148531187027 0.042551644188; 0.123710785142 0.944407715529 0.123710785142; 1.079575108835 1.261286085876 1.079575108835; 0.492703854423 0.168505531387 0.492703854423; 0.087770506619 0.06066521855 0.087770506619; 0.034461978527 0.030825188879 0.034461978527; 0.018214262744 0.018585357642 0.018214262744; 0.011230027228 0.012410289203 0.011230027228; 0.007607320239 0.008868510563 0.007607320239; 0.005490942587 0.006651305827 0.005490942587; 0.004148615687 0.005172181047 0.004148615687]
-    isapprox(res.data, data_ref; atol=1e-9)
+    @test isapprox(res.data, data_ref; atol=1e-9)
 end
 
 

test/test_scga.jl

@@ -1,7 +1,23 @@
-@testitem "diamond_lattice" begin
-    # test against JuliaSCGA (S. Gao)
-    res_jscga = [0.7046602277469309, 0.8230846832863896, 0.23309034250417973, 0.40975668535137943, 0.8474163786642979,
-                 0.8230846832694241, 0.723491683211756, 0.5939752161027589, 0.6506966347286152, 0.8012263819500781]
+@testitem "Square lattice" begin
+    latvecs = lattice_vectors(1, 1, 10, 90, 90, 90)
+    cryst = Crystal(latvecs, [[0, 0, 0]])
+    s = 1
+    sys = System(cryst, [1 => Moment(; s, g=1)], :dipole; seed=0)
+    set_exchange!(sys, -1, Bond(1, 1, [1, 0, 0]))
+    set_exchange!(sys, 0.5, Bond(1, 1, [1, 1, 0]))
+    set_exchange!(sys, 0.25, Bond(1, 1, [2, 0, 0]))
+    measure = ssf_perp(sys)
+    kT = 27.5*meV_per_K
+    scga = SCGA(sys; measure, kT, dq=1/8)
+    path = q_space_path(cryst, [[-1, -1, 0], [-0.04, -1, 0]], 9)
+    res = Sunny.intensities_static(scga, path)
+    ref_jscga = [0.3578208506624103, 0.3850668587212228, 0.4211633125595451, 0.3930558742921638, 0.3586715762816389,
+                 0.3775147329089877, 0.4188431376589759, 0.4009835864744404, 0.36119385315852837]
+    ref_jscga *= s^2 * Sunny.natoms(cryst)
+    @test isapprox(vec(res.data)/2, ref_jscga; rtol=1e-5)
+end
+
+@testitem "Diamond lattice" begin
     a = 8.5031 # (Å)
     latvecs = lattice_vectors(a, a, a, 90, 90, 90)
     cryst = Crystal(latvecs, [[0, 0, 0]], 227; choice="1")
@@ -13,34 +29,18 @@
     measure = ssf_perp(sys)
     kT = 15*meV_per_K
     scga = SCGA(sys; measure, kT, dq=1/8)
-    γ = s^2 * Sunny.natoms(cryst)
     grid = q_space_grid(cryst, [1, 0, 0], range(0, 3.6, 5), [0, 1, 0], range(0, 0.9, 2))
     res = intensities_static(scga, grid)
-    @test isapprox(vec(res.data)/γ, res_jscga; rtol=1e-8)
-end
-
-@testitem "square_lattice" begin
-    # test against JuliaSCGA (S. Gao)
-    res_jscga = [0.3578208506624103, 0.3850668587212228, 0.4211633125595451, 0.3930558742921638, 0.3586715762816389,
-                 0.3775147329089877, 0.4188431376589759, 0.4009835864744404, 0.36119385315852837]
-    a = 1 # (Å)
-    latvecs = lattice_vectors(a, a, 10a, 90, 90, 90)
-    cryst = Crystal(latvecs, [[0, 0, 0]])
-    sys = System(cryst, [1 => Moment(; s=1, g=1)], :dipole; seed=0)
-    set_exchange!(sys, -1, Bond(1, 1, [1, 0, 0]))
-    set_exchange!(sys, 0.5, Bond(1, 1, [1, 1, 0]))
-    set_exchange!(sys, 0.25, Bond(1, 1, [2, 0, 0]))
-    measure = ssf_perp(sys)
-    kT = 27.5*meV_per_K
-    scga = SCGA(sys; measure, kT, dq=1/8)
-    path = q_space_path(cryst, [[-1, -1, 0], [-0.04, -1, 0]], 9)
-    res = Sunny.intensities_static(scga, path)
-    @test isapprox(vec(res.data)/2, res_jscga; rtol=1e-5)
+    ref_jscga = [0.7046602277469309, 0.8230846832863896, 0.23309034250417973, 0.40975668535137943, 0.8474163786642979,
+                 0.8230846832694241, 0.723491683211756, 0.5939752161027589, 0.6506966347286152, 0.8012263819500781]
+    ref_jscga *= s^2 * Sunny.natoms(cryst)
+    @test isapprox(vec(res.data), ref_jscga; rtol=1e-8)
 end
 
 @testitem "MgCr2O4" begin
     using LinearAlgebra
-    # Reproduce calculation in Conlon and Chalker, PRL 102, 237206 (2009)
+    # Reproduce calculation in Bai et al., Phys. Rev. Lett. 122, 097201 (2019).
+    # See also Conlon and Chalker, Phys. Rev. B 81, 224413 (2010).
     latvecs = lattice_vectors(8.3342, 8.3342, 8.3342, 90, 90, 90)
     positions = [[1/2, 1/2, 1/2]]
     cryst = Crystal(latvecs, positions, 227)
@@ -59,18 +59,15 @@ end
     scga = SCGA(sys; measure, kT, dq=1/4)
     grid = q_space_grid(cryst, [1, 0, 0], range(-1.5, 1.5, 4), [0, 1, 0], range(-1.5, 1.5, 4))
     res = Sunny.intensities_static(scga, grid)
-    res_cc = [2.4168819 1.237733 1.237733 2.4168819;
-              1.237733 0.16242284 0.16242284 1.237733;
-              1.237733 0.16242284 0.16242284 1.237733;
-              2.4168819 1.237733 1.237733 2.4168819]
-    γ = s*(s+1) / det(primitive_cell(cryst)) # Adapt to normalization convention of C+C
-    @test isapprox(vec(res.data), γ*vec(res_cc); rtol=1e-3)
+    ref = [2.4168819 1.237733 1.237733 2.4168819;
+           1.237733 0.16242284 0.16242284 1.237733;
+           1.237733 0.16242284 0.16242284 1.237733;
+           2.4168819 1.237733 1.237733 2.4168819]
+    ref *= s*(s+1) / det(primitive_cell(cryst)) # Differing normalization convention
+    @test isapprox(vec(res.data), vec(ref); rtol=1e-3)
 end
 
 @testitem "Arbitrary Anisotropy" begin
-    # test against JuliaSCGA (S. Gao)
-    res_jscga = [0.700944539713289, 0.6175127864675868, 0.5205893697530085, 0.48172879530096047,
-                 0.5219040226511135, 0.6218544838522482, 0.7110417061581527, 0.738269833048121]
     latvecs = lattice_vectors(1, 1, 10, 90, 90, 90)
     cryst = Crystal(latvecs, [[0, 0 ,0]], 1)
     sys = System(cryst, [1 => Moment(; s=1, g=1)], :dipole; seed=0)
@@ -91,7 +88,9 @@ end
     scga = SCGA(sys; measure, kT, dq=1/8)
     path = q_space_path(cryst, [[0.125, 0.625, 0], [1, 0.625, 0]], 8)
     res = Sunny.intensities_static(scga, path)
-    @test isapprox(vec(res.data), res_jscga; rtol=1e-6)
+    ref_jscga = [0.700944539713289, 0.6175127864675868, 0.5205893697530085, 0.48172879530096047,
+                 0.5219040226511135, 0.6218544838522482, 0.7110417061581527, 0.738269833048121]
+    @test isapprox(vec(res.data), ref_jscga; rtol=1e-6)
 end
 
 @testitem "Ferrimagnetic chain sum rule" begin