Fix intensity scale for SampledCorrelations with reshaped system

examples/spinw_tutorials/SW05_Simple_kagome_FM.jl

@@ -69,6 +69,6 @@ res2 = powder_average(cryst, radii, 1000) do qs
 end
 
 fig = Figure(size=(768, 800))
-plot_intensities!(fig[1, 1], res1; units, colorrange=(0,10))
-plot_intensities!(fig[2, 1], res2; units, colorrange=(0,10))
+plot_intensities!(fig[1, 1], res1; units, colorrange=(0,10), title="FWHM 0.02 meV")
+plot_intensities!(fig[2, 1], res2; units, colorrange=(0,10), title="FWHM 0.25 meV")
 fig

examples/spinw_tutorials/SW18_Distorted_kagome.jl

@@ -24,8 +24,8 @@ view_crystal(cryst)
 
 # Define the interactions.
 
-spininfos = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
-sys = System(cryst, spininfos, :dipole, seed=0)
+moments = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
+sys = System(cryst, moments, :dipole, seed=0)
 J   = -2
 Jp  = -1
 Jab = 0.75

ext/PlottingExt.jl

@@ -609,9 +609,8 @@ function draw_atoms_or_dipoles(; ax, full_crystal_toggle, dipole_menu, cryst, sy
             # Show dipoles. Mostly consistent with code in plot_spins.
             if !isnothing(sys) && ismagnetic
                 sites = Sunny.position_to_site.(Ref(sys), rs)
-                L = length(eachsite(sys))
-                g0 = norm(sys.gs) / sqrt(L * 3)
-                N0 = norm(sys.Ns) / sqrt(L)
+                g0 = norm(sys.gs) / sqrt(length(sys.gs) * 3)
+                N0 = norm(sys.Ns) / sqrt(length(sys.Ns))
                 s0 = (N0 - 1) / 2
                 spin_dipoles = sys.dipoles[sites] / s0
                 magn_dipoles = magnetic_moment.(Ref(sys), sites) / (s0*g0)

src/Integrators.jl

@@ -210,7 +210,7 @@ function suggest_timestep_aux(sys::System{N}, integrator; tol) where N
     # information about the frequency of oscillation. Consider, e.g., a spin
     # approximately aligned with an external field: the precession frequency is
     # given by |∇E| = |B|.
-    drift_rms = sqrt(acc/length(eachsite(sys)))
+    drift_rms = sqrt(acc/nsites(sys))
     if iszero(drift_rms)
         error("Cannot suggest a timestep without an energy scale!")
     end

src/KPM/SpinWaveTheoryKPM.jl

@@ -70,7 +70,7 @@ function intensities!(data, swt_kpm::SpinWaveTheoryKPM, qpts; energies, kernel::
     @assert eltype(data) == eltype(measure)
     @assert size(data) == (length(energies), length(qpts.qs))
 
-    Na = length(eachsite(sys))
+    Na = nsites(sys)
     Ncells = Na / natoms(cryst)
     Nf = nflavors(swt)
     L = Nf*Na

src/Measurements/IntensitiesTypes.jl

@@ -43,7 +43,12 @@ struct PowderIntensities{T} <: AbstractIntensities
 end
 
 function Base.show(io::IO, res::AbstractIntensities)
-    sz = join([size(res.data, 1), sizestr(res.qpts)], "×")
+    sz = string(size(res.data, 1)) * "×" * sizestr(res.qpts)
+    print(io, string(typeof(res)) * " ($sz elements)")
+end
+
+function Base.show(io::IO, res::StaticIntensities)
+    sz = sizestr(res.qpts)
     print(io, string(typeof(res)) * " ($sz elements)")
 end
 

src/MonteCarlo/WangLandau.jl

@@ -31,8 +31,7 @@ function drive_system_to_energy_bounds!(sys::System{N}, bounds, propose::PR, max
     sampler = LocalSampler(; kT, nsweeps=0.01, propose)
     for _ in 1:max_sweeps
         step!(sys, sampler)
-        nsites = length(eachsite(sys))
-        E = E0 + sampler.ΔE/nsites
+        E = E0 + sampler.ΔE/nsites(sys)
         inbounds(E, bounds) && return
     end
 
@@ -64,17 +63,17 @@ end
 # sample for specified number of sweeps
 function step_ensemble!(WL::WangLandau, nsweeps::Int64)
     (; sys) = WL
-    nsites = length(eachsite(sys))
-    WL.E = energy(sys) / nsites
+    N = nsites(sys)
+    WL.E = energy_per_site(sys)
     @assert inbounds(WL.E, WL.bounds)
 
     # try to fulfill the histogram criterion in set number of MC sweeps
     accepts = 0
-    for _ in 1:nsweeps*nsites
+    for _ in 1:nsweeps*N
         # perform single-spin update and calculate proposal energy
         site = rand(sys.rng, eachsite(sys))
         state = WL.propose(sys, site)
-        E_next = WL.E + local_energy_change(sys, site, state) / nsites
+        E_next = WL.E + local_energy_change(sys, site, state) / N
 
         if inbounds(E_next, WL.bounds)
             iszero(WL.ln_g[E_next]) && add_new!(WL, E_next)

src/Operators/Rotation.jl

@@ -118,7 +118,8 @@ function unitary_for_rotation(R::Mat3, gen)
     return exp(-im*θ*(n'*gen))
 end
 
-# Unitary for a rotation matrix in the N-dimensional irrep of SU(2).
+# Unitary for a rotation matrix in the N-dimensional irrep of SU(2). TODO: Use
+# polynomial expansions for speed: http://www.emis.de/journals/SIGMA/2014/084/.
 function unitary_irrep_for_rotation(R::Mat3; N::Int)
     gen = spin_matrices_of_dim(; N)
     unitary_for_rotation(R, gen)

src/SampledCorrelations/CorrelationUtils.jl

@@ -1,13 +1,15 @@
 """
-    available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1))
+    available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1))
 
-Returns all wave vectors for which `sc` contains exact values. `bsize` specifies
-the number of Brillouin zones to be included.
+Returns the grid of wave vectors for which `sc` contains exact values.
+Optionally extend by a given number of `counts` along each grid axis. If the
+system was not reshaped, then the number of Brillouin zones included is
+`prod(counts)`.
 """
-function available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1))
+function available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1))
     Ls = sc.sys_dims
     offsets = map(L -> isodd(L) ? 1 : 0, Ls)
-    up = Ls .* bzsize
+    up = Ls .* counts
     hi = map(L -> L - div(L, 2), up) .- offsets
     lo = map(L -> 1 - div(L, 2), up) .- offsets
 

src/SampledCorrelations/DataRetrieval.jl

@@ -70,7 +70,8 @@ end
 
 contains_dynamic_correlations(sc) = !isnan(sc.Δω)
 
-# Documented under intensities function for LSWT.
+# Documented under intensities function for LSWT. TODO: As a hack, this function
+# is also being used as the back-end to intensities_static.
 function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT)
     if !isnothing(kT) && kT <= 0
         error("Positive `kT` required for classical-to-quantum corrections, or set `kT=nothing` to disable.")
