Proposed modifications to the interface (#86)

* Proposed modifications to the interface

* Standardize pretty-printing of System and StructureFactor

* Add IntensityFormula types

* Fix things in docs without tests

Project.toml

@@ -14,6 +14,7 @@ Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 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/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

@@ -48,7 +48,7 @@ end
 # 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); 
+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
@@ -90,14 +90,16 @@ params = unit_resolution_binning_parameters(dsf)
 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.
-# The arguments to [`intensities_binned`](@ref) beyond `params` specify which dipole, temperature,
-# and atomic form factor corrections should be applied during the intensity calculation.
-intensity,counts = intensities_binned(dsf, params, :perp; kT, formfactors);
+# 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, :perp; kT, formfactors) #hide
+intensity,counts = intensities_binned(dsf, params; formula)#hide
 normalized_intensity = intensity ./ counts;#hide
 bin_centers = axes_bincenters(params);
 
@@ -150,7 +152,7 @@ params = one_dimensional_cut_binning_parameters(dsf,[0,0,0],[1,1,0],x_axis_bin_c
 # 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, :perp; kT, formfactors)
+intensity,counts = intensities_binned(dsf, params; formula)
 log_intensity = log10.(intensity ./ counts);
 bin_centers = axes_bincenters(params);#hide
 fig = Figure()#hide
@@ -162,7 +164,7 @@ 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, :perp; kT, formfactors)#hide
+intensity,counts = intensities_binned(dsf, params; formula)#hide
 log_intensity = log10.(intensity ./ counts);#hide
 bin_centers = axes_bincenters(params);#hide
 fig = Figure()#hide

examples/fei2_tutorial.jl

@@ -352,7 +352,7 @@ end
 # resolve. For the time step, twice the value used for the Langevin integrator
 # is usually a good choice.
 
-sf = DynamicStructureFactor(sys_large; Δt=2Δt, nω=120, ωmax=7.5);
+sf = DynamicStructureFactor(sys_large; Δt=2Δt, nω=120, ωmax=7.5)
 
 # `sf` currently contains dynamical structure data generated from a single
 # sample. Additional samples can be added by generating a new spin configuration
@@ -366,18 +366,21 @@ for _ in 1:2
 end
 
 # ## Accessing structure factor data 
-# The basic function for accessing intensity data is [`intensities`](@ref),
-# which, in addition to the structure factor data itself, takes a list of wave
-# vectors and a mode parameter. The options for the mode parameter are `:trace`,
-# `:perp` and `:full` which return, respectively, the trace, the unpolarized
-# intensity, and the full set of matrix elements (correlations of spin
-# components) at the specified wave vectors. An optional keyword argument `kT`
-# enables a classical-to-quantum rescaling.
-# 
-# Below, we plot two single-$q$ slices.
+# The basic functions for accessing intensity data are [`intensities_interpolated`](@ref)
+# and [`intensities_binned`](@ref). Both functions accept an [`intensity_formula`](@ref)
+# to specify how to combine the correlations recorded in the `StructureFactor` into
+# intensity data. By default, a formula computing the unpolarized intensity is used,
+# but alternative formulas can be specified.
+#
+# By way of example, we will use a formula which computes the trace of the structure
+# factor and applies a classical-to-quantum temperature-dependent rescaling `kT`.
+
+formula = intensity_formula(sf, :trace; kT = kT);
+
+# Using the formula, we plot single-$q$ slices at (0,0,0) and (π,π,π):
 
 qs = [[0, 0, 0], [0.5, 0.5, 0.5]]
-is = intensities(sf, qs, :trace; kT)
+is = intensities_interpolated(sf, qs; interpolation = :round, formula = formula)
 
 fig = lines(ωs(sf), is[1,:]; axis=(xlabel="meV", ylabel="Intensity"), label="(0,0,0)")
 lines!(ωs(sf), is[2,:]; label="(π,π,π)")
@@ -398,14 +401,14 @@ density = 40
 path, markers = connected_path(sf, points, density);
 
 # Calculate and plot the intensities along this path using [`FormFactor`](@ref)
-# corrections appropriate for `Fe2` magnetic ions.
+# corrections appropriate for `Fe2` magnetic ions, and apply a dipole correction `:perp`.
 
 formfactors = [FormFactor(1, "Fe2"; g_lande=3/2)]
+new_formula = intensity_formula(sf, :perp; kT = kT, formfactors = formfactors);
 
-is = intensities(sf, path, :perp; 
+is = intensities_interpolated(sf, path;
     interpolation = :linear,       # Interpolate between available wave vectors
-    kT,                            # Temperature for intensity correction
-    formfactors,                   # Form factor information 
+    formula = new_formula
 )
 is = broaden_energy(sf, is, (ω, ω₀)->lorentzian(ω-ω₀, 0.05))  # Add artificial broadening
 
@@ -421,7 +424,10 @@ heatmap(1:size(is,1), ωs(sf), is;
     )
 )
 
-# For comparison, we can make the same plot using histogram bins:
+# Whereas [`intensities_interpolated`](@ref) either rounds or linearly interpolates
+# between the discrete ``(Q,ω)`` points Sunny calculates correlations at, [`intensities_binned`](@ref)
+# performs histogram binning analgous to what is done in experiments.
+# The graphs produced by each method are similar.
 cut_width = 0.3
 density = 15
 paramsList, markers, ranges = connected_path_bins(sf,points,density,cut_width)
@@ -431,7 +437,7 @@ energy_bins = paramsList[1].numbins[4]
 is = zeros(Float64,total_bins,energy_bins)
 integrated_kernel = integrated_lorentzian(0.05) # Lorentzian broadening
 for k in 1:length(paramsList)
-    h,c = intensities_binned(sf,paramsList[k],:perp,kT=kT,formfactors=formfactors,integrated_kernel = integrated_kernel)
+    h,c = intensities_binned(sf,paramsList[k];formula = new_formula,integrated_kernel = integrated_kernel)
     is[ranges[k],:] = h[:,1,1,:] ./ c[:,1,1,:]
 end
 
@@ -453,11 +459,7 @@ npoints = 60
 qvals = range(-2, 2, length=npoints)
 qs = [[a, b, 0] for a in qvals, b in qvals]
 
-is = intensities(sf, qs, :perp;
-    interpolation = :linear,
-    kT,
-    formfactors,
-);
+is = intensities_interpolated(sf, qs; formula = new_formula,interpolation = :linear);
 
 ωidx = 30
 hm = heatmap(is[:,:,ωidx]; axis=(title="ω=$(ωs(sf)[ωidx]) meV", aspect=true))
@@ -480,25 +482,17 @@ A = [0.5  1  0;
      0    0  1]
 ks = [A*q for q in qs]
 
-is_ortho = intensities(sf, ks, :perp;
-    interpolation = :linear,
-    kT,
-    formfactors,
-)
+is_ortho = intensities_interpolated(sf, ks; formula = new_formula, interpolation = :linear)
 
 hm = heatmap(is_ortho[:,:,ωidx]; axis=(title="ω=$(ωs(sf)[ωidx]) meV", aspect=true))
 Colorbar(hm.figure[1,2], hm.plot)
 hidedecorations!(hm.axis); hidespines!(hm.axis)
 hm
 
 # Finally, we note that instantaneous structure factor data, ``𝒮(𝐪)``, can be
-# obtained from a dynamic structure factor with [`instant_intensities`](@ref).
+# obtained from a dynamic structure factor with [`instant_intensities_interpolated`](@ref).
 
-is_static = instant_intensities(sf, ks, :perp;
-    interpolation = :linear,
-    kT,
-    formfactors,
-)
+is_static = instant_intensities_interpolated(sf, ks; formula = new_formula, interpolation = :linear)
 
 hm = heatmap(is_static; axis=(title="Instantaneous Structure Factor", aspect=true))
 Colorbar(hm.figure[1,2], hm.plot)

examples/powder_averaging.jl

@@ -37,7 +37,7 @@ sf = DynamicStructureFactor(sys;
     nω=100,
     ωmax=5.5,
     process_trajectory=:symmetrize
-);
+)
 
 # To get some intuition about the expected results, we first look at the "single
 # crystal" results along a high-symmetry path in the first Brillouin zone. While
