Update comments

examples/03_LLD_CoRh2O4.jl

@@ -205,5 +205,5 @@ heatmap(radii, ωs, output;
         ylabel="Energy Transfer (meV)", 
         aspect = 1.4,
     ),
-    colorrange = (0, 20.0) 
+    colorrange = (0, 20.0)
 )

src/Intensities/Interpolation.jl

@@ -150,15 +150,6 @@ function intensities_interpolated(sc::SampledCorrelations, qs, formula::Classica
     # Call type-stable version of the function
     intensities_interpolated!(intensities, sc, qs, ωvals, interp, formula, stencil_info, Val(return_type))
 
-    # `intensities` is in units of S²/BZ/fs, but we want to convert it to S²/R.L.U./meV.
-    # The conversion factor from BZ→R.L.U. is 1 because one BZ of the Bravais lattice is [0,1]³.
-    if !isnan(sc.Δω)
-      # If there is a frequency axis, we need to convert 1/fs→1/meV.
-      # The conversion factor is:
-      meV_per_fs = 1/(size(sc.samplebuf,6) * sc.Δω)
-      intensities .*= meV_per_fs
-    end
-
     return intensities
 end
 
@@ -186,6 +177,19 @@ function intensities_interpolated!(intensities, sc::SampledCorrelations, q_targe
             end
         end
     end
+
+    # If Δω is nan then this is an "instantaneous" structure factor, and no
+    # processing of the ω axis is needed.
+    if !isnan(sc.Δω)
+        n_all_ω = size(sc.samplebuf,6)
+        # The usual discrete fourier transform (implemented by FFT) produces an
+        # extensive result. Dividing by n_all_ω makes the result intensive
+        # (i.e., intensity scale becomes independent of the number of ω values).
+        # Additionally dividing by Δω produces a density in ω space. Note that
+        # intensities is already a density in q space.
+        intensities ./= (n_all_ω * sc.Δω)
+    end
+
     return intensities
 end
 

src/SampledCorrelations/CorrelationSampling.jl

@@ -60,14 +60,6 @@ function new_sample!(sc::SampledCorrelations, sys::System)
     return nothing
 end
 
-
-function subtract_mean!(sc::SampledCorrelations)
-    (; samplebuf) = sc
-    nsteps = size(samplebuf, 6)
-    meanvals = sum(samplebuf, dims=6) ./ nsteps
-    samplebuf .-= meanvals
-end
-
 function no_processing(::SampledCorrelations)
     nothing
 end

src/SampledCorrelations/SampledCorrelations.jl

@@ -7,20 +7,20 @@ initialized by calling either [`dynamical_correlations`](@ref) or
 """
 struct SampledCorrelations{N}
     # 𝒮^{αβ}(q,ω) data and metadata
-    data           :: Array{ComplexF64, 7}                 # Raw SF data for 1st BZ (numcorrelations × natoms × natoms × latsize × energy)
+    data           :: Array{ComplexF64, 7}                 # Raw SF with sublattice indices (ncorrs × natoms × natoms × latsize × nω)
     M              :: Union{Nothing, Array{Float64, 7}}    # Running estimate of (nsamples - 1)*σ² (where σ² is the variance of intensities)
     crystal        :: Crystal                              # Crystal for interpretation of q indices in `data`
     origin_crystal :: Union{Nothing,Crystal}               # Original user-specified crystal (if different from above) -- needed for FormFactor accounting
     Δω             :: Float64                              # Energy step size (could make this a virtual property)  
     observables    :: ObservableInfo
 
     # Specs for sample generation and accumulation
-    samplebuf    :: Array{ComplexF64, 6}   # New sample buffer
-    corrbuf      :: Array{ComplexF64, 4}   # Buffer for real-time correlations 
-    space_fft!   :: FFTW.AbstractFFTs.Plan # Pre-planned FFT
-    time_fft!    :: FFTW.AbstractFFTs.Plan # Pre-planned FFT
-    corr_fft!    :: FFTW.AbstractFFTs.Plan # Pre-planned FFT
-    corr_ifft!   :: FFTW.AbstractFFTs.Plan # Pre-planned FFT
+    samplebuf    :: Array{ComplexF64, 6}   # Buffer for observables (nobservables × latsize × natoms × nsnapshots)
+    corrbuf      :: Array{ComplexF64, 4}   # Buffer for correlations (latsize × nω)
+    space_fft!   :: FFTW.AbstractFFTs.Plan # Pre-planned FFT for samplebuf
+    time_fft!    :: FFTW.AbstractFFTs.Plan # Pre-planned FFT for samplebuf
+    corr_fft!    :: FFTW.AbstractFFTs.Plan # Pre-planned FFT for corrbuf
+    corr_ifft!   :: FFTW.AbstractFFTs.Plan
     measperiod   :: Int                    # Steps to skip between saving observables (downsampling for dynamical calcs)
     apply_g      :: Bool                   # Whether to apply the g-factor
     dt           :: Float64                # Step size for trajectory integration 
@@ -58,9 +58,8 @@ Base.getproperty(sc::SampledCorrelations, sym::Symbol) = sym == :latsize ? size(
 
 function clone_correlations(sc::SampledCorrelations{N}) where N
     dims = size(sc.data)[2:4]
-    n_all_ω = size(sc.data, 7)
     # Avoid copies/deep copies of C-generated data structures
-    space_fft! = 1/√(prod(dims)) * FFTW.plan_fft!(sc.samplebuf, (2,3,4))
+    space_fft! = 1/√prod(dims) * FFTW.plan_fft!(sc.samplebuf, (2,3,4))
     time_fft! = FFTW.plan_fft!(sc.samplebuf, 6)
     corr_fft! = FFTW.plan_fft!(sc.corrbuf, 4)
     corr_ifft! = FFTW.plan_ifft!(sc.corrbuf, 4)
@@ -183,7 +182,7 @@ function dynamical_correlations(sys::System{N}; dt=nothing, Δt=nothing, nω, ω
     # the same way because the number of estimates varies with Δt. These conventions
     # ensure consistency with this spec:
     # https://sunnysuite.github.io/Sunny.jl/dev/structure-factor.html
-    space_fft! = 1/√(prod(sys.latsize)) * FFTW.plan_fft!(samplebuf, (2,3,4))
+    space_fft! = 1/√prod(sys.latsize) * FFTW.plan_fft!(samplebuf, (2,3,4))
     time_fft! = FFTW.plan_fft!(samplebuf, 6)
     corr_fft! = FFTW.plan_fft!(corrbuf, 4)
     corr_ifft! = FFTW.plan_ifft!(corrbuf, 4)

test/test_correlation_sampling.jl

@@ -37,7 +37,7 @@
     Δω = sc.Δω
     add_sample!(sc, sys)
     qgrid = available_wave_vectors(sc)
-    Δq³ = 1/prod(sys.latsize) # Spacing between available wave vectors
+    Δq³ = 1/prod(sys.latsize) # Fraction of a BZ for one cell
     formula = intensity_formula(sc,:trace)
     Sqw = intensities_interpolated(sc, qgrid, formula; negative_energies=true)
     expected_sum_rule = Sunny.norm2(sys.dipoles[1]) # S^2 classical sum rule
@@ -161,6 +161,6 @@ end
     # println(round.(data; digits=10))
 
     # Compare with reference
-    data_golden = (1/sc.Δω) .* [0.9264336057 1.6408721819 -0.0 -0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0; 0.0252093265 0.0430096149 0.6842894708 2.9288686 4.4296891436 5.546520481 5.2442499151 1.4849591286 -0.1196709189 0.0460970304]
+    data_golden = [1.5688306301 2.7786670553 0.0 -0.0 -0.0 -0.0 0.0 0.0 0.0 0.0; 0.0426896901 0.0728328515 1.158781671 4.9597712594 7.5012736668 9.3925254521 8.880657878 2.5146425508 -0.2026517626 0.0780611074]
     @test data ≈ data_golden
 end

test/test_intensities_setup.jl

@@ -96,10 +96,15 @@ end
     add_sample!(sc, sys)
 
     sum_rule_ixs = Sunny.Trace(sc.observables).indices
-    sub_lat_sum_rules = sum(sc.data[sum_rule_ixs,:,:,:,:,:,:],dims = [1,4,5,6,7])[1,:,:,1,1,1,1]
-    # Since sc.data is in units of S²/BZ/fs, each entry is already an estimate of S².
-    # To get the frequency-averaged estimate, we need to divide by the number of estimators:
-    sub_lat_sum_rules ./= prod(size(sc.data)[[4,5,6,7]])
+    sub_lat_sum_rules = sum(sc.data[sum_rule_ixs,:,:,:,:,:,:], dims=[1,4,5,6,7])[1,:,:,1,1,1,1]
+
+    Δq³ = 1/prod(sys.latsize) # Fraction of a BZ for one cell
+    n_all_ω = size(sc.data, 7)
+    # Intensities in sc.data are a density in q, but already integrated over dω
+    # bins, and then scaled by n_all_ω. Therefore, we need the factor below to
+    # convert the previous sum to an integral. (See same logic in
+    # intensities_interpolated function.)
+    sub_lat_sum_rules .*= Δq³ / n_all_ω
 
     # SU(N) sum rule for S = 1/2:
     # ⟨∑ᵢSᵢ²⟩ = 3/4 on every site, but because we're classical, we
@@ -111,7 +116,7 @@ end
     expected_sum = gS_squared
     # This sum rule should hold for each sublattice, independently, and only
     # need to be taken over a single BZ (which is what sc.data contains) to hold:
-    @test [sub_lat_sum_rules[i,i] for i = 1:Sunny.natoms(sc.crystal)] ≈ expected_sum * ones(ComplexF64,Sunny.natoms(sc.crystal))
+    @test [sub_lat_sum_rules[i,i] for i in 1:Sunny.natoms(sc.crystal)] ≈ expected_sum * ones(ComplexF64,Sunny.natoms(sc.crystal))
 
     formula = intensity_formula(sc,:trace)
     # The polyatomic sum rule demands going out 4 BZ's for the diamond crystal