Fix intensity scale for SampledCorrelations with reshaped system

examples/spinw_tutorials/SW18_Distorted_kagome.jl

@@ -24,8 +24,8 @@ view_crystal(cryst)
 
 # Define the interactions.
 
-spininfos = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
-sys = System(cryst, spininfos, :dipole, seed=0)
+moments = [1 => Moment(s=1/2, g=2), 3 => Moment(s=1/2, g=2)]
+sys = System(cryst, moments, :dipole, seed=0)
 J   = -2
 Jp  = -1
 Jab = 0.75

src/Operators/Rotation.jl

@@ -118,7 +118,8 @@ function unitary_for_rotation(R::Mat3, gen)
     return exp(-im*θ*(n'*gen))
 end
 
-# Unitary for a rotation matrix in the N-dimensional irrep of SU(2).
+# Unitary for a rotation matrix in the N-dimensional irrep of SU(2). TODO: Use
+# polynomial expansions for speed: http://www.emis.de/journals/SIGMA/2014/084/.
 function unitary_irrep_for_rotation(R::Mat3; N::Int)
     gen = spin_matrices_of_dim(; N)
     unitary_for_rotation(R, gen)

src/SampledCorrelations/CorrelationUtils.jl

@@ -1,13 +1,15 @@
 """
-    available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1))
+    available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1))
 
-Returns all wave vectors for which `sc` contains exact values. `bsize` specifies
-the number of Brillouin zones to be included.
+Returns the grid of wave vectors for which `sc` contains exact values.
+Optionally extend by a given number of `counts` along each grid axis. If the
+system was not reshaped, then the number of Brillouin zones included is
+`prod(counts)`.
 """
-function available_wave_vectors(sc::SampledCorrelations; bzsize=(1,1,1))
+function available_wave_vectors(sc::SampledCorrelations; counts=(1,1,1))
     Ls = sc.sys_dims
     offsets = map(L -> isodd(L) ? 1 : 0, Ls)
-    up = Ls .* bzsize
+    up = Ls .* counts
     hi = map(L -> L - div(L, 2), up) .- offsets
     lo = map(L -> 1 - div(L, 2), up) .- offsets
 

src/SampledCorrelations/DataRetrieval.jl

@@ -70,7 +70,8 @@ end
 
 contains_dynamic_correlations(sc) = !isnan(sc.Δω)
 
-# Documented under intensities function for LSWT.
+# Documented under intensities function for LSWT. TODO: As a hack, this function
+# is also being used as the back-end to intensities_static.
 function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT)
     if !isnothing(kT) && kT <= 0
         error("Positive `kT` required for classical-to-quantum corrections, or set `kT=nothing` to disable.")
@@ -92,16 +93,21 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
 
     # Prepare memory and configuration variables for actual calculation
     qpts = Base.convert(AbstractQPoints, qpts)
+    qs_reshaped = [to_reshaped_rlu(sc, q) for q in qpts.qs]
+
     ffs = sc.measure.formfactors[1, :] # FIXME
     intensities = zeros(eltype(sc.measure), isnan(sc.Δω) ? 1 : length(ωs), length(qpts.qs)) # N.B.: Inefficient indexing order to mimic LSWT
-    q_idx_info = pruned_wave_vector_info(sc, qpts.qs)
+    q_idx_info = pruned_wave_vector_info(sc, qs_reshaped)
     crystal = @something sc.origin_crystal sc.crystal
     NCorr  = Val{size(sc.data, 1)}()
     NAtoms = Val{size(sc.data, 2)}()
 
     # Intensities calculation
     intensities_rounded!(intensities, sc.data, sc.crystal, sc.measure, ffs, q_idx_info, ωidcs, NCorr, NAtoms)
 
+    # Convert to a q-space density in original (not reshaped) RLU.
+    intensities .*= det(sc.crystal.recipvecs) / det(crystal.recipvecs)
+
     # Post-processing steps for dynamical correlations 
     if contains_dynamic_correlations(sc) 
         # Convert time axis to a density.
@@ -120,7 +126,6 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
 
     intensities = reshape(intensities, length(ωs), size(qpts.qs)...)
 
-    # TODO: Refactor this logic so that return value is uniform.
     return if contains_dynamic_correlations(sc) 
         Intensities(crystal, qpts, collect(ωs), intensities)
     else
@@ -202,6 +207,3 @@ end
 =#
 
 
-
-
-

src/SampledCorrelations/SampledCorrelations.jl

@@ -3,7 +3,7 @@ mutable struct SampledCorrelations
     const data           :: Array{ComplexF64, 7}                 # Raw SF with sublattice indices (ncorrs × natoms × natoms × sys_dims × nω)
     const M              :: Union{Nothing, Array{Float64, 7}}    # Running estimate of (nsamples - 1)*σ² (where σ² is the variance of intensities)
     const crystal        :: Crystal                              # Crystal for interpretation of q indices in `data`
-    const origin_crystal :: Union{Nothing,Crystal}               # Original user-specified crystal (if different from above) -- needed for FormFactor accounting
+    const origin_crystal :: Union{Nothing,Crystal}               # Original user-specified crystal (if different from above)
     const Δω             :: Float64                              # Energy step size (could make this a virtual property)  
     measure              :: MeasureSpec                          # Observable, correlation pairs, and combiner
 

src/Sunny.jl

@@ -126,7 +126,7 @@ export BinningParameters, load_nxs, generate_mantid_script_from_binning_paramete
 include("deprecated.jl")
 export set_external_field!, set_external_field_at!, meV_per_K,
     dynamic_correlations, instant_correlations, intensity_formula, reciprocal_space_path,
-    set_spiral_order_on_sublattice!, set_spiral_order!
+    delta_function_kernel, set_spiral_order_on_sublattice!, set_spiral_order!
 
 isloaded(pkg::String) = any(k -> k.name == pkg, keys(Base.loaded_modules))
 

src/System/System.jl

@@ -227,6 +227,19 @@ An iterator over all [`Site`](@ref)s in the system.
 """
 @inline eachsite(sys::System) = CartesianIndices(sys.dipoles)
 
+"""
+nsites(sys::System) = length(eachsite(sys))
+"""
+nsites(sys::System) = length(eachsite(sys))
+
+"""
+    ncells(sys::System)
+
+Number of chemical cells in the system using the original crystal.
+"""
+ncells(sys::System) = nsites(sys) / natoms(orig_crystal(sys))
+
+
 """
     global_position(sys::System, site::Site)
 

src/deprecated.jl

@@ -85,7 +85,7 @@ Base.@deprecate instant_correlations(sys; opts...) let
 end
 
 Base.@deprecate intensity_formula(opts1...; opts2...) let
-    error("Formulas have been removed. Call `intensities` directly.")
+    error("Formulas have been removed. See revised Sunny docs for new interface.")
 end
 
 Base.@deprecate reciprocal_space_path(cryst::Crystal, qs, density) let
@@ -110,6 +110,9 @@ Base.@deprecate SpinInfo(i; S, g) let
     i => Moment(; s=S, g)
 end
 
+Base.@deprecate_binding delta_function_kernel nothing
+
+
 # REMEMBER TO ALSO DELETE:
 #
 # * view_crystal(cryst, max_dist)

test/test_correlation_sampling.jl

@@ -26,7 +26,6 @@
             step!(sys, langevin)
         end
     end
-
 
     # Test classical sum rule in SU(N) mode
     sys = simple_model_fcc(; mode=:SUN)
@@ -35,13 +34,13 @@
     sc = SampledCorrelations(sys; dt=0.08, energies, measure=ssf_trace(sys; apply_g=false))
     Δω = sc.Δω
     add_sample!(sc, sys)
-    qgrid = Sunny.QPoints(Sunny.available_wave_vectors(sc)[:])
+    qgrid = Sunny.available_wave_vectors(sc)[:]
     Δq³ = 1/prod(sys.dims) # Fraction of a BZ
     Sqw = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing)
     expected_sum_rule = Sunny.norm2(sys.dipoles[1]) # S^2 classical sum rule
     @test isapprox(sum(Sqw.data) * Δq³ * Δω, expected_sum_rule; atol=1e-12)
 
