Candidate main

Project.toml

@@ -4,6 +4,7 @@ authors = ["The Sunny team"]
 version = "0.4.4"
 
 [deps]
+CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CrystalInfoFramework = "6007d9b0-c6b2-11e8-0510-1d10e825f3f1"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -12,8 +13,10 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FilePathsBase = "48062228-2e41-5def-b9a4-89aafe57970f"
 Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
+JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearMaps = "7a12625a-238d-50fd-b39a-03d52299707e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"

docs/make.jl

@@ -4,7 +4,7 @@ using Literate, Documenter, Sunny
 
 execute = true # set `false` to disable cell evaluation
 
-example_names = ["fei2_tutorial", "powder_averaging", "ising2d"]
+example_names = ["fei2_tutorial", "powder_averaging", "ising2d", "binning_tutorial"]
 example_sources = [joinpath(@__DIR__, "..", "examples", name*".jl") for name in example_names]
 example_destination = joinpath(@__DIR__, "src", "examples")
 example_doc_paths = ["examples/$name.md" for name in example_names]

docs/src/structure-factor.md

@@ -134,15 +134,15 @@ the dynamic case. As was true there, it is important to ensure that the spins in
 ### Extracting information from structure factors
 
 The basic function for extracting information from a dynamic `StructureFactor`
-at a particular wave vector, $𝐪$, is [`intensities`](@ref). It takes a
+at a particular wave vector, $𝐪$, is [`intensities_interpolated`](@ref). It takes a
 `StructureFactor`, a list of wave vectors, and a contraction mode. For example,
-`intensities(sf, [[0.0, 0.5, 0.5]], :trace)` will calculate intensities for the
-wavevector $𝐪 = (𝐛_2 + 𝐛_3)/2$. The option `:trace` will contract spin
-indices, returning $𝒮^{αα}(𝐪,ω)$, summing over ``α``. The option `:perp` will
-instead perform a contraction that includes polarization corrections. The option
-`:full` will return data for the full tensor $𝒮^{αβ}(𝐪,ω)$. `intensities`
-returns a list of `nω` elements. The corresponding $ω$ values can be retrieved
-by calling [`ωs`](@ref).
+`intensities_interpolated(sf, [[0.0, 0.5, 0.5]])` will calculate intensities for the
+wavevector $𝐪 = (𝐛_2 + 𝐛_3)/2$. The keyword argument `formula` can be used to
+specify an [`intensity_formula`](@ref) for greater control over the intensity calculation.
+The default formula performs a contraction of $𝒮^{αβ}(𝐪,ω)$ that includes
+polarization corrections. `intensities_interpolated`
+returns a list of `nω` elements at each wavevector. The corresponding $ω$ values can be retrieved
+by calling [`ωs`](@ref) on `sf`.
 
 Since Sunny currently only calculates the structure factor on a finite lattice,
 it is important to realize that exact information is only available at a
@@ -151,32 +151,35 @@ information at $q_i = \frac{n}{L_i}$, where $n$ runs from $(\frac{-L_i}{2}+1)$
 to $\frac{L_i}{2}$ and $L_i$ is the linear dimension of the lattice used for the
 calculation. If you request a wave vector that does not fall into this set,
 Sunny will automatically round to the nearest $𝐪$ that is available. If
-`intensities` is given the keyword argument `interpolation=:linear`, Sunny will
-use trilinear interpolation to determine the results at the requested wave
+`intensities_interpolated` is given the keyword argument `interpolation=:linear`, Sunny will
+use trilinear interpolation to determine a result at the requested wave
 vector. 
 
 To retrieve the intensities at all wave vectors for which there is exact data,
 first call the function [`all_exact_wave_vectors`](@ref) to generate a list of
 `qs`. This takes an optional keyword argument `bzsize`, which must be given a
 tuple of three integers specifying the number of Brillouin zones to calculate,
 e.g., `bzsize=(2,2,2)`. The resulting list of wave vectors may then be passed to
-`intensities`.
+`intensities_interpolated`.
+
+Alternatively, [`intensities_binned`](@ref) can be used to place the exact data
+into histogram bins for comparison with experiment.
 
 The convenience function [`connected_path`](@ref) returns a list of wavevectors
 sampled along a path that connects specified $𝐪$ points. This list can be used
 as an input to `intensities`. Another convenience method,
 [`spherical_shell`](@ref) will provide a list of wave vectors on a sphere of a
 specified radius. This is useful for powder averaging. 
 
-A number of keyword arguments are available which modify the calculation of
-structure factor intensity. See the documentation of [`intensities`](@ref) for a
-full list. It is generally recommended to provide a value to `kT` corresponding
-to the temperature of sampled configurations. Given `kT`, Sunny will apply an
-energy- and temperature-dependent classical-to-quantum rescaling of intensities. 
+A number of arguments for [`intensity_formula`](@ref) are available which
+modify the calculation of structure factor intensity. It is generally recommended
+to provide a value of `kT` corresponding to the temperature of sampled configurations.
+Given `kT`, Sunny will include an energy- and temperature-dependent classical-to-quantum 
+rescaling of intensities in the formula.
 
 To retrieve intensity data from a instantaneous structure factor, use
