Add intensities_static with integration bounds (#298)

* intensities_static with integration bounds for SWT

* add integration bounds to intensities_static for SampledCorrelations

---------

Co-authored-by: ddahlbom <david.dahlbom@gmail.com>

docs/src/structure-factor.md

@@ -450,7 +450,7 @@ the classical Boltzmann distribution.
 
 ### The instantaneous structure factor
 
-Use [`intensities_instant`](@ref) to calculate ``\mathcal{S}(𝐪)``, i.e.,
+Use [`intensities_static`](@ref) to calculate ``\mathcal{S}(𝐪)``, i.e.,
  correlations that are "instantaneous" in real-time. Mathematically,
 ``\mathcal{S}(𝐪)`` denotes an integral of the full dynamical structure factor
 ``\mathcal{S}(𝐪, ω)``, taken over all energies ``ω``. In
@@ -460,7 +460,7 @@ be applied within [`intensities`](@ref) prior to energy integration.
 
 Sunny also supports a mechanism to calculate static correlations without any
 spin dynamics. To collect such statistics, construct a
-`SampledCorrelationsStatic` object. In this case, `intensities_instant` will
+`SampledCorrelationsStatic` object. In this case, `intensities_static` will
 return static correlations sampled from the classical Boltzmann distribution.
-This dynamics-free approach is faster, but may miss important quantum mechanical
-features of the structure factors.
+This dynamics-free approach is faster, but may miss important features that
+derive from the quantum mechanical excitation spectrum.

docs/src/versions.md

@@ -7,7 +7,7 @@ This **major release** introduces several breaking changes:
 
 * The interface for calculating intensities has been revised to unify
   functionality across backends. The functions [`intensities_bands`](@ref),
-  [`intensities`](@ref), and [`intensities_instant`](@ref) no longer expect a
+  [`intensities`](@ref), and [`intensities_static`](@ref) no longer expect a
   "formula", and instead take keyword arguments directly. Pair correlations are
   now specified using [`ssf_perp`](@ref) and related functions. The constructors
   [`SampledCorrelations`](@ref) and [`SampledCorrelationsStatic`](@ref) replace

examples/01_LSWT_SU3_FeI2.jl

@@ -194,7 +194,7 @@ plot_spins(sys; color=[s[3] for s in sys.dipoles])
 # A different understanding of the magnetic ordering can be obtained by moving
 # to Fourier space. The 'instantaneous' structure factor ``𝒮(𝐪)`` is an
 # experimental observable. To investigate ``𝒮(𝐪)`` as true 3D data, Sunny
-# provides [`intensities_instant`](@ref). Here, however, we will use
+# provides [`SampledCorrelationsStatic`](@ref). Here, however, we will use
 # [`print_wrapped_intensities`](@ref), which gives average intensities for the
 # individual Bravais sublattices (in effect, all wavevectors are wrapped to the
 # first Brillouin zone).

examples/03_LLD_CoRh2O4.jl

@@ -99,7 +99,7 @@ grid = q_space_grid(cryst, [1, 0, 0], range(-10, 10, 200), [0, 1, 0], (-10, 10))
 # [`FormFactor`](@ref) for Co2⁺.
 
 formfactors = [FormFactor("Co2")]
-res = intensities_instant(sc, grid; formfactors)
+res = intensities_static(sc, grid; formfactors)
 plot_intensities(res; saturation=0.995)
 
 
@@ -142,7 +142,7 @@ qpts = q_space_path(cryst, qs, 500)
 
 # Calculate ``I(𝐪, ω)`` intensities along this path and plot.
 
-res = intensities(sc, qpts; energies, kT)
+res = intensities(sc, qpts; energies, formfactors, kT)
 plot_intensities(res; units, saturation=0.85, colormap=:viridis)
 
 # ### Powder averaged intensity

examples/spinw_tutorials/SW10_Energy_cut.jl

@@ -37,18 +37,8 @@ swt = SpinWaveTheory(sys; measure=ssf_perp(sys))
 res = intensities(swt, grid; energies=[3.75], kernel=gaussian(fwhm=0.2))
 plot_intensities(res; units)
 
-# Sunny does not include Horace integration. For now, one can write custom code
-# to manually integrate intensity bands between 3.5 and 4 meV.
+# Integrate intensities between 3.5 and 4 meV using [`intensities_static`](@ref)
+# with the `bounds` option.
 
-res = intensities_bands(swt, grid)
-summed_intensity = zeros(size(res.disp)[2:3])
-for ci in CartesianIndices(res.disp)
-    if 3.5 < res.disp[ci] < 4.0
-        summed_intensity[ci[2], ci[3]] += res.data[ci]
-    end
-end
-
-# To make the plot look nice, for now we use a private function internal to
-# Sunny.
-
-plot_intensities(Sunny.InstantIntensities(cryst, grid, summed_intensity))
+res = intensities_static(swt, grid; bounds=(3.5, 4.0))
+plot_intensities(res; units)

examples/spinw_tutorials/SW11_La2CuO4.jl

@@ -51,5 +51,5 @@ plot_intensities(res; units)
 
 # Plot instantaneous itensities, integrated over ω.
 
-res = intensities_instant(swt, path)
+res = intensities_static(swt, path)
 plot_intensities(res; colorrange=(0,20), units)

ext/PlottingExt.jl

@@ -1127,7 +1127,7 @@ function Sunny.plot_intensities!(panel, res::Sunny.Intensities{Float64}; colorma
     end
 end
 
-function Sunny.plot_intensities!(panel, res::Sunny.InstantIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, axisopts=Dict())
+function Sunny.plot_intensities!(panel, res::Sunny.StaticIntensities{Float64}; colormap=nothing, colorrange=nothing, saturation=0.9, allpositive=true, units=nothing, into=nothing, axisopts=Dict())
     (; crystal, qpts, data) = res
 
     colorrange_suggest = colorrange_from_data(; data=reshape(data, 1, size(data)...), saturation, sensitivity=0, allpositive)

src/MagneticOrdering.jl

@@ -13,7 +13,7 @@ this function will typically be used as input to
 
 Because this function does not incorporate phase information in its averaging
 over sublattices, the printed weights are not directly comparable with
-experiment. For that purpose, use [`intensities_instant`](@ref) instead.
+experiment. For that purpose, use [`SampledCorrelationsStatic`](@ref) instead.
 """
 function print_wrapped_intensities(sys::System{N}; nmax=10) where N
     sys.crystal == orig_crystal(sys) || error("Cannot perform this analysis on reshaped system.")

src/Measurements/IntensitiesTypes.jl

@@ -22,7 +22,7 @@ struct Intensities{T, Q <: AbstractQPoints, D} <: AbstractIntensities
     data :: Array{T, D} # (nω × nq...)
 end
 
-struct InstantIntensities{T, Q <: AbstractQPoints, D} <: AbstractIntensities
+struct StaticIntensities{T, Q <: AbstractQPoints, D} <: AbstractIntensities
     # Original chemical cell
     crystal :: Crystal
     # Wavevectors in RLU

src/SampledCorrelations/DataRetrieval.jl

@@ -124,7 +124,7 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, fo
     return if contains_dynamic_correlations(sc) 
         Intensities(crystal, qpts, collect(ωs), intensities)
     else
-        InstantIntensities(crystal, qpts, dropdims(intensities; dims=1))
+        StaticIntensities(crystal, qpts, dropdims(intensities; dims=1))
     end
 end
 
@@ -154,14 +154,20 @@ function intensities_rounded!(intensities, data, crystal, measure::MeasureSpec{O
 end
 
 
-function intensities_instant(sc::SampledCorrelations, qpts; formfactors=nothing, kT)
-    res = intensities(sc, qpts; formfactors, kT, energies=:available_with_negative)
+function intensities_static(sc::SampledCorrelations, qpts; bounds = (-Inf, Inf), formfactors=nothing, kT)
+    ωs = available_energies(sc; negative_energies=true)
+    ωidcs = findall(x -> bounds[1] <= x < bounds[2], ωs)
+    if iszero(length(ωidcs))
+        error("No information available within specified energy `bounds`. Try a larger interval.")
+    end
+    energies = sort(ωs[ωidcs])
+    res = intensities(sc, qpts; formfactors, kT, energies)
     data_new = dropdims(sum(res.data, dims=1), dims=1) * sc.Δω
-    InstantIntensities(res.crystal, res.qpts, data_new)
+    StaticIntensities(res.crystal, res.qpts, data_new)
 end
 
-function intensities_instant(sc::SampledCorrelationsStatic, qpts; formfactors=nothing)
-    # Contrary to the name, this will actually return an InstantIntensities
+function intensities_static(sc::SampledCorrelationsStatic, qpts; formfactors=nothing)
+    # Contrary to the name, this will actually return an StaticIntensities
     intensities(sc.parent, qpts; formfactors, kT=nothing, energies=:available)
 end
 

src/SampledCorrelations/SampledCorrelations.jl

@@ -124,7 +124,7 @@ equilibrium. A classical spin dynamics trajectory will be simulated of
 sufficient length to achieve the target energy resolution. The resulting data
 can can then be extracted as pair-correlation [`intensities`](@ref) with
 appropriate classical-to-quantum correction factors. See also
-[`intensities_instant`](@ref), which integrates over the available energy range.
+[`intensities_static`](@ref), which integrates over energy.
 """
 function SampledCorrelations(sys::System; measure, energies, dt, calculate_errors=false)
     if isnothing(energies)
@@ -182,13 +182,13 @@ end
 """
     SampledCorrelationsStatic(sys::System; measure)
 
-An object to accumulate samples of static pair correlations. Similar to
+An object to accumulate samples of static pair correlations. It is similar to
 [`SampledCorrelations`](@ref), but no time-integration will be performed on
-calls to [`add_sample!`](@ref). As a result, dynamical [`intensities`](@ref)
-data will be unavailable for `SampledCorrelationsStatic`. Furthermore,
-[`intensities_instant`](@ref) data is associated with the classical Boltzmann
-distribution, and misses classical-to-quantum corrections that can be captured
-by `SampledCorrelations`.
+calls to [`add_sample!`](@ref). The resulting object can be used with
+[`intensities_static`](@ref) to calculate statistics from the classical
+Boltzmann distribution. Dynamical [`intensities`](@ref) data, however, will be
+unavailable. Similarly, classical-to-quantum corrections that rely on the
+excitation spectrum cannot be performed.
 """
 struct SampledCorrelationsStatic
     parent :: SampledCorrelations

src/SpinWaveTheory/DispersionAndIntensities.jl

@@ -141,12 +141,13 @@ end
     intensities_bands(swt::SpinWaveTheory, qpts; formfactors=nothing, kT=0)
 
 Calculate spin wave excitation bands for a set of q-points in reciprocal space.
-This calculation is analogous [`intensities`](@ref), but does not perform
-broadening of the bands.
+This calculation is analogous to [`intensities`](@ref), but does not perform
+line broadening of the bands.
 """
-function intensities_bands(swt::SpinWaveTheory, qpts; formfactors=nothing, kT=0)
+function intensities_bands(swt::SpinWaveTheory, qpts; formfactors=nothing, kT=0, with_negative=false)
     (; sys, measure) = swt
     isempty(measure.observables) && error("No observables! Construct SpinWaveTheory with a `measure` argument.")
+    with_negative && error("Option `with_negative=true` not yet supported.")
 
     qpts = convert(AbstractQPoints, qpts)
     cryst = orig_crystal(sys)
@@ -249,7 +250,7 @@ end
     intensities(swt::SpinWaveTheory, qpts; energies, kernel, formfactors=nothing, kT=0)
     intensities(swt::SampledCorrelations, qpts; energies, kernel=nothing, formfactors=nothing, kT)
 
-Calculate pair correlation intensities for a set of q-points in reciprocal
+Calculates pair correlation intensities for a set of ``𝐪``-points in reciprocal
 space.
 
 Traditional spin wave theory calculations are performed with an instance of
@@ -277,24 +278,33 @@ function intensities(swt::AbstractSpinWaveTheory, qpts; energies, kernel::Abstra
 end
 
 """
-    intensities_instant(sc::SpinWaveTheory, qpts; formfactors=nothing, kT=0)
-    intensities_instant(sc::SampledCorrelations, qpts; formfactors=nothing, kT)
-    intensities_instant(sc::SampledCorrelationsStatic, qpts; formfactors=nothing)
-
-Calculate the instantaneous (equal-time) correlations for a set of
-``𝐪``-points. This is the integral of ``\\mathcal{S}(𝐪, ω)`` over all energies
-``ω``.
-
-In [`SpinWaveTheory`](@ref) the integral can be realized as a discrete sum over
-bands. In [`SampledCorrelations`](@ref) there is an analogous integral over the
-available energies.
+    intensities_static(sc::SpinWaveTheory, qpts; bounds=(-Inf, Inf), formfactors=nothing, kT=0)
+    intensities_static(sc::SampledCorrelations, qpts; bounds=(-Inf, Inf), formfactors=nothing, kT)
+    intensities_static(sc::SampledCorrelationsStatic, qpts; formfactors=nothing)
+
+Like [`intensities`](@ref), but integrates the dynamical correlations
+``\\mathcal{S}(𝐪, ω)`` over a range of energies ``ω``. By default, the
+integration `bounds` are ``(-∞, ∞)``, yielding the instantaneous (equal-time)
+correlations.
+
+In [`SpinWaveTheory`](@ref) the integral will be realized as a sum over discrete
+bands. A [`SampledCorrelations`](@ref) object will have a finite grid of
+available `energies`, which will constrain the domain of integration. A
+[`SampledCorrelationsStatic`](@ref) object stores no dynamical data; here, the
+return value represents instantaneous correlations for the _classical_ Boltzmann
+distribution.
 
 The parameter `kT` can be used to account for the quantum thermal occupation of
 excitations at finite temperature. For details, see the documentation in
 [`intensities`](@ref).
 """
-function intensities_instant(swt::AbstractSpinWaveTheory, qpts; formfactors=nothing, kT=0)
-    res = intensities_bands(swt, qpts; formfactors, kT)
-    data_new = dropdims(sum(res.data, dims=1), dims=1)
-    InstantIntensities(res.crystal, res.qpts, data_new)
+function intensities_static(swt::AbstractSpinWaveTheory, qpts; bounds=(-Inf, Inf), formfactors=nothing, kT=0)
+    res = intensities_bands(swt, qpts; formfactors, kT)  # TODO: with_negative=true
+    data_reduced = zeros(eltype(res.data), size(res.data)[2:end])
+    for ib in axes(res.data, 1), iq in CartesianIndices(data_reduced)
+        if bounds[1] <= res.disp[ib, iq] < bounds[2]
+            data_reduced[iq] += res.data[ib, iq]
+        end
+    end
+    StaticIntensities(res.crystal, res.qpts, data_reduced)
 end

src/Sunny.jl

@@ -91,7 +91,7 @@ include("SpinWaveTheory/HamiltonianSUN.jl")
 include("SpinWaveTheory/DispersionAndIntensities.jl")
 include("SpinWaveTheory/LSWTCorrections.jl")
 export SpinWaveTheory, excitations, excitations!, dispersion, intensities, intensities_bands,
-    intensities_instant # TODO
+    intensities_static
 
 include("Spiral/LuttingerTisza.jl")
 include("Spiral/SpiralEnergy.jl")

test/test_correlation_sampling.jl

@@ -79,12 +79,12 @@
 
 
     # Test static from dynamic intensities working
-    is_static = intensities_instant(sc, qgrid; kT=nothing)
+    is_static = intensities_static(sc, qgrid; kT=nothing)
     total_intensity_static = sum(is_static.data)
     @test isapprox(total_intensity_static, total_intensity_trace * sc.Δω; atol=1e-9)  # Order of summation can lead to very small discrepancies
 
     # Test quantum-to-classical increases intensity
-    is_static_c2q = intensities_instant(sc, qgrid; kT=0.1)
+    is_static_c2q = intensities_static(sc, qgrid; kT=0.1)
     total_intensity_static_c2q = sum(is_static_c2q.data)
     @test total_intensity_static_c2q > total_intensity_static 
 
@@ -93,7 +93,7 @@
     thermalize_simple_model!(sys; kT=0.1)
     ic = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=false))
     add_sample!(ic, sys)
-    true_static_vals = intensities_instant(ic, qgrid)
+    true_static_vals = intensities_static(ic, qgrid)
     true_static_total = sum(true_static_vals.data)
     @test isapprox(true_static_total / prod(sys.latsize), 1.0; atol=1e-12)
 end

test/test_intensities_setup.jl

@@ -16,7 +16,7 @@
         # Dipole magnetization observables (classical)
         sc = SampledCorrelationsStatic(sys; measure=ssf_custom((q, ssf) -> ssf, sys; apply_g=true))
         add_sample!(sc, sys)
-        is = intensities_instant(sc, Sunny.QPoints([[0,0,0]]))
+        is = intensities_static(sc, Sunny.QPoints([[0,0,0]]))
         corr_mat = is.data[1]
 
         # Compute magnetization correlation "by hand", averaging over sites