-    # Test classical sum rule in dipole mode 
+    # Test classical sum rule in dipole mode
     sys = simple_model_fcc(; mode=:dipole)
     thermalize_simple_model!(sys; kT=0.1)
     sc = SampledCorrelations(sys; dt=0.08, energies, measure=ssf_trace(sys; apply_g=false))
@@ -57,44 +56,41 @@
     total_intensity_unpolarized = sum(Sqw_perp.data)
     @test total_intensity_unpolarized < total_intensity_trace
 
-
     # Test diagonal elements are approximately real (at one wave vector)
     function ssf_diag_imag(sys::System{N}; apply_g=false) where N
         return ssf_custom(sys; apply_g) do _, ssf
             tr(imag(ssf))
         end
     end
     sc.measure = ssf_diag_imag(sys; apply_g=false)
-    is = intensities(sc, Sunny.QPoints([[0.25, 0.5, 0]]); energies=:available, kT=nothing)
-    is.data
-    @test sum(is.data) < 1e-14
-
+    res = intensities(sc, [[0.25, 0.5, 0]]; energies=:available, kT=nothing)
+    res.data
+    @test sum(res.data) < 1e-14
 
     # Test form factor correction works and is doing something.
     formfactors = [1 => FormFactor("Fe2")]
     sc.measure = ssf_trace(sys; apply_g=false, formfactors)
-    is = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing)
-    total_intensity_ff = sum(is.data)
+    res = intensities(sc, qgrid; energies=:available_with_negative, kT=nothing)
+    total_intensity_ff = sum(res.data)
     @test total_intensity_ff != total_intensity_trace
 
-
     # Test static from dynamic intensities working
     sc.measure = ssf_trace(sys; apply_g=false)
