Fix wraparound bug (#246)

Avoids artifacts in structure factor intensities estimated from `dynamical_correlations`.

The option `process_trajectory=:symmetrize` is no longer needed, and now deprecated.

Clarify and improve structure factor documentation.

---------

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

docs/src/structure-factor.md

@@ -334,13 +334,12 @@ and describes dynamical correlations of magnetic moments. It will differ
 nontrivially from the spin-spin correlations if the $g_j$-tensor varies with
 sublattice $j$.
 
-## Calculating the structure factor in Sunny
+## Conventions for the Sunny-calculated structure factor
 
 Calculating the structure factor involves several steps, with various possible
-settings. Sunny provides a number of tools to facilitate this calculation and to
-extract information from the results. These tools are briefly outlined below.
-Please see the documentation tutorials for example usage. Detailed function
-information is available in the [Library API](@ref).
+settings. Sunny provides tools to facilitate this calculation and to extract
+information from the results. For details, please see our [tutorials](@ref "2.
+Spin wave simulations of CoRh₂O₄") as well as the complete [Library API](@ref).
 
 Sunny will calculate the structure factor in dimensionless, intensive units,
 
@@ -364,38 +363,62 @@ through normalization by $N_\mathrm{cells}$. Its numerical value _is_ dependent
 on the size of the chemical cell (larger chemical cells lead to greater
 intensities reported by Sunny).
 
+In most cases, users will calculate the structure factor within linear
+[`SpinWaveTheory`](@ref), whereby magnetic excitations are approximated as
+Holstein-Primakoff bosons. This calculation technique is relatively
+straightforward and efficient. Linear spin wave theory has two primary
+limitations, however. It cannot account for thermal fluctuations beyond the
+harmonic approximation, and it scales poorly in the size of the magnetic cell
+size (e.g., as needed to study systems with chemical disorder). The efficiency
+limitation can be overcome with [recently proposed
+algorithms](https://arxiv.org/abs/2312.08349) that are planned for
+[implementation in Sunny](https://github.com/SunnySuite/Sunny.jl/pull/92).
+However, the study of finite temperature fluctuations requires a calculation
+method that is entirely different from linear spin wave theory.
+
 ## Estimating stucture factors with classical dynamics
 
-Classical dynamics may be used to estimate structure factor data by analyzing
-the spin-spin correlations of dynamical trajectories. This is fundamentally a
-Monte Carlo approach, as the trajectories must be started from an initial spin
-configuration that is sampled at thermal equilibrium. (Note that it is not
-possible to estimate a true T=0 dynamical structure factor using this method,
-but the temperature may be very low.) Samples are accumulated into a
+Finite temperature structure factor intensities can be estimated from the
+dynamical correlations of classical spin dynamics (e.g. Landau-Lifshitz, or its
+SU($N$) generalization). This is fundamentally a Monte Carlo approach, as the
+trajectories must be initialized to a spin configuration that is sampled from
+the finite-temperature thermal equilibrium. Samples are accumulated into a
 `SampledCorrelations`, from which intensity information may be extracted. The
 user does not typically build their own `SampledCorrelations` but instead
 initializes one by calling either `dynamical_correlations` or
 `instant_correlations`, as described below.
 
 ### Estimating a dynamical structure factor: ``𝒮(𝐤,ω)``
 
-A `SampledCorrelations` for estimating the dynamical structure factor,
-$𝒮^{αβ}(𝐤,ω)$, may be created by calling [`dynamical_correlations`](@ref). This
-requires three keyword arguments. These will determine the dynamics used to
-calculate samples and, consequently, the $ω$ information that will be available. 
-
-1. `dt`: Determines the step size used for simulating the dynamics. Typically
-   this will be limited by numerical stability. The function
-   [`suggest_timestep`](@ref) can recommend a value.
-2. `ωmax`: Sets the maximum resolved energy. Very large `ωmax` may require
-   smaller `dt`. 
-3. `nω`: Determines the number of energy bins to resolve. A larger number will
-   require more calculation time.
+A `SampledCorrelations` object for estimating the dynamical structure factor is
+created by calling [`dynamical_correlations`](@ref). This requires three keyword
+arguments. These will determine the dynamics used to calculate samples and,
+consequently, the $ω$ information that will be available. 
+
+1. `ωmax`: Sets the maximum resolved energy.
+2. `nω`: Sets the number of discrete energy values to resolve. The corresponding
+   energy resolution is approximately `Δω ≈ ωmax / nω`. To estimate the
+   structure factor with resolution `Δω`, Sunny must integrate a classical spin
+   dynamic trajectory over a time-scale of order `1 / Δω`. Computational cost
+   therefore scales approximately linearly in `nω`.
+3. `dt`: Determines the step size for dynamical time-integration. Larger is more
+   efficient, but the choice will be limited by the stability and accuracy
+   requirements of the [`ImplicitMidpoint`](@ref) integration method. The
+   function [`suggest_timestep`](@ref) can recommend a good value. The inverse
+   of `ωmax` also imposes an upper bound on `dt`.
+
+!!! warning "Intensity scale"
+
+    As of Sunny v0.5, the ``𝒮(𝐤,ω)`` intensities calculated with
+    `dynamical_correlations` will be scaled by a prefactor of `Δω ≈ ωmax / nω`,
+    the discretization in energy space. This prefactor will be removed in a
+    future Sunny version. See
+    [Issue 264](https://github.com/SunnySuite/Sunny.jl/issues/264).
 
 A sample may be added by calling `add_sample!(sc, sys)`. The input `sys` must be
 a spin configuration in good thermal equilibrium, e.g., using the continuous
 [`Langevin`](@ref) dynamics or using single spin flip trials with
-[`LocalSampler`](@ref). The statistical quality of the $𝒮^{αβ}(𝐤,ω)$ can be
+[`LocalSampler`](@ref). The statistical quality of the $𝒮(𝐤,ω)$ can be
 improved by repeatedly generating decorrelated spin configurations in `sys` and
 calling `add_sample!` on each configuration.
 
@@ -413,22 +436,21 @@ end
 The calculation may be configured in a number of ways; see the
 [`dynamical_correlations`](@ref) documentation for a list of all keywords.
 
-
 ### Estimating an instantaneous ("static") structure factor: ``𝒮(𝐤)``
 
-Sunny provides two methods for calculating instantaneous, or static, structure
-factors: $𝒮^{αβ}(𝐤)$. The first involves calculating spatial spin-spin
-correlations at single time slices. The second involves calculating a dynamic
-structure factor first and integrating out the $ω$ information. The advantage of
-the latter approach is that it enables application of an $ω$-dependent
-classical-to-quantum rescaling of structure factor intensities, a method that
-should be preferred whenever comparing results to experimental data or spin wave
-calculations. A disadvantage of this approach is that it is computationally more
-expensive. There are also many cases when it is not straightforward to calculate
-a meaningful dynamics, as when working with Ising spins. In this section we will
-discuss how to calculate instantaneous structure factors from static spin
-configurations. Information about calculating instantaneous data from a
-dynamical correlations can be found in the following section.
+Sunny provides two methods for calculating instantaneous structure factors
+$𝒮(𝐤)$. The first involves calculating spatial spin-spin correlations at
+single time slices. The second involves calculating a dynamic structure factor
+first and then integrating over $ω$. The advantage of the latter approach is
+that it enables application of an $ω$-dependent classical-to-quantum rescaling
+of structure factor intensities, a method that should be preferred whenever
+comparing results to experimental data or spin wave calculations. A disadvantage
+of this approach is that it is computationally more expensive. There are also
+many cases when it is not straightforward to calculate a meaningful dynamics, as
+when working with Ising spins. In this section we will discuss how to calculate
+instantaneous structure factors from static spin configurations. Information
+about calculating instantaneous data from a dynamical correlations can be found
+in the following section.
 
 The basic usage for the instantaneous case is very similar to the dynamic case,
 except one calls [`instant_correlations`](@ref) instead of
@@ -448,15 +470,18 @@ out the energy axis. An approach to doing this is described in the next section.
 The basic function for extracting information from a `SampledCorrelations` at a
 particular wave vector, $𝐤$, is [`intensities_interpolated`](@ref). It takes a
 `SampledCorrelations`, a _list_ of wave vectors, and an
-[`intensity_formula`](@ref). The `intensity_formula` specifies how to contract and correct
-correlation data to arrive at a physical intensity.
-A simple example is `formula = intensity_formula(sc, :perp)`, which will
-instruct Sunny apply polarization corrections: $\sum_{αβ}(I-q_α q_β) 𝒮^{αβ}(𝐤,ω)$.
-An intensity at the wave vector $𝐤 = (𝐛_2 + 𝐛_3)/2$
-may then be retrieved with  `intensities_interpolated(sf, [[0.0, 0.5, 0.5]], formula)` . 
+[`intensity_formula`](@ref). The `intensity_formula` specifies how to contract
+and correct correlation data to arrive at a physical intensity. A simple example
+is `formula = intensity_formula(sc, :perp)`, which will instruct Sunny apply
+polarization corrections: $\sum_{αβ}(I-q_α q_β) 𝒮^{αβ}(𝐤,ω)$. An intensity at
+the wave vector $𝐤 = (𝐛_2 + 𝐛_3)/2$ may then be retrieved with
+`intensities_interpolated(sf, [[0.0, 0.5, 0.5]], formula)` .
 `intensities_interpolated` returns a list of `nω` elements at each wavevector.
 The corresponding $ω$ values can be retrieved by calling
-[`available_energies`](@ref) on `sf`.
+[`available_energies`](@ref) on `sf`. Note that there will always be some amount
+of "blurring" between neighboring energy values. This blurring originates from
+the finite-length dynamical trajectories following the algorithm [specified
+here](https://github.com/SunnySuite/Sunny.jl/pull/246#issuecomment-2119294846).
 
 Since Sunny only calculates the structure factor on a finite lattice when
 performing classical simulations, it is important to realize that exact

docs/src/versions.md

@@ -3,15 +3,19 @@
 ## v0.5.10
 (In development)
 
-* In dipole mode, fix bugs to support the case that spin-``S`` varies between
-  sites. In SU(``N``) mode, however, there is still no support for varying
-  the Hilbert space dimension ``N`` between sites.
 * [`view_crystal`](@ref) called on a [`System`](@ref) now optionally shows
   spin or magnetic dipoles.
 * Interactions for [`enable_dipole_dipole!`](@ref) are now supported in linear
   spin wave theory, with proper Ewald summation. For a faster alternative, the
   experimental function [`modify_exchange_with_truncated_dipole_dipole!`](@ref)
   will accept a real-space cutoff.
+* Intensities calculated with [`dynamical_correlations`](@ref) now avoid
+  "smearing artifacts" at low-energy (long-timescale) modes. See [PR
+  246](https://github.com/SunnySuite/Sunny.jl/pull/246) for details. This
+  eliminates the need for `process_trajectory=:symmetrize`.
+* In dipole mode, having spin-``S`` vary between sites was previously broken,
+  and is now fixed. In SU(``N``) mode, however, there is still no support for
+  varying the Hilbert space dimension ``N`` between sites.
 * Long-range dipole-dipole was previously broken for systems with multiple
   cells, but is now fixed.
 * General biquadratic interactions (beyond scalar) in dipole mode are fixed in

docs/src/writevtk.md

@@ -44,7 +44,7 @@ end
 
 ωmax=10.
 
-dsf = dynamical_correlations(sys; dt, nω=48, ωmax, process_trajectory=:symmetrize)
+dsf = dynamical_correlations(sys; dt, nω=48, ωmax)
 
 nsamples = 10
 for _ in 1:nsamples

examples/03_LLD_CoRh2O4.jl

@@ -132,7 +132,7 @@ heatmap(q1s, q2s, iq;
 dt = 2*langevin.dt
 ωmax = 6.0  # Maximum energy to resolve (meV)
 nω = 50     # Number of energies to resolve
-sc = dynamical_correlations(sys; dt, nω, ωmax, process_trajectory=:symmetrize)
+sc = dynamical_correlations(sys; dt, nω, ωmax)
 
 # Use Langevin dynamics to sample spin configurations from thermal equilibrium.
 # For each sample, use [`add_sample!`](@ref) to run a classical spin dynamics

src/Intensities/Binning.jl

@@ -208,7 +208,7 @@ function unit_resolution_binning_parameters(ωvals,latsize,args...)
     params
 end
 
-unit_resolution_binning_parameters(sc::SampledCorrelations; kwargs...) = unit_resolution_binning_parameters(available_energies_including_zero(sc),sc.latsize,sc;kwargs...)
+unit_resolution_binning_parameters(sc::SampledCorrelations; negative_energies=false,kwargs...) = unit_resolution_binning_parameters(available_energies_including_zero(sc;negative_energies),sc.latsize,sc;kwargs...)
 
 function unit_resolution_binning_parameters(ωvals::AbstractVector{Float64})
     if !all(abs.(diff(diff(ωvals))) .< 1e-12)
@@ -382,7 +382,7 @@ function intensities_binned(sc::SampledCorrelations, params::BinningParameters,
     return_type = typeof(formula).parameters[1]
     output_intensities = zeros(return_type,numbins...)
     output_counts = zeros(Float64,numbins...)
-    ωvals = available_energies_including_zero(sc)
+    ωvals = available_energies_including_zero(sc;negative_energies=true)
 
     # Find an axis-aligned bounding box containing the histogram.
     # The AABB needs to be in r.l.u for the (possibly reshaped) crystal
@@ -475,11 +475,28 @@ function intensities_binned(sc::SampledCorrelations, params::BinningParameters,
             end
         end
     end
-    return output_intensities, output_counts
+
+    # `output_intensities` is in units of S²/BZ/fs, and includes the sum
+    # of `output_counts`-many individual scattering intensities.
+    # To give the value integrated over the bin, we need to multiply by
+    # the binwidth, Δω×ΠᵢΔqᵢ. But Δq/BZ = 1/N, where N is the number of
+    # bins covering one BZ, which is itself equal to the latsize.
+    #
+    # To find the number of bins covering one BZ, we first compute the volume of the BZ
+    # in histogram label space: it is det(covectors[1:3,1:3]). Next, we compute the volume
+    # of one bin in label space: it is prod(binwidth[1:3]).
+    # Then, N = det(covectors[1:3,1:3])/prod(binwidth[1:3]).
+    #
+    # For the time axis, Δω/fs = 1/N, where N is the number of frequency bins
+    # 
+    # So the division by N here makes it so the result has units of
+    # raw S² (to be summed over M-many BZs to recover M-times the sum rule)
+    N_bins_in_BZ = abs(det(covectors[1:3,1:3])) / prod(binwidth[1:3])
+    return output_intensities ./ N_bins_in_BZ ./ length(ωvals), output_counts
 end
 
-function available_energies_including_zero(x)
-    ωs = available_energies(x)
+function available_energies_including_zero(x;kwargs...)
+    ωs = available_energies(x;kwargs...)
     # Special case due to NaN definition of instant_correlations
     (length(ωs) == 1 && isnan(ωs[1])) ? [0.] : ωs
 end