@@ -51,7 +51,7 @@ sf = DynamicStructureFactor(sys;
 qpoints = [[0.0, 0.0, 0.0], [0.5, 0.0, 0.0], [0.5, 0.5, 0.0], [0.0, 0.0, 0.0]]
 qs, markers = connected_path(sf, qpoints, 50) 
 
-is = intensities(sf, qs, :trace; interpolation=:none)
+is = intensities_interpolated(sf, qs; interpolation=:round, formula = intensity_formula(sf,:trace))
 is_broad = broaden_energy(sf, is, (ω, ω₀) -> lorentzian(ω-ω₀, 0.1))
 
 ## Plot results
@@ -75,14 +75,13 @@ fig
 # factor, a list of radius values (Å⁻¹), and a density parameter (Å⁻²) that will
 # control the number of wave vectors to sample at each radius. For each radius
 # `r`, the function will generate wavevectors on a sphere of this radius and
-# retrieve their [`intensities`](@ref). These intensities will be broadened, as
+# retrieve their [`intensities_interpolated`](@ref). These intensities will be broadened, as
 # just demonstrated above, and then averaged to produce a single vector of
 # energy-intensities for each `r`. Note that our `powder_average` function
-# passes most of its keywords through to [`intensities`](@ref), so it can be
-# given `kT`, `formfactors`, etc., and these parameters will be applied to the
-# calculation.
+# passes most of its keywords through to [`intensities_interpolated`](@ref), so it can be
+# given an [`intensity_formula`](@ref).
 
-function powder_average(sf, rs, density; η=0.1, mode=:perp, kwargs...)
+function powder_average(sf, rs, density; η=0.1, kwargs...)
     nω = length(ωs(sf))
     output = zeros(Float64, length(rs), nω)
 
@@ -91,7 +90,7 @@ function powder_average(sf, rs, density; η=0.1, mode=:perp, kwargs...)
         if length(qs) == 0                    
             qs = [[0., 0., 0.]]  # If no points (r is too small), just look at 0 vector
         end
-        vals = intensities(sf, qs, mode; kwargs...)  # Retrieve energy intensities
+        vals = intensities_interpolated(sf, qs; kwargs...)  # Retrieve energy intensities
         vals[:,1] .*= 0.0  # Remove elastic peaks before broadening
         vals = broaden_energy(sf, vals, (ω,ω₀)->lorentzian(ω-ω₀, η))  # Apply Lorentzian broadening
         output[i,:] = reshape(mean(vals, dims=1), (nω,))  # Average single radius results and save
@@ -106,7 +105,8 @@ rs = range(0, 6π, length=55)  # Set of radius values
 η = 0.05                      # Lorentzian broadening parameter
 density = 0.15                # Number of samples in Å⁻²
 
-pa = powder_average(sf, rs, density; η, kT);
+formula = intensity_formula(sf,:perp)
+pa = powder_average(sf, rs, density; η, formula);
 
 # and plot the results.
 
@@ -137,7 +137,7 @@ fig1
 radial_binning_parameters = (0,6π,6π/55)
 integrated_kernel = integrated_lorentzian(0.05) # Lorentzian broadening
 
-pa_intensities, pa_counts = powder_averaged_bins(sf,radial_binning_parameters,:perp,kT=kT,integrated_kernel = integrated_kernel)
+pa_intensities, pa_counts = powder_averaged_bins(sf,radial_binning_parameters;integrated_kernel = integrated_kernel,formula)
 
 pa_normalized_intensities = pa_intensities ./ pa_counts
 
@@ -147,10 +147,13 @@ rs_bincenters = axes_bincenters(radial_binning_parameters...)
 heatmap!(ax, rs_bincenters[1], ωs(sf), pa_normalized_intensities; colorrange=(0,3.0))
 fig
 
+# The striations in the graph tell us that the simulation is under resolved for this bin size.
+# We should increase the size of either the periodic lattice, or the bins.
+#
 # Using the `bzsize` option, we can even resolve the contribution from each brillouin zone:
-intensity_firstBZ, counts_firstBZ = powder_averaged_bins(sf,radial_binning_parameters,:perp,kT=kT,integrated_kernel = integrated_kernel, bzsize=(1,1,1))
+intensity_firstBZ, counts_firstBZ = powder_averaged_bins(sf,radial_binning_parameters;integrated_kernel = integrated_kernel, bzsize=(1,1,1),formula)
 #md #intensity_secondBZ, counts_secondBZ = powder_averaged_bins(..., bzsize=(2,2,2))
-intensity_secondBZ, counts_secondBZ = powder_averaged_bins(sf,radial_binning_parameters,:perp,kT=kT,integrated_kernel = integrated_kernel, bzsize=(2,2,2))#hide
+intensity_secondBZ, counts_secondBZ = powder_averaged_bins(sf,radial_binning_parameters;integrated_kernel = integrated_kernel, bzsize=(2,2,2),formula)#hide
 #md #intensity_thirdBZ, counts_thirdBZ = powder_averaged_bins(..., bzsize=(3,3,3))
 intensity_thirdBZ = pa_intensities;#hide
 counts_thirdBZ = pa_counts;#hide