-    is_static = intensities_static(sc, qgrid; kT=nothing)
-    total_intensity_static = sum(is_static.data)
+    res_static = intensities_static(sc, qgrid; kT=nothing)
+    total_intensity_static = sum(res_static.data)
     @test isapprox(total_intensity_static, total_intensity_trace * sc.Δω; atol=1e-9)  # Order of summation can lead to very small discrepancies
 
     # Test quantum-to-classical increases intensity
-    is_static_c2q = intensities_static(sc, qgrid; kT=0.1)
-    total_intensity_static_c2q = sum(is_static_c2q.data)
-    @test total_intensity_static_c2q > total_intensity_static 
+    res_static_c2q = intensities_static(sc, qgrid; kT=0.1)
+    total_intensity_static_c2q = sum(res_static_c2q.data)
+    @test total_intensity_static_c2q > total_intensity_static
 
     # Test static intensities working
     sys = simple_model_fcc(; mode=:dipole)
     thermalize_simple_model!(sys; kT=0.1)
-    ic = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=false))
-    add_sample!(ic, sys)
-    true_static_vals = intensities_static(ic, qgrid)
+    res = SampledCorrelationsStatic(sys; measure=ssf_trace(sys; apply_g=false))
+    add_sample!(res, sys)
+    true_static_vals = intensities_static(res, qgrid)
     true_static_total = sum(true_static_vals.data)
     @test isapprox(true_static_total / prod(sys.dims), 1.0; atol=1e-12)
 end
@@ -106,7 +102,7 @@ end
     randomize_spins!(sys)
 
     # Set up Langevin sampler.
-    dt_langevin = 0.07 
+    dt_langevin = 0.07
     langevin = Langevin(dt_langevin; damping=0.1, kT=0.1723)
 
     measure = ssf_trace(sys)
@@ -141,10 +137,10 @@ end
 
     sc = SampledCorrelations(sys; energies=range(0.0, 5.5, 10), dt=0.12, measure=ssf_perp(sys))
     add_sample!(sc, sys)
-    qs = Sunny.QPoints([[0.0, 0.0, 0.0], [-0.2, 0.4, -0.1]])
-    is = intensities(sc, qs; energies=:available, kT=nothing)
+    qs = [[0.0, 0.0, 0.0], [-0.2, 0.4, -0.1]]
+    res = intensities(sc, qs; energies=:available, kT=nothing)
 
     # Compare with reference
-    data_golden = [33.52245944537883 31.523781055757002; 16.76122972268928 16.188214427928443; 1.6337100022968968e-14 5.3112747876921045; 1.590108516238891e-13 1.8852219123773621; -1.5194916032483394e-14 0.06400012935688847; -6.217248937900877e-15 -0.01803103943766904; 7.863406895322359e-14 -0.04445301974061088; 7.014086281484114e-14 0.0025512102338097653; 3.195591939212742e-14 -0.02515685630480813; 3.6201681383269686e-14 0.023924996100518413] 
-    @test is.data ≈ data_golden
+    data_golden = [33.52245944537883 31.523781055757002; 16.76122972268928 16.188214427928443; 1.6337100022968968e-14 5.3112747876921045; 1.590108516238891e-13 1.8852219123773621; -1.5194916032483394e-14 0.06400012935688847; -6.217248937900877e-15 -0.01803103943766904; 7.863406895322359e-14 -0.04445301974061088; 7.014086281484114e-14 0.0025512102338097653; 3.195591939212742e-14 -0.02515685630480813; 3.6201681383269686e-14 0.023924996100518413]
+    @test res.data ≈ data_golden
 end

test/test_ewald.jl

@@ -17,14 +17,14 @@
     latvecs = lattice_vectors(1,1,1,90,90,90)
     positions = [[0,0,0]]
     cryst = Crystal(latvecs, positions)
-    infos = [1 => Moment(s=1, g=1)]
-    sys = System(cryst, infos, :dipole)
+    moments = [1 => Moment(s=1, g=1)]
+    sys = System(cryst, moments, :dipole)
     enable_dipole_dipole!(sys, 1.0)
     @test ewalder_energy(sys) ≈ -1/6
     @test isapprox(energy(sys), -1/6; atol=1e-13)
 
     # Same thing, with multiple unit cells
-    sys = System(cryst, infos, :dipole; dims=(2, 3, 4))
+    sys = System(cryst, moments, :dipole; dims=(2, 3, 4))
     enable_dipole_dipole!(sys, 1.0)
     @test isapprox(energy_per_site(sys), -1/6; atol=1e-13)
 
@@ -33,12 +33,12 @@
     positions = [[0,0,0], [0.1,0,0], [0.6,0.4,0.5]]
     cryst = Crystal(latvecs, positions)
     Random.seed!(0) # Don't have sys.rng yet
-    infos = [
+    moments = [
         1 => Moment(s=1, g=rand(3,3)),
         2 => Moment(s=3/2, g=rand(3,3)),
         3 => Moment(s=2, g=rand(3,3)),
     ]
-    sys = System(cryst, infos, :dipole)
+    sys = System(cryst, moments, :dipole)
     enable_dipole_dipole!(sys, 1.0)
     randomize_spins!(sys)