@@ -92,16 +93,21 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
 
     # Prepare memory and configuration variables for actual calculation
     qpts = Base.convert(AbstractQPoints, qpts)
+    qs_reshaped = [to_reshaped_rlu(sc, q) for q in qpts.qs]
+
     ffs = sc.measure.formfactors[1, :] # FIXME
     intensities = zeros(eltype(sc.measure), isnan(sc.Δω) ? 1 : length(ωs), length(qpts.qs)) # N.B.: Inefficient indexing order to mimic LSWT
-    q_idx_info = pruned_wave_vector_info(sc, qpts.qs)
+    q_idx_info = pruned_wave_vector_info(sc, qs_reshaped)
     crystal = @something sc.origin_crystal sc.crystal
     NCorr  = Val{size(sc.data, 1)}()
     NAtoms = Val{size(sc.data, 2)}()
 
     # Intensities calculation
     intensities_rounded!(intensities, sc.data, sc.crystal, sc.measure, ffs, q_idx_info, ωidcs, NCorr, NAtoms)
 
+    # Convert to a q-space density in original (not reshaped) RLU.
+    intensities .*= det(sc.crystal.recipvecs) / det(crystal.recipvecs)
+
     # Post-processing steps for dynamical correlations 
     if contains_dynamic_correlations(sc) 
         # Convert time axis to a density.
