Skip to content

Commit

Permalink
Classical intensities as density (#272)
Browse files Browse the repository at this point in the history 
The result of `intensities_interpolated` is now a proper density in both q and ω variables. Also pre-plan some FFTs and fix color scales in example heatmaps.

---------

Co-authored-by: Kipton Barros <kbarros@gmail.com>
Co-authored-by: Sam Quinn <samquinn0451@gmail.com>
  • Loading branch information
@ddahlbom @kbarros @Lazersmoke
3 people authored Jun 7, 2024
1 parent 290a3ff commit b169d68
Show file tree
Hide file tree
Showing 10 changed files with 90 additions and 64 deletions.
15 changes: 7 additions & 8 deletions docs/src/versions.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
8 changes: 4 additions & 4 deletions examples/04_GSD_FeI2.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -475,4 +475,4 @@ hm = heatmap(as, bs, is_static;
)
)
Colorbar(hm.figure[1,2], hm.plot)
hm
hm
2 changes: 1 addition & 1 deletion examples/08_Scattering_conventions.jl → examples/08_Momentum_Conventions.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
13 changes: 13 additions & 0 deletions src/Intensities/Interpolation.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
25 changes: 9 additions & 16 deletions src/SampledCorrelations/CorrelationSampling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/SampledCorrelations/DataRetrieval.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
38 changes: 24 additions & 14 deletions src/SampledCorrelations/SampledCorrelations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,17 +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
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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -166,24 +170,30 @@ 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]

# 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
Expand Down
10 changes: 5 additions & 5 deletions test/test_binning.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -79,20 +79,20 @@
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

# 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.)

Expand Down
22 changes: 14 additions & 8 deletions test/test_correlation_sampling.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Loading

0 comments on commit b169d68

Please sign in to comment.