Classical intensities as density

docs/src/versions.md

@@ -3,14 +3,13 @@
 ## v0.6.0
 (In progress)
 
-* Fix direction of magnetic moment and direction of time ([PR
-  247](https://github.com/SunnySuite/Sunny.jl/pull/247)).
-* Fix various sign conventions: direction of magnetic moment, ([Issue
-  190](https://github.com/SunnySuite/Sunny.jl/issues/190)), direction of time in
-  Landau-Lifshitz spin dynamics (resolved in [PR
-  247](https://github.com/SunnySuite/Sunny.jl/pull/247)), and direction of
-  momentum transfer when calculating the structure factor from classical
-  dynamics ([Issue 270](https://github.com/SunnySuite/Sunny.jl/issues/270)).
+* Various correctness fixes. Direction of magnetic moment and direction of time
+  in LL dynamics ([PR 247](https://github.com/SunnySuite/Sunny.jl/pull/247)).
+  Direction of momentum transfer ([PR
+  271](https://github.com/SunnySuite/Sunny.jl/pull/271)). Dynamical structure
+  factor intensity now a proper energy density ([PR
+  272](https://github.com/SunnySuite/Sunny.jl/pull/272)); consequently, color
+  ranges in plots may need to be rescaled.
 
 ## v0.5.11
 (June 2, 2024)

examples/04_GSD_FeI2.jl

@@ -223,7 +223,7 @@ fig
 # for `Fe2` magnetic ions, and a dipole polarization correction `:perp`.
 
 formfactors = [FormFactor("Fe2"; g_lande=3/2)]
-new_formula = intensity_formula(sc, :perp; kT, formfactors = formfactors)
+new_formula = intensity_formula(sc, :perp; kT, formfactors)
 
 # Frequently, one wants to extract energy intensities along lines that connect
 # special wave vectors--a so-called "spaghetti plot". The function
@@ -283,11 +283,11 @@ ax_top = Axis(fig[1,1],ylabel = "meV",xticklabelrotation=π/8,xticklabelsize=12;
 ax_bottom = Axis(fig[2,1],ylabel = "meV",xticks = (markers, string.(points)),xticklabelrotation=π/8,xticklabelsize=12)
 
 heatmap!(ax_top,1:size(is_interpolated,1), ωs, is_interpolated;
-    colorrange=(0.0,0.07),
+    colorrange=(0.0, 1.0),
 )
 
 heatmap!(ax_bottom,1:size(is_binned,1), ωs, is_binned;
-    colorrange=(0.0,0.05),
+    colorrange=(0.0, 1.5e-6),
 )
 
 fig

examples/08_Scattering_conventions.jl → examples/08_Momentum_Conventions.jl

@@ -89,6 +89,6 @@ disp_swt, intens_swt = intensities_bands(swt, path, formula);
 fig = Figure()
 ax = Axis(fig[1,1]; aspect=1.4, ylabel="ω (meV)", xlabel="𝐪 (r.l.u.)",
           xticks, xticklabelrotation=π/10)
-heatmap!(ax, 1:size(data, 1), available_energies(sc), data;  colorrange=(0.0, 0.8))
+heatmap!(ax, 1:size(data, 1), available_energies(sc), data;  colorrange=(0.0, 50.0))
 lines!(ax, disp_swt[:,1]; color=:magenta, linewidth=2)
 fig

src/Intensities/Interpolation.jl

@@ -150,6 +150,9 @@ 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))
 
+    # Make density
+    intensities ./= isnan(sc.Δω) ? 1.0 : sc.Δω
+
     return intensities
 end
 

src/SampledCorrelations/CorrelationSampling.jl

@@ -73,7 +73,7 @@ function no_processing(::SampledCorrelations)
 end
 
 function accum_sample!(sc::SampledCorrelations; window)
-    (; data, M, observables, samplebuf, nsamples, space_fft!, time_fft!) = sc
+    (; data, M, observables, samplebuf, corrbuf, nsamples, space_fft!, time_fft!, corr_fft!, corr_ifft!) = sc
     natoms = size(samplebuf)[5]
 
     num_time_offsets = size(samplebuf,6)
@@ -122,22 +122,24 @@ function accum_sample!(sc::SampledCorrelations; window)
         # According to Sunny convention, the correlation is between
         # α† and β. This conjugation implements both the dagger on the α
         # as well as the appropriate spacetime offsets of the correlation.
-        corr = FFTW.ifft(conj.(sample_α) .* sample_β,4) ./ n_contrib
+        @. corrbuf = conj(sample_α) * sample_β
+        corr_ifft! * corrbuf
+        corrbuf ./= n_contrib
 
         if window == :cosine
           # Apply a cosine windowing to force the correlation at Δt=±(T-1) to be zero
           # to force periodicity. In terms of the spectrum S(ω), this applys a smoothing
           # with a characteristic lengthscale of O(1) frequency bins.
           window_func = cos.(range(0,π,length = num_time_offsets + 1)[1:end-1]).^2
-          corr .*= reshape(window_func,1,1,1,num_time_offsets)
+          corrbuf .*= reshape(window_func,1,1,1,num_time_offsets)
         end
 
-        FFTW.fft!(corr,4)
+        corr_fft! * corrbuf
 
         if isnothing(M)
             for k in eachindex(databuf)
                 # Store the diff for one complex number on the stack.
-                diff = corr[k] - databuf[k]
+                diff = corrbuf[k] - databuf[k]
 
                 # Accumulate into running average
                 databuf[k] += diff * (1/count)
@@ -149,15 +151,14 @@ function accum_sample!(sc::SampledCorrelations; window)
                 μ_old = databuf[k]
 
                 # Update running mean.
-                matrixelem = sample_α[k] * conj(sample_β[k])
-                databuf[k] += (matrixelem - databuf[k]) * (1/count)
+                databuf[k] += (corrbuf[k] - databuf[k]) / count
                 μ = databuf[k]
 
                 # Update variance estimate.
                 # Note that the first term of `diff` is real by construction
                 # (despite appearances), but `real` is explicitly called to
                 # avoid automatic typecasting errors caused by roundoff.
-                Mbuf[k] += real((matrixelem - μ_old)*conj(matrixelem - μ))
+                Mbuf[k] += real((corrbuf[k] - μ_old)*conj(corrbuf[k] - μ))
             end
         end
     end

src/SampledCorrelations/DataRetrieval.jl

@@ -212,7 +212,7 @@ function broaden_energy(sc::SampledCorrelations, is, kernel::Function; negative_
     for (ω₀i, ω₀) in enumerate(ωvals)
         for (ωi, ω) in enumerate(ωvals)
             for qi in CartesianIndices(dims[1:end-1])
-                out[qi,ωi] += is[qi,ω₀i]*kernel(ω, ω₀)
+                out[qi,ωi] += is[qi,ω₀i]*kernel(ω, ω₀)*sc.Δω
             end
         end
     end

src/SampledCorrelations/SampledCorrelations.jl

@@ -16,8 +16,11 @@ struct SampledCorrelations{N}
 
     # 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
     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 
@@ -59,9 +62,11 @@ function clone_correlations(sc::SampledCorrelations{N}) where N
     # Avoid copies/deep copies of C-generated data structures
     space_fft! = 1/√(prod(dims)) * FFTW.plan_fft!(sc.samplebuf, (2,3,4))
     time_fft! = 1/√(n_all_ω) * FFTW.plan_fft!(sc.samplebuf, 6)
+    corr_fft! = FFTW.plan_fft!(sc.corrbuf, 4)
+    corr_ifft! = FFTW.plan_ifft!(sc.corrbuf, 4)
     M = isnothing(sc.M) ? nothing : copy(sc.M)
     return SampledCorrelations{N}(copy(sc.data), M, sc.crystal, sc.origin_crystal, sc.Δω,
-        deepcopy(sc.observables), copy(sc.samplebuf), space_fft!, time_fft!, sc.measperiod, sc.apply_g, sc.dt,
+        deepcopy(sc.observables), copy(sc.samplebuf), copy(sc.corrbuf), space_fft!, time_fft!, corr_fft!, corr_ifft!, sc.measperiod, sc.apply_g, sc.dt,
         copy(sc.nsamples), sc.processtraj!)
 end
 
@@ -166,6 +171,7 @@ function dynamical_correlations(sys::System{N}; dt=nothing, Δt=nothing, nω, ω
 
     # The sample buffer holds n_non_neg_ω measurements, and the rest is a zero buffer
     samplebuf = zeros(ComplexF64, num_observables(observables), sys.latsize..., na, n_all_ω)
+    corrbuf = zeros(ComplexF64, sys.latsize..., n_all_ω)
 
     # The output data has n_all_ω many (positive and negative and zero) frequencies
     data = zeros(ComplexF64, num_correlations(observables), na, na, sys.latsize..., n_all_ω)
@@ -176,14 +182,16 @@ function dynamical_correlations(sys::System{N}; dt=nothing, Δt=nothing, nω, ω
     # https://github.com/SunnySuite/Sunny.jl/issues/264 (subject to change).
     space_fft! = 1/√(prod(sys.latsize)) * FFTW.plan_fft!(samplebuf, (2,3,4))
     time_fft! = 1/√(n_all_ω) * FFTW.plan_fft!(samplebuf, 6)
+    corr_fft! = FFTW.plan_fft!(corrbuf, 4)
+    corr_ifft! = FFTW.plan_ifft!(corrbuf, 4)
 
     # Other initialization
     nsamples = Int64[0]
 
     # Make Structure factor and add an initial sample
     origin_crystal = isnothing(sys.origin) ? nothing : sys.origin.crystal
     sc = SampledCorrelations{N}(data, M, sys.crystal, origin_crystal, Δω, observables,
-                                samplebuf, space_fft!, time_fft!, measperiod, apply_g, dt, nsamples, process_trajectory)
+                                samplebuf, corrbuf, space_fft!, time_fft!, corr_fft!, corr_ifft!, measperiod, apply_g, dt, nsamples, process_trajectory)
 
     return sc
 end

test/test_correlation_sampling.jl

@@ -34,6 +34,7 @@
     S = spin_matrices(1/2)
     observables = Dict(:Sx => S[1], :Sy => S[2], :Sz => S[3])
     sc = dynamical_correlations(sys; nω=100, ωmax=10.0, dt=0.1, apply_g=false, observables)
+    Δω = sc.Δω
     add_sample!(sc, sys)
     qgrid = available_wave_vectors(sc)
     formula = intensity_formula(sc,:trace)
@@ -42,7 +43,7 @@
     # factor is already included in the Sqw intensity data. This convention will
     # change per https://github.com/SunnySuite/Sunny.jl/issues/264
     expected_sum_rule = prod(sys.latsize) * Sunny.norm2(sys.dipoles[1]) # NS^2 classical sum rule
-    @test isapprox(sum(Sqw), expected_sum_rule; atol=1e-12)
+    @test isapprox(sum(Sqw) * Δω, expected_sum_rule; atol=1e-12)
 
     # Test sum rule with default observables in dipole mode 
     sys = simple_model_fcc(; mode=:dipole)
@@ -53,7 +54,7 @@
     Sqw = intensities_interpolated(sc, qgrid, trace_formula; negative_energies=true)
     total_intensity_trace = sum(Sqw)
     expected_sum_rule = prod(sys.latsize) * Sunny.norm2(sys.dipoles[1]) # NS^2 classical sum rule
-    @test isapprox(sum(Sqw), expected_sum_rule; atol=1e-12)
+    @test isapprox(sum(Sqw) * Δω, expected_sum_rule; atol=1e-12)
 
     # Test perp reduces intensity
     perp_formula = intensity_formula(sc,:perp)
@@ -155,6 +156,6 @@ end
     # println(round.(data; digits=10))
 
     # Compare with reference
-    data_golden = [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/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]
     @test data ≈ data_golden
 end