@@ -120,7 +126,6 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
 
     intensities = reshape(intensities, length(ωs), size(qpts.qs)...)
 
-    # TODO: Refactor this logic so that return value is uniform.
     return if contains_dynamic_correlations(sc) 
         Intensities(crystal, qpts, collect(ωs), intensities)
     else
@@ -202,6 +207,3 @@ end
 =#
 
 
-
-
-

src/SampledCorrelations/SampledCorrelations.jl

@@ -3,7 +3,7 @@ mutable struct SampledCorrelations
     const data           :: Array{ComplexF64, 7}                 # Raw SF with sublattice indices (ncorrs × natoms × natoms × sys_dims × nω)
     const M              :: Union{Nothing, Array{Float64, 7}}    # Running estimate of (nsamples - 1)*σ² (where σ² is the variance of intensities)
     const crystal        :: Crystal                              # Crystal for interpretation of q indices in `data`
-    const origin_crystal :: Union{Nothing,Crystal}               # Original user-specified crystal (if different from above) -- needed for FormFactor accounting
+    const origin_crystal :: Union{Nothing,Crystal}               # Original user-specified crystal (if different from above)
     const Δω             :: Float64                              # Energy step size (could make this a virtual property)  
     measure              :: MeasureSpec                          # Observable, correlation pairs, and combiner
 

src/SpinWaveTheory/DispersionAndIntensities.jl

@@ -154,7 +154,7 @@ function intensities_bands(swt::SpinWaveTheory, qpts; kT=0, with_negative=false)
 
     # Number of (magnetic) atoms in magnetic cell
     @assert sys.dims == (1,1,1)
-    Na = length(eachsite(sys))
+    Na = nsites(sys)
     # Number of chemical cells in magnetic cell
     Ncells = Na / natoms(cryst)
     # Number of quasiparticle modes

src/SpinWaveTheory/SpinWaveTheory.jl

