Sampled correlations (#104)

`StructureFactor` renamed `SampledCorrelations. `DynamicStructureFactor` -> `dynamic_correlations` and `InstantStructureFactor` -> `instant_correlations`.

docs/src/structure-factor.md

@@ -49,21 +49,25 @@ extract information from the results. These tools are briefly outlined below.
 Please see the Examples for a "real life" use case. Detailed function
 information is available in the Library API.
 
-## Basic Usage
-
-The basic data type for calculating, storing and retrieving structure factor
-data is `StructureFactor`. Rather than creating a `StructureFactor` directly,
-one should call either [`DynamicStructureFactor`](@ref), for $𝒮^{αβ}(𝐪,ω)$, or
-[`InstantStructureFactor`](@ref), for $𝒮^{αβ}(𝐪)$. These functions will
-configure and return an appropriate `StructureFactor`.
-
-### Calculating a dynamical stucture factor: ``𝒮(𝐪,ω)``
-
-Calling `DynamicStructureFactor(sys; Δt, ωmax, nω)` will create a
-`StructureFactor` for the user and calculate an initial sample. There are three
-keywords that must be specified. These keywords will determine the dynamics
-used to calculate the sample and, consequently, the $ω$ information that will be
-available after the calculation has completed.
+## Estimating stucture factors with classical dynamics
+
+Classical dynamics may be used to estimate structure factor data by analyzing
+the spin-spin correlations of dynamical trajectories. This is fundamentally a
+Monte Carlo approach, as the trajectories must be started from an initial spin
+configuration that is sampled at thermal equilibrium. (Note that it is not
+possible to estimate a true T=0 dynamical structure factor using this method,
+but the temperature may be very low.) Samples are accumulated into a
+`SampledCorrelations`, from which intensity information may be extracted. The
+user does not typically build their own `SampledCorrelations` but instead
+initializes one by calling either `dynamical_correlations` or
+`instant_correlations`, as described below.
+
+### Estimating a dynamical structure factor: ``𝒮(𝐪,ω)``
+
+A `SampledCorrelations` for estimating the dynamical structure factor,
+$𝒮^{αβ}(𝐪,ω)$, may be created by calling [`dynamical_correlations`](@ref). This
+requires three keyword arguments. These will determine the dynamics used to
+calculate samples and, consequently, the $ω$ information that will be available. 
 
 1. `Δt`: Determines the step size used for simulating the dynamics. A smaller
    number will require proportionally more calculation time. While a smaller
@@ -81,61 +85,62 @@ available after the calculation has completed.
 3. `nω`: Determines the number of energy bins to resolve. A larger number will
    require more calculation time.
 
-Upon constructing `DynamicStructureFactor`, classical spin dynamics will be
-performed, and spin correlation data will be accumulated into $𝒮^{αβ}(𝐪,ω)$.
-The input `sys` must be a spin configuration in good thermal equilibrium, e.g.,
-using the continuous [`Langevin`](@ref) dynamics or using single spin flip
-trials with [`LocalSampler`](@ref). The statistical quality of the
-$𝒮^{αβ}(𝐪,ω)$ can be improved by generating a decorrelated spin configuration
-in `sys`, and then calling [`add_sample!`](@ref).
+A sample may be added by calling `add_sample!(sc, sys)`. The input `sys` must be
+a spin configuration in good thermal equilibrium, e.g., using the continuous
+[`Langevin`](@ref) dynamics or using single spin flip trials with
+[`LocalSampler`](@ref). The statistical quality of the $𝒮^{αβ}(𝐪,ω)$ can be
+improved by repeatedly generating decorrelated spin configurations in `sys` and
+calling `add_sample!` on each configuration.
 
 The outline of typical use case might look like this:
 ```
-# Make a `StructureFactor` and calculate an initial sample
-sf = DynamicStructureFactor(sys; Δt=0.05, ωmax=10.0, nω=100) 
+# Make a `SampledCorrelations`
+sc = dynamical_correlations(sys; Δt=0.05, ωmax=10.0, nω=100) 
 
-# Add additional samples
+# Add samples
 for _ in 1:nsamples
    decorrelate_system(sys) # Perform some type of Monte Carlo simulation
-   add_sample!(sf, sys)    # Use spins to calculate and accumulate new sample of 𝒮(𝐪,ω)
+   add_sample!(sc, sys)    # Use spins to calculate trajectory and accumulate new sample of 𝒮(𝐪,ω)
 end
 ```
-
 The calculation may be configured in a number of ways; see the
-[`DynamicStructureFactor`](@ref) documentation for a list of all keywords.
+[`dynamical_correlations`](@ref) documentation for a list of all keywords.
 
 
-### Calculating an instantaneous ("static") structure factor: ``𝒮(𝐪)``
+### Estimating an instantaneous ("static") structure factor: ``𝒮(𝐪)``
 
 Sunny provides two methods for calculating instantaneous, or static, structure
-factors: $𝒮^{αβ}(𝐪)$. The first involves calculating spin-spin correlations at
-single instances of time. The second involves calculating a dynamic structure
-factor first and integrating out the $ω$ information. The advantage of the
-latter approach is that it enables application of an $ω$-dependent
+factors: $𝒮^{αβ}(𝐪)$. The first involves calculating spatial spin-spin
+correlations at single time slices. The second involves calculating a dynamic
+structure factor first and integrating out the $ω$ information. The advantage of
+the latter approach is that it enables application of an $ω$-dependent
 classical-to-quantum rescaling of structure factor intensities, a method that
 should be preferred whenever comparing results to experimental data or spin wave
 calculations. A disadvantage of this approach is that it is computationally more
 expensive. There are also many cases when it is not straightforward to calculate
 a meaningful dynamics, as when working with Ising spins. In this section we will
-discuss how to calculate instantaneous structure factors from fixed spin
-configurations. Information about calculating instantaneous data from a dynamic
-structure factor can be found in the following section.
+discuss how to calculate instantaneous structure factors from static spin
+configurations. Information about calculating instantaneous data from a
+dynamical correlations can be found in the following section.
 
 The basic usage for the instantaneous case is very similar to the dynamic case,
-except one calls `InstantStructureFactor` instead of `DynamicStructureFactor`.
-Note that there are no required keywords as there is no need to specify any
-dynamics. `InstantStructureFactor` will immediately calculate a sample of
-$𝒮(𝐪)$ using the spin configuration contained in `sys`. It is therefore
-important that `sys` be properly thermalized before calling this function.
-Additional samples may be added with `add_sample!(sf, sys)`, just as was done in
-the dynamic case. As was true there, it is important to ensure that the spins in
-`sys` represents a new equilibrium sample before calling `add_sample!`.
-
-### Extracting information from structure factors
-
-The basic function for extracting information from a dynamic `StructureFactor`
+except one calls [`instant_correlations`](@ref) instead of
+`dynamical_correlations` to configure a `SampledCorrelations`. Note that there
+are no required keywords as there is no need to specify any dynamics.
+`instant_correlations` will return a `SampledCorrelations` containing no data.
+Samples may be added by calling `add_sample!(sc, sys)`, where `sc` is the
+`SampledCorrelations`. When performing a finite-temperature calculation, it is
+important to ensure that the spin configuration in the `sys` represents a good
+equilibrium sample, as in the dynamical case. Note, however, that we recommend
+calculating instantaneous correlations at finite temperature calculations by
+using full dynamics (i.e., using `dynamical_correlations`) and then integrating
+out the energy axis. An approach to doing this is described in the next section.
+
+### Extracting information from correlation data 
+
+The basic function for extracting information from a `SampledCorrelations`
 at a particular wave vector, $𝐪$, is [`intensities_interpolated`](@ref). It takes a
-`StructureFactor`, a list of wave vectors, and a contraction mode. For example,
+`SampledCorrelations` and a list of wave vectors. For example,
 `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.
@@ -151,9 +156,9 @@ 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_interpolated` is given the keyword argument `interpolation=:linear`, Sunny will
-use trilinear interpolation to determine a result at the requested wave
-vector. 
+`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
@@ -179,7 +184,8 @@ rescaling of intensities in the formula.
 
 To retrieve intensity data from a instantaneous structure factor, use
 [`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. 
+`intensities_interpolated`. This function may also be used to calculate
+instantaneous information from a dynamical structure factor, i.e. from a
+`SampledCorrelations` created with `dynamical_correlations`. 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. 

docs/src/versions.md

@@ -4,11 +4,20 @@ This version includes many **breaking changes**.
 
 Added support for Dipole-mode Linear Spin Wave Theory. (Thanks Hao Zhang!)
 
-Split `intensities` into calculation ([`intensity_formula`](@ref)) and presentation ([`intensities_interpolated`](@ref), [`intensities_binned`](@ref)). This is a **breaking change**, see the docs to migrate your code.
+Split `intensities` into calculation ([`intensity_formula`](@ref)) and
+presentation ([`intensities_interpolated`](@ref), [`intensities_binned`](@ref)).
+This is a **breaking change**, see the docs to migrate your code.
 
-Broadened support for custom observables in `StructureFactor` for use in `intensity_formula`.
+`StructureFactor` type renamed [`SampledCorrelations`](@ref). An appropriate
+`SampledCorrelations` is created by calling either
+[`dynamical_correlations`](@ref) or `instant_correlations`(@ref) instead of
+`DynamicStructureFactor` or `InstantStructureFactor`.
 
-Added function [`load_nxs`](@ref) to load experimental neutron scattering data to compare with `intensities_binned`.
+Broadened support for custom observables in `SampledCorrelations` for use in
+`intensity_formula`.
+
+Added function [`load_nxs`](@ref) to load experimental neutron scattering data
+to compare with `intensities_binned`.
 
 Replace `set_anisotropy!` with a new function [`set_onsite_coupling!`](@ref)
 (and similarly [`set_onsite_coupling_at!`](@ref)). The latter expects an
@@ -46,7 +55,7 @@ can be resolved by passing an explicit `offset`.
 The function [`remove_periodicity!`](@ref) disables periodicity along specified
 dimensions.
 
-Rename `StaticStructureFactor` to [`InstantStructureFactor`](@ref).
+Rename `StaticStructureFactor` to [`InstantStructureFactor`].
 
 
 # Version 0.4.2

docs/src/writevtk.md

@@ -14,7 +14,7 @@ ParaView supports volumetric rendering:
 First, generate some correlation data in Sunny.
 We will use a 2D lattice, since the correlation data ``S(Q_x,Q_y,\omega)`` is 
 3D and can be exported in its entirety.
-The following code sets up the system, thermalizes it, and records the correlation data in a [`DynamicalStructureFactor`](@ref) called `dsf`.
+The following code sets up the system, thermalizes it, and records the correlation data in a `SampledCorrelations` called `dsf`.
 
 ```julia
 using Sunny
@@ -44,11 +44,11 @@ end
 
 ωmax=10.
 
-dsf = DynamicStructureFactor(sys
-                            ;Δt=Δt
-                            ,nω=48
-                            ,ωmax=ωmax
-                            ,process_trajectory=:symmetrize)
+dsf = dynamical_correlations(sys
+                             ;Δt=Δt
+                             ,nω=48
+                             ,ωmax=ωmax
+                             ,process_trajectory=:symmetrize)
 
 nsamples = 10
 for _ in 1:nsamples

examples/binning_tutorial.jl

@@ -57,23 +57,24 @@ 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)
+sc = dynamical_correlations(sys; Δt=Δt, nω=nω, ωmax=ωmax, process_trajectory=:symmetrize)
+add_sample!(sc, sys)
 
 # 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)
+    add_sample!(sc, 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:
 
-params = unit_resolution_binning_parameters(dsf)#hide
-bin_rlu_as_absolute_units!(params,dsf)#hide
+params = unit_resolution_binning_parameters(sc)#hide
+bin_rlu_as_absolute_units!(params,sc)#hide
 lower_aabb_q, upper_aabb_q = Sunny.binning_parameters_aabb(params)#hide
 lower_aabb_cell = floor.(Int64,lower_aabb_q .* latsize .+ 1)#hide
 upper_aabb_cell = ceil.(Int64,upper_aabb_q .* latsize .+ 1)#hide
@@ -93,7 +94,7 @@ 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)
+params = unit_resolution_binning_parameters(sc)
 
 # Since this is a 4D histogram,
 # it further has to be integrated over two of those directions in order to be displayed.
@@ -104,13 +105,13 @@ integrate_axes!(params;axes = [2,4]) # Integrate over Qy (2) and E (4)
 # 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)
+formula = intensity_formula(sc, :perp; kT, formfactors)
+intensity,counts = intensities_binned(sc, 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
+intensity,counts = intensities_binned(sc, params; formula)#hide
 normalized_intensity = intensity ./ counts;#hide
 bin_centers = axes_bincenters(params);
 
@@ -127,14 +128,14 @@ 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
+params = unit_resolution_binning_parameters(sc)#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
+params = unit_resolution_binning_parameters(sc)#hide
 integrate_axes!(params;axes = [2,4])#hide
 params.binwidth[1] *= 2
 params.binwidth[3] *= 2
@@ -144,7 +145,7 @@ plot_data(params)#hide
 # over-resolved in energy:
 x = zeros(Float64,0)#hide
 y = zeros(Float64,0)#hide
-for omega = ωs(dsf), qx = unique(Qx)#hide
+for omega = ωs(sc), qx = unique(Qx)#hide
     push!(x,qx)#hide
     push!(y,omega)#hide
 end#hide
@@ -156,14 +157,14 @@ scatter!(ax,x,y)
 # See [`slice_2D_binning_parameters`](@ref).
 x_axis_bin_count = 10
 cut_width = 0.3
-params = slice_2D_binning_parameters(dsf,[0,0,0],[1,1,0],x_axis_bin_count,cut_width)
+params = slice_2D_binning_parameters(sc,[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)
+intensity,counts = intensities_binned(sc, params; formula)
 log_intensity = log10.(intensity ./ counts);
 bin_centers = axes_bincenters(params);#hide
 fig = Figure()#hide
@@ -175,7 +176,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; formula)#hide
+intensity,counts = intensities_binned(sc, params; formula)#hide
 log_intensity = log10.(intensity ./ counts);#hide
 bin_centers = axes_bincenters(params);#hide
 fig = Figure()#hide