-[`instant_intensities`](@ref), which shares keyword arguments with
-`intensities`. This function may also be used to calculate instantaneous
+[`instant_intensities_interpolated`](@ref), which accepts similar arguments to
+`intensities_interpolated`. This function may also be used to calculate instantaneous
 information from a dynamical structure factor. Note that it is important to
 supply a value to `kT` to reap the benefits of this approach over simply
-calculating a static structure factor at the outset. 
+calculating a static structure factor at the outset. 

examples/binning_tutorial.jl

@@ -0,0 +1,175 @@
+# # Binning Tutorial
+
+using Sunny, GLMakie
+
+# We specify the crystal lattice structure of CTFD using the lattice parameters specified by
+#
+# W Wan et al 2020 J. Phys.: Condens. Matter 32 374007 DOI 10.1088/1361-648X/ab757a
+latvecs = lattice_vectors(8.113,8.119,12.45,90,100,90)
+positions = [[0,0,0]]
+types = ["Cu"]
+formfactors  = [FormFactor(1,"Cu2")]
+xtal = Crystal(latvecs,positions;types);
+
+# We will use a somewhat small periodic lattice size of 6x6x4 in order
+# to showcase the effect of a finite lattice size.
+latsize = (6,6,4);
+
+# In this system, the magnetic lattice is the same as the chemical lattice,
+# and there is a spin-1/2 dipole on each site.
+magxtal = xtal;
+valS = 1/2;
+sys = System(magxtal, latsize, [SpinInfo(1;S = valS)], :dipole; seed=1);
+
+# Quoted value of J = +6.19(2) meV (antiferromagnetic) between
+# nearest neighbors on the square lattice
+J = 6.19 # meV
+characteristic_energy_scale = abs(J * valS)
+set_exchange!(sys,J,Bond(1,1,[1,0,0]))
+set_exchange!(sys,J,Bond(1,1,[0,1,0]))
+
+# Thermalize the system using the Langevin intergator.
+# The timestep and the temperature are roughly based off the characteristic
+# energy scale of the problem.
+Δt = 0.05 / characteristic_energy_scale
+kT = 0.01 * characteristic_energy_scale
+langevin = Langevin(Δt; λ=0.1, kT=kT)
+randomize_spins!(sys);
+for i in 1:10_000 # Long enough to reach equilibrium
+    step!(sys, langevin)
+end 
+
+# The neutron spectrometer used in the experiment had an incident neutron energy of 36.25 meV.
+# Since this is the most amount of energy that can be deposited by the neutron into the sample, we
+# don't need to consider energies higher than this.
+ωmax = 36.25;
+
+# We choose the resolution in energy (specified by the number of nω modes resolved)
+# to be ≈20× better than the experimental resolution in order to demonstrate
+# the effect of over-resolving in energy.
+nω = 480; 
+dsf = DynamicStructureFactor(sys; Δt=Δt, nω=nω, ωmax=ωmax, process_trajectory=:symmetrize)
+
+# We re-sample from the thermal equilibrium distribution several times to increase our sample size
+nsamples = 3
+for _ in 1:nsamples
+    for _ in 1:8000 
+        step!(sys, langevin)
+    end
+    add_sample!(dsf, sys)
+end
+
+# Since the SU(N)NY crystal has only finitely many lattice sites, there are finitely
+# many ways for a neutron to scatter off of the sample.
+# We can visualize this discreteness by plotting each possible Qx and Qz, for example:
+
+lower_aabb_q, upper_aabb_q = Sunny.binning_parameters_aabb(unit_resolution_binning_parameters(dsf))#hide
+lower_aabb_cell = floor.(Int64,lower_aabb_q .* latsize .+ 1)#hide
+upper_aabb_cell = ceil.(Int64,upper_aabb_q .* latsize .+ 1)#hide
+
+Qx = zeros(Float64,0)#hide
+Qz = zeros(Float64,0)#hide
+for cell in CartesianIndices(Tuple(((:).(lower_aabb_cell,upper_aabb_cell))))#hide
+    q = (cell.I .- 1) ./ latsize # q is in R.L.U.#hide
+    push!(Qx,q[1])#hide
+    push!(Qz,q[3])#hide
+end#hide
+fig = Figure()#hide
+ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Qz [r.l.u.]")#hide
+## Compute some scattering vectors at and around the first BZ...
+scatter!(ax,Qx,Qz)
+fig#hide
+
+# One way to display the structure factor is to create a histogram with
+# one bin centered at each discrete scattering possibility using [`unit_resolution_binning_parameters`](@ref) to create a set of [`BinningParameters`](@ref).
+params = unit_resolution_binning_parameters(dsf)
+
+# Since this is a 4D histogram,
+# it further has to be integrated over two of those directions in order to be displayed.
+# Here, we integrate over Qy and Energy using [`integrate_axes!`](@ref):
+integrate_axes!(params;axes = [2,4]) # Integrate over Qy (2) and E (4)
+
+# Now that we have parameterized the histogram, we can bin our data.
+# In addition to the [`BinningParameters`](@ref), an [`intensity_formula`](@ref) needs to be
+# provided to specify which dipole, temperature, and atomic form factor
+# corrections should be applied during the intensity calculation.
+formula = intensity_formula(dsf, :perp; kT, formfactors)
+intensity,counts = intensities_binned(dsf, params; formula)
+normalized_intensity = intensity ./ counts;
+
+# With the data binned, we can now plot it. The axes labels give the bin centers of each bin, as given by [`axes_bincenters`](@ref).
+function plot_data(params) #hide
+intensity,counts = intensities_binned(dsf, params; formula)#hide
+normalized_intensity = intensity ./ counts;#hide
+bin_centers = axes_bincenters(params);
+
+fig = Figure()#hide
+ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Qz [r.l.u.]")#hide
+heatmap!(ax,bin_centers[1],bin_centers[3],normalized_intensity[:,1,:,1])
+scatter!(ax,Qx,Qz)
+xlims!(ax,params.binstart[1],params.binend[1])
+ylims!(ax,params.binstart[3],params.binend[3])
+return fig#hide
+
+end#hide
+
+plot_data(params)#hide
+
+# Note that some bins have no scattering vectors at all when the bin size is made too small:
+params = unit_resolution_binning_parameters(dsf)#hide
+integrate_axes!(params;axes = [2,4])#hide
+params.binwidth[1] /= 1.2
+params.binwidth[3] /= 2.5
+plot_data(params)#hide
+
+# Conversely, making the bins bigger doesn't break anything, but loses resolution:
+params = unit_resolution_binning_parameters(dsf)#hide
+integrate_axes!(params;axes = [2,4])#hide
+params.binwidth[1] *= 2
+params.binwidth[3] *= 2
+plot_data(params)#hide
+
+# Recall that while we under-resolved in Q by choosing a small lattice, we 
+# over-resolved in energy:
+x = zeros(Float64,0)#hide
+y = zeros(Float64,0)#hide
+for omega = ωs(dsf), qx = unique(Qx)#hide
+    push!(x,qx)#hide
+    push!(y,omega)#hide
+end#hide
+ax = Axis(fig[1,1];xlabel="Qx [r.l.u]",ylabel="Energy [meV]")#hide
+scatter!(ax,x,y)
+
+# Let's make a new histogram which includes the energy axis.
+# The x-axis of the histogram will be a 1D cut from `Q = [0,0,0]` to `Q = [1,1,0]`.
+# See [`one_dimensional_cut_binning_parameters`](@ref).
+x_axis_bin_count = 10
+cut_width = 0.3
+params = one_dimensional_cut_binning_parameters(dsf,[0,0,0],[1,1,0],x_axis_bin_count,cut_width)
+
+# There are no longer any scattering vectors exactly in the plane of the cut. Instead,
+# as described in the `BinningParameters` output above, the transverse Q
+# directions are integrated over, so slightly out of plane points are included.
+#
+# We plot the intensity on a log-scale to improve visibility.
+intensity,counts = intensities_binned(dsf, params; formula)
+log_intensity = log10.(intensity ./ counts);
+bin_centers = axes_bincenters(params);#hide
+fig = Figure()#hide
+ax = Axis(fig[1,1];xlabel="Progress along cut [r.l.u]",ylabel="Energy [meV]")#hide
+heatmap!(ax,bin_centers[1],bin_centers[4],log_intensity[:,1,1,:])
+xlims!(ax,params.binstart[1],params.binend[1])#hide
+ylims!(ax,params.binstart[4],params.binend[4])#hide
+fig#hide
+
+# By reducing the number of energy bins to be closer to the number of bins on the x-axis, we can make the dispersion curve look nicer:
+params.binwidth[4] *= 20
+intensity,counts = intensities_binned(dsf, params; formula)#hide
+log_intensity = log10.(intensity ./ counts);#hide
+bin_centers = axes_bincenters(params);#hide
+fig = Figure()#hide
+ax = Axis(fig[1,1];xlabel="Progress along cut [r.l.u]",ylabel="Energy [meV]")#hide
+heatmap!(ax,bin_centers[1],bin_centers[4],log_intensity[:,1,1,:])
+xlims!(ax,params.binstart[1],params.binend[1])#hide
+ylims!(ax,params.binstart[4],params.binend[4])#hide
+fig#hide