@@ -43,7 +43,7 @@ function SpinWaveTheory(sys::System; measure::Union{Nothing, MeasureSpec}, regul
     end
 
     measure = @something measure empty_measurespec(sys)
-    if length(eachsite(sys)) != prod(size(measure.observables)[2:5])
+    if nsites(sys) != prod(size(measure.observables)[2:5])
         error("Size mismatch. Check that measure is built using consistent system.")
     end
 
@@ -162,7 +162,7 @@ end
 # rotating these into the local reference frame determined by the ground state.
 function swt_data(sys::System{N}, measure) where N
     # Calculate transformation matrices into local reference frames
-    Na = length(eachsite(sys))
+    Na = nsites(sys)
     Nobs = size(measure.observables, 1)
     observables = reshape(measure.observables, Nobs, Na)
 
@@ -228,7 +228,7 @@ end
 
 # Compute Stevens coefficients in the local reference frame
 function swt_data(sys::System{0}, measure)
-    Na = length(eachsite(sys))
+    Na = nsites(sys)
     Nobs = size(measure.observables, 1)
 
     # Operators for rotating vectors into local frame

src/Spiral/SpinWaveTheorySpiral.jl

@@ -189,9 +189,9 @@ function intensities_bands(sswt::SpinWaveTheorySpiral, qpts; kT=0) # TODO: branc
 
     # Number of atoms in magnetic cell
     @assert sys.dims == (1,1,1)
-    Na = length(eachsite(sys))
+    Na = nsites(sys)
     # Number of chemical cells in magnetic cell
-    Ncells = Na / natoms(cryst) # TODO check invariance
+    Ncells = Na / natoms(cryst)
     # Number of quasiparticle modes
     L = nbands(swt)
     # Number of wavevectors

src/Spiral/SpiralEnergy.jl

@@ -33,7 +33,7 @@ The [`spiral_energy`](@ref) divided by the number of sites in `sys`. The special
 case ``𝐤 = 0`` yields a result identical to [`energy_per_site`](@ref).
 """
 function spiral_energy_per_site(sys::System{0}; k, axis)
-    return spiral_energy(sys; k, axis) / length(eachsite(sys))
+    return spiral_energy(sys; k, axis) / nsites(sys)
 end
 
 function spiral_energy_and_gradient_aux!(dEds, sys::System{0}; k, axis)

src/Sunny.jl

@@ -126,7 +126,7 @@ export BinningParameters, load_nxs, generate_mantid_script_from_binning_paramete
 include("deprecated.jl")
 export set_external_field!, set_external_field_at!, meV_per_K,
     dynamic_correlations, instant_correlations, intensity_formula, reciprocal_space_path,
-    set_spiral_order_on_sublattice!, set_spiral_order!
+    delta_function_kernel, set_spiral_order_on_sublattice!, set_spiral_order!
 
 isloaded(pkg::String) = any(k -> k.name == pkg, keys(Base.loaded_modules))
 

src/System/Interactions.jl

@@ -234,7 +234,7 @@ end
 The total system [`energy`](@ref) divided by the number of sites.
 """
 function energy_per_site(sys::System{N}) where N
-    return energy(sys) / length(eachsite(sys))
+    return energy(sys) / nsites(sys)
 end
 
 """

src/System/System.jl

@@ -194,7 +194,7 @@ const Site = Union{NTuple{4, Int}, CartesianIndex{4}}
 end
 
 # Offset a `cell` by `ncells`
-@inline offsetc(cell::CartesianIndex{3}, ncells, dims) = CartesianIndex(altmod1.(Tuple(cell) .+ Tuple(ncells), dims))
+@inline offsetc(cell::CartesianIndex{3}, delta, dims) = CartesianIndex(altmod1.(Tuple(cell) .+ Tuple(delta), dims))
 
 # Split a site `site` into its cell and sublattice parts
 @inline to_cell(site) = CartesianIndex((site[1],site[2],site[3]))
@@ -227,6 +227,15 @@ An iterator over all [`Site`](@ref)s in the system.
 """
 @inline eachsite(sys::System) = CartesianIndices(sys.dipoles)
 
+"""
+nsites(sys::System) = length(eachsite(sys))
+"""
+nsites(sys::System) = length(eachsite(sys))
+
+# Number of (original) crystal cells in the system
+ncells(sys::System) = nsites(sys) / natoms(orig_crystal(sys))
+
+
 """
     global_position(sys::System, site::Site)
 

src/deprecated.jl

@@ -85,7 +85,7 @@ Base.@deprecate instant_correlations(sys; opts...) let
 end
 
 Base.@deprecate intensity_formula(opts1...; opts2...) let
-    error("Formulas have been removed. Call `intensities` directly.")
+    error("Formulas have been removed. See revised Sunny docs for new interface.")
 end
 
 Base.@deprecate reciprocal_space_path(cryst::Crystal, qs, density) let
@@ -110,6 +110,9 @@ Base.@deprecate SpinInfo(i; S, g) let
     i => Moment(; s=S, g)
 end
 
+Base.@deprecate_binding delta_function_kernel nothing
+
+
 # REMEMBER TO ALSO DELETE:
 #
 # * view_crystal(cryst, max_dist)