diff --git a/docs/src/versions.md b/docs/src/versions.md index f983ac2d3..8e8d53763 100644 --- a/docs/src/versions.md +++ b/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) diff --git a/examples/04_GSD_FeI2.jl b/examples/04_GSD_FeI2.jl index a6569b476..712cbe36b 100644 --- a/examples/04_GSD_FeI2.jl +++ b/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, 0.5e-3), ) fig @@ -475,4 +475,4 @@ hm = heatmap(as, bs, is_static; ) ) Colorbar(hm.figure[1,2], hm.plot) -hm +hm \ No newline at end of file diff --git a/examples/08_Scattering_conventions.jl b/examples/08_Momentum_Conventions.jl similarity index 99% rename from examples/08_Scattering_conventions.jl rename to examples/08_Momentum_Conventions.jl index 32c01ff9c..c0e8342a1 100644 --- a/examples/08_Scattering_conventions.jl +++ b/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 diff --git a/src/Intensities/Interpolation.jl b/src/Intensities/Interpolation.jl index 4ef8ba927..e0b972e96 100644 --- a/src/Intensities/Interpolation.jl +++ b/src/Intensities/Interpolation.jl @@ -177,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 diff --git a/src/SampledCorrelations/CorrelationSampling.jl b/src/SampledCorrelations/CorrelationSampling.jl index 6c2fd12fd..4eef5908f 100644 --- a/src/SampledCorrelations/CorrelationSampling.jl +++ b/src/SampledCorrelations/CorrelationSampling.jl @@ -60,20 +60,12 @@ 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 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 +114,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 +143,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 diff --git a/src/SampledCorrelations/DataRetrieval.jl b/src/SampledCorrelations/DataRetrieval.jl index 3fbfa1336..ca1b99bd5 100644 --- a/src/SampledCorrelations/DataRetrieval.jl +++ b/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 diff --git a/src/SampledCorrelations/SampledCorrelations.jl b/src/SampledCorrelations/SampledCorrelations.jl index f1badd4da..d4f31fe13 100644 --- a/src/SampledCorrelations/SampledCorrelations.jl +++ b/src/SampledCorrelations/SampledCorrelations.jl @@ -7,7 +7,7 @@ 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 @@ -15,9 +15,12 @@ struct SampledCorrelations{N} observables :: ObservableInfo # Specs for sample generation and accumulation - samplebuf :: Array{ComplexF64, 6} # New sample buffer - space_fft! :: FFTW.AbstractFFTs.Plan # Pre-planned FFT - time_fft! :: 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 @@ -55,13 +58,14 @@ 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)) - time_fft! = 1/√(n_all_ω) * FFTW.plan_fft!(sc.samplebuf, 6) + 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) 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,16 +170,22 @@ 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_ω) M = calculate_errors ? zeros(Float64, size(data)...) : nothing - # The normalization is defined so that structure factor estimates of form - # ifft(fft * fft) carry an overall scaling factor consistent with this spec: - # 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) + # The normalization is defined so that the prod(sys.latsize)-many estimates + # of the structure factor produced by the correlation conj(space_fft!) * space_fft! + # are correctly averaged over. The corresponding time-average can't be applied in + # 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)) + time_fft! = FFTW.plan_fft!(samplebuf, 6) + corr_fft! = FFTW.plan_fft!(corrbuf, 4) + corr_ifft! = FFTW.plan_ifft!(corrbuf, 4) # Other initialization nsamples = Int64[0] @@ -183,7 +193,7 @@ function dynamical_correlations(sys::System{N}; dt=nothing, Δt=nothing, nω, ω # 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 diff --git a/test/test_binning.jl b/test/test_binning.jl index 570aebd66..cb5a4d835 100644 --- a/test/test_binning.jl +++ b/test/test_binning.jl @@ -61,14 +61,14 @@ # TODO: Test broadening formula = intensity_formula(sc, :perp) is, counts = intensities_binned(sc, params, formula) - is_golden = [0.2452071781; 0.0864959953; 0.1158561543; 0.0299947084;;;; 0.1226035891; 0.0432479977; 0.0579280772; 0.0149973542;;;; -0.0; -0.0; 0.0; -0.0] + is_golden = size(sc.samplebuf,6) * [0.2452071781; 0.0864959953; 0.1158561543; 0.0299947084;;;; 0.1226035891; 0.0432479977; 0.0579280772; 0.0149973542;;;; -0.0; -0.0; 0.0; -0.0] @test is ≈ is_golden @test all(counts .== 1.) # SQTODO: need to implement correct bin widths in powder_average_binned! is, counts = powder_average_binned(sc, (0,6π,6π/4), formula) - is_golden = [8.9511865548 4.4755932774 -0.0; 35.9054210445 17.9527105222 -0.0; 102.2777600371 51.1388800185 -0.0; 91.4466208136 45.7233104068 -0.0] + is_golden = size(sc.samplebuf,6) * [8.9511865548 4.4755932774 -0.0; 35.9054210445 17.9527105222 -0.0; 102.2777600371 51.1388800185 -0.0; 91.4466208136 45.7233104068 -0.0] counts_golden = [3.0 3.0 3.0; 15.0 15.0 15.0; 28.0 28.0 28.0; 39.0 39.0 39.0] @test is ≈ is_golden @test counts == counts_golden @@ -79,12 +79,12 @@ golden_correlations = [(sz,sz),(sy,sz),(sx,sz),(sy,sy),(sx,sy),(sx,sx)] formula = intensity_formula((k,ω,c) -> c,sc,golden_correlations; kT = 4.7, formfactors = [FormFactor("Fe2")], return_type = StaticArrays.SVector{6,ComplexF64}) is, counts = intensities_binned(sc, params, formula) - is_golden = [0.0180419588 -0.0; -0.0 -0.0; -0.0225463643 0.0; 0.0 0.0; -0.0384472716 0.0; 0.0 -0.0; 0.0281753523 -0.0; -0.0 0.0; 0.0480461242 0.0; -0.0 -0.0; 0.0819308319 -0.0; -0.0 0.0; 0.0235568126 -0.0; 0.0 0.0; -0.00045821 0.0; -0.0150200679 -0.0; -0.0069056941 -0.0; 0.0143894403 -0.0; 0.0095858637 -0.0; -0.0 -0.0; -0.0090405317 -0.0; -0.0046830351 -0.0; 0.0108140522 -0.0; -0.0 0.0] + is_golden = size(sc.samplebuf,6) * [0.0180419588 -0.0; -0.0 -0.0; -0.0225463643 0.0; 0.0 0.0; -0.0384472716 0.0; 0.0 -0.0; 0.0281753523 -0.0; -0.0 0.0; 0.0480461242 0.0; -0.0 -0.0; 0.0819308319 -0.0; -0.0 0.0; 0.0235568126 -0.0; 0.0 0.0; -0.00045821 0.0; -0.0150200679 -0.0; -0.0069056941 -0.0; 0.0143894403 -0.0; 0.0095858637 -0.0; -0.0 -0.0; -0.0090405317 -0.0; -0.0046830351 -0.0; 0.0108140522 -0.0; -0.0 0.0] @test reinterpret(Float64, is[1:2, 1, 1, 2:3]) ≈ is_golden @test all(counts .== 1.) is, counts = intensities_binned(sc, params, formula; integrated_kernel=integrated_lorentzian(fwhm=1.0)) - is_golden = [0.0112630945 0.005987424; -0.0 -0.0; -0.0140750699 -0.0074822609; 0.0 0.0; -0.0240015653 -0.0127591532; 0.0 0.0; 0.0175890909 0.0093503029; -0.0 -0.0; 0.0299938627 0.0159446388; -0.0 -0.0; 0.0511471459 0.0271896546; -0.0 -0.0; 0.0147058647 0.0078175893; 0.0 0.0; -0.0002860478 -0.0001520621; -0.0093766118 -0.004984576; -0.0043110333 -0.0022917311; 0.0089829284 0.0047752952; 0.0059841889 0.0031811751; -0.0 -0.0; -0.0056437532 -0.0030002006; -0.0029234889 -0.0015541171; 0.0067509129 0.0035887631; -0.0 -0.0] + is_golden = size(sc.samplebuf,6) * [0.0112630945 0.005987424; -0.0 -0.0; -0.0140750699 -0.0074822609; 0.0 0.0; -0.0240015653 -0.0127591532; 0.0 0.0; 0.0175890909 0.0093503029; -0.0 -0.0; 0.0299938627 0.0159446388; -0.0 -0.0; 0.0511471459 0.0271896546; -0.0 -0.0; 0.0147058647 0.0078175893; 0.0 0.0; -0.0002860478 -0.0001520621; -0.0093766118 -0.004984576; -0.0043110333 -0.0022917311; 0.0089829284 0.0047752952; 0.0059841889 0.0031811751; -0.0 -0.0; -0.0056437532 -0.0030002006; -0.0029234889 -0.0015541171; 0.0067509129 0.0035887631; -0.0 -0.0] counts_golden = [0.7164129516; 0.7164129516; 0.7164129516; 0.7164129516;;;; 0.6814242206; 0.6814242206; 0.6814242206; 0.6814242206;;;; 0.5437318669; 0.5437318669; 0.5437318669; 0.5437318669] @test reinterpret(Float64, is[1:2, 1, 1, 2:3]) ≈ is_golden @test counts ≈ counts_golden @@ -92,7 +92,7 @@ # Test all components using :full formula = intensity_formula(sc, :full; kT = 4.7, formfactors = [FormFactor("Fe2")]) is, counts = intensities_binned(sc, params, formula) - is2_golden = [0.0206923266 + 0.0im -0.0172987545 - 0.0089608308im -0.0132138143 + 0.0275337118im; -0.0172987545 + 0.0089608308im 0.0183422291 + 0.0im -0.0008767696 - 0.0287403967im; -0.0132138143 - 0.0275337118im -0.0008767696 + 0.0287403967im 0.0450751717 + 0.0im] + is2_golden = size(sc.samplebuf,6) * [0.0206923266 + 0.0im -0.0172987545 - 0.0089608308im -0.0132138143 + 0.0275337118im; -0.0172987545 + 0.0089608308im 0.0183422291 + 0.0im -0.0008767696 - 0.0287403967im; -0.0132138143 - 0.0275337118im -0.0008767696 + 0.0287403967im 0.0450751717 + 0.0im] @test is[2] ≈ is2_golden @test all(counts .== 1.) diff --git a/test/test_correlation_sampling.jl b/test/test_correlation_sampling.jl index e5e1780b2..0ff1202d4 100644 --- a/test/test_correlation_sampling.jl +++ b/test/test_correlation_sampling.jl @@ -34,15 +34,17 @@ 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) + Δq³ = 1/prod(sys.latsize) # Fraction of a BZ formula = intensity_formula(sc,:trace) Sqw = intensities_interpolated(sc, qgrid, formula; negative_energies=true) - # A simple sum (rather than energy-integral) is possible because a `Δω` - # 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) + expected_sum_rule = Sunny.norm2(sys.dipoles[1]) # S^2 classical sum rule + @test isapprox(sum(Sqw) * Δq³ * Δω, expected_sum_rule; atol=1e-12) + params = unit_resolution_binning_parameters(sc;negative_energies=true) + Sqw_integrated, counts = intensities_binned(sc, params, formula) + @test isapprox(sum(Sqw_integrated), expected_sum_rule; atol=1e-12) # Test sum rule with default observables in dipole mode sys = simple_model_fcc(; mode=:dipole) @@ -52,8 +54,12 @@ trace_formula = intensity_formula(sc,:trace) 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) + expected_sum_rule = Sunny.norm2(sys.dipoles[1]) # S^2 classical sum rule + @test isapprox(sum(Sqw) * Δq³ * Δω, expected_sum_rule; atol=1e-12) + # Test binned version doesn't require ΔqΔω measure + params = unit_resolution_binning_parameters(sc;negative_energies=true) + Sqw_integrated, counts = intensities_binned(sc, params, trace_formula) + @test isapprox(sum(Sqw_integrated), expected_sum_rule; atol=1e-12) # Test perp reduces intensity perp_formula = intensity_formula(sc,:perp) @@ -155,6 +161,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.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 diff --git a/test/test_intensities_setup.jl b/test/test_intensities_setup.jl index 7206db3fa..b13dc1d93 100644 --- a/test/test_intensities_setup.jl +++ b/test/test_intensities_setup.jl @@ -95,11 +95,16 @@ end sc = dynamical_correlations(sys; dt=1, nω=3, ωmax=1) add_sample!(sc, sys) - # A simple sum (rather than energy-integral) is possible because a `Δω` - # factor is already included in the intensity data. This convention will - # change per https://github.com/SunnySuite/Sunny.jl/issues/264 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] + 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 + 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 @@ -108,10 +113,10 @@ end # need to include the g factor. This is the equal-space-and-time correlation value: gS_squared = (2 * 1/2)^2 - expected_sum = gS_squared * prod(sys.latsize) + 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 @@ -132,5 +137,5 @@ end is, counts = intensities_binned(sc,params_pasr,formula) expected_multi_BZ_sum = gS_squared * prod(nbzs) * nfs # ⟨S⋅S⟩ expected_multi_BZ_sum_times_natoms = expected_multi_BZ_sum * Sunny.natoms(sc.crystal) # Nₐ×⟨S⋅S⟩ - @test sum(is ./ counts) ≈ expected_multi_BZ_sum_times_natoms / size(sc.data,7) + @test sum(is) ≈ expected_multi_BZ_sum_times_natoms end