diff --git a/Project.toml b/Project.toml index c45eaeca3..4fe05b74c 100644 --- a/Project.toml +++ b/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" diff --git a/docs/src/structure-factor.md b/docs/src/structure-factor.md index d739b3b40..049f04189 100644 --- a/docs/src/structure-factor.md +++ b/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,8 +151,8 @@ 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, @@ -160,7 +160,10 @@ 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 @@ -168,15 +171,15 @@ 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. \ No newline at end of file +calculating a static structure factor at the outset. diff --git a/examples/binning_tutorial.jl b/examples/binning_tutorial.jl index b771f0a33..e50dbbee9 100644 --- a/examples/binning_tutorial.jl +++ b/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 diff --git a/examples/fei2_tutorial.jl b/examples/fei2_tutorial.jl index 6da3cb83f..7f877b424 100644 --- a/examples/fei2_tutorial.jl +++ b/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,11 +482,7 @@ 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) @@ -492,13 +490,9 @@ 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) diff --git a/examples/powder_averaging.jl b/examples/powder_averaging.jl index 2c3f5acf2..4e1cdb946 100644 --- a/examples/powder_averaging.jl +++ b/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 diff --git a/src/SpinWaveTheory/Intensities.jl b/src/SpinWaveTheory/Intensities.jl new file mode 100644 index 000000000..3da464713 --- /dev/null +++ b/src/SpinWaveTheory/Intensities.jl @@ -0,0 +1,160 @@ +#DD this will become more similar to the existing intensities. +""" + intensities(swt::SpinWaveTheory, qs, ωvals, η::Float64) + +Computes the unpolarized inelastic neutron scattering intensities given a +`SpinWaveTheory`, an array of wave vectors `qs`, a list of energies `ωvals`, and +a Lorentzian broadening parameter `η`. + +Note that `qs` is an array of wave vectors of arbitrary dimension. Each element +``q`` of `qs` must be a 3-vector in reciprocal lattice units. I.e., ``qᵢ`` is +given in ``2π/|aᵢ|`` with ``|aᵢ|`` the lattice constant of the chemical lattice. + +The output will be an array with indices identical to `qs`. Each entry of the +array will be an unpolarized intensity. +""" +# DD: incorporate existing SF utilties (e.g., form factor) +function intensities(swt::SpinWaveTheory, qs, ωvals, η::Float64; formula = intensity_formula(swt) :: SpinWaveIntensityFormula) + (; sys) = swt + qs = Vec3.(qs) + Nm, Ns = length(sys.dipoles), sys.Ns[1] # number of magnetic atoms and dimension of Hilbert space + Nf = sys.mode == :SUN ? Ns-1 : 1 + nmodes = Nf * Nm + + num_ω = length(ωvals) + is = zeros(Float64, size(qs)..., num_ω) + if typeof(formula).parameters[1] != Float64 + error("Intensity formulae with complicated return types are not yet supported for SpinWaveTheory") + end + + for qidx in CartesianIndices(qs) + disp, intensity = formula.calc_intensity(swt,qs[qidx]) + + for band = 1:nmodes + # At a Goldstone mode, where the intensity is divergent, use a delta-function for the intensity. + if (disp[band] < 1.0e-3) && (intensity[band] > 1.0e3) + is[qidx, 1] += intensity[band] + else + for index_ω = 1:num_ω + is[qidx, index_ω] += intensity[band] * lorentzian(ωvals[index_ω]-disp[band], η) + end + end + end + end + return is +end + +struct SpinWaveIntensityFormula{T} + kT :: Float64 + formfactors + string_formula :: String + calc_intensity :: Function +end + +function Base.show(io::IO, formula::SpinWaveIntensityFormula{T}) where T + print(io,"SpinWaveIntensityFormula{$T}") +end + +function Base.show(io::IO, ::MIME"text/plain", formula::SpinWaveIntensityFormula{T}) where T + printstyled(io, "Quantum Scattering Intensity Formula\n";bold=true, color=:underline) + + formula_lines = split(formula.string_formula,'\n') + + intensity_equals = " Intensity(Q,ω) = ∑ᵢ δ(ω-ωᵢ) " + println(io,"At any Q and for each band ωᵢ = εᵢ(Q), with S = S(Q,ωᵢ):") + println(io) + println(io,intensity_equals,formula_lines[1]) + for i = 2:length(formula_lines) + precursor = repeat(' ', textwidth(intensity_equals)) + println(io,precursor,formula_lines[i]) + end + println(io) + println(io,"Intensity :: $(T) and dispersion reported for each band") +end + +intensity_formula(swt::SpinWaveTheory; kwargs...) = intensity_formula(swt, :perp; kwargs...) +function intensity_formula(swt::SpinWaveTheory, mode::Symbol; kwargs...) + if mode == :trace + contractor = Trace(swt) + string_formula = "Tr S" + elseif mode == :perp + contractor = DipoleFactor(swt) + string_formula = "∑_ij (I - Q⊗Q){i,j} S{i,j}\n\n(i,j = Sx,Sy,Sz)" + elseif mode == :full + contractor = FullTensor(swt) + string_formula = "S{α,β}" + end + intensity_formula(swt,contractor;string_formula,kwargs...) +end + +function intensity_formula(swt::SpinWaveTheory, contractor::Contraction; kwargs...) + return_type = contraction_return_type(contractor) + intensity_formula(swt,required_correlations(contractor); return_type = return_type,kwargs...) do k,ω,correlations + intensity = contract(correlations, k, contractor) + end +end + +function intensity_formula(f::Function,swt::SpinWaveTheory,corr_ix::AbstractVector{Int64}; kT = Inf, formfactors = nothing, return_type = Float64, string_formula = "f(Q,ω,S{α,β}[ix_q,ix_ω])") + (; sys, positions_chem, s̃_mat) = swt + Nm, Ns = length(sys.dipoles), sys.Ns[1] # number of magnetic atoms and dimension of Hilbert space + Nf = sys.mode == :SUN ? Ns-1 : 1 + N = Nf + 1 + nmodes = Nf * Nm + M = sys.mode == :SUN ? 1 : (Ns-1) # scaling factor (=1) if in the fundamental representation + sqrt_M = √M + sqrt_Nm_inv = 1.0 / √Nm + + # Preallocation + Hmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) + Vmat = zeros(ComplexF64, 2*nmodes, 2*nmodes) + Avec_pref = zeros(ComplexF64, Nm) + disp = zeros(Float64, nmodes) + intensity = zeros(Float64, nmodes) + + # Calculate DSSF + formula = function(swt::SpinWaveTheory,q::Vec3) + _, qmag = chemical_to_magnetic(swt, q) + + swt_hamiltonian!(swt, qmag, Hmat) + bogoliubov!(disp, Vmat, Hmat, swt.energy_tol) + + for site = 1:Nm + # note that d is the chemical coordinates + chemical_coor = positions_chem[site] + phase = exp(-2im * π * dot(q, chemical_coor)) + Avec_pref[site] = sqrt_Nm_inv * phase * sqrt_M + end + + for band = 1:nmodes + v = Vmat[:, band] + Avec = zeros(ComplexF64, 3) + for site = 1:Nm + @views tS_μ = s̃_mat[:, :, :, site] + for μ = 1:3 + for α = 2:N + Avec[μ] += Avec_pref[site] * (tS_μ[α, 1, μ] * v[(site-1)*(N-1)+α-1+nmodes] + tS_μ[1, α, μ] * v[(site-1)*(N-1)+α-1]) + end + end + end + + # DD: Generalize this based on list of arbitrary operators, optimize out symmetry, etc. + Sαβ = Matrix{ComplexF64}(undef,3,3) + Sαβ[1,1] = real(Avec[1] * conj(Avec[1])) + Sαβ[1,2] = Avec[1] * conj(Avec[2]) + Sαβ[1,3] = Avec[1] * conj(Avec[3]) + Sαβ[2,2] = real(Avec[2] * conj(Avec[2])) + Sαβ[2,3] = Avec[2] * conj(Avec[3]) + Sαβ[3,3] = real(Avec[3] * conj(Avec[3])) + Sαβ[2,1] = conj(Sαβ[1,2]) + Sαβ[3,1] = conj(Sαβ[3,1]) + Sαβ[3,2] = conj(Sαβ[2,3]) + + k = swt.recipvecs_chem * q + intensity[band] = f(k,disp[band],Sαβ[corr_ix]) + end + return disp, intensity + end + SpinWaveIntensityFormula{return_type}(kT,formfactors,string_formula,formula) +end + + diff --git a/src/SpinWaveTheory/SWTCalculations.jl b/src/SpinWaveTheory/SWTCalculations.jl index 1a205dfd0..8ed6fb8fd 100644 --- a/src/SpinWaveTheory/SWTCalculations.jl +++ b/src/SpinWaveTheory/SWTCalculations.jl @@ -470,48 +470,4 @@ function dssf(swt::SpinWaveTheory, qs) end -#DD this will become more similar to the existing intensities. -""" - intensities(swt::SpinWaveTheory, qs, ωvals, η::Float64) - -Computes the unpolarized inelastic neutron scattering intensities given a -`SpinWaveTheory`, an array of wave vectors `qs`, a list of energies `ωvals`, and -a Lorentzian broadening parameter `η`. - -Note that `qs` is an array of wave vectors of arbitrary dimension. Each element -``q`` of `qs` must be a 3-vector in reciprocal lattice units. I.e., ``qᵢ`` is -given in ``2π/|aᵢ|`` with ``|aᵢ|`` the lattice constant of the chemical lattice. - -The output will be an array with indices identical to `qs`. Each entry of the -array will be an unpolarized intensity. -""" -# DD: incorporate existing SF utilties (e.g., form factor, polarization correction) -function intensities(swt::SpinWaveTheory, qs, ωvals, η::Float64) - (; sys) = swt - qs = Vec3.(qs) - Nm, Ns = length(sys.dipoles), sys.Ns[1] # number of magnetic atoms and dimension of Hilbert space - Nf = sys.mode == :SUN ? Ns-1 : 1 - nmodes = Nf * Nm - disp, Sαβs = dssf(swt, qs) - - num_ω = length(ωvals) - is = zeros(Float64, size(qs)..., num_ω) - - for qidx in CartesianIndices(qs) - polar_mat = polarization_matrix(swt.recipvecs_chem * qs[qidx]) - - for band = 1:nmodes - band_intensity = real(sum(polar_mat .* Sαβs[qidx,band])) - # At a Goldstone mode, where the intensity is divergent, use a delta-function for the intensity. - if (disp[qidx, band] < 1.0e-3) && (band_intensity > 1.0e3) - is[qidx, 1] += band_intensity - else - for index_ω = 1:num_ω - is[qidx, index_ω] += band_intensity * lorentzian(ωvals[index_ω]-disp[qidx,band], η) - end - end - end - end - return is -end \ No newline at end of file diff --git a/src/StructureFactors/DataRetrieval.jl b/src/StructureFactors/DataRetrieval.jl index a973ebb7b..73cb5c0c8 100644 --- a/src/StructureFactors/DataRetrieval.jl +++ b/src/StructureFactors/DataRetrieval.jl @@ -2,27 +2,6 @@ # Basic functions for retrieving 𝒮(𝐪,ω) values ################################################################################ -# Internal function for getting a single 𝒮(q, ω) intensity -function calc_intensity(sf::StructureFactor, k, cell, ω, iω, contractor, kT, ffdata, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms} - elems = phase_averaged_elements(view(sf.data,:,:,:,cell,iω), k, sf, ffdata, Val(NCorr), Val(NAtoms)) - intensity = contract(elems, k, contractor) - return intensity * classical_to_quantum(ω, kT) # DDTodo: Probably make this a post-processing step -end - -classical_to_quantum(ω, kT::Float64) = ω > 0 ? ω/(kT*(1 - exp(-ω/kT))) : iszero(ω) ? 1.0 : -ω*exp(ω/kT)/(kT*(1 - exp(ω/kT))) -classical_to_quantum(ω, ::Nothing) = 1.0 - - -function prepare_form_factors(sf, formfactors) - if isnothing(formfactors) - cryst = isnothing(sf.origin_crystal) ? sf.crystal : sf.origin_crystal - class_indices = [findfirst(==(class_label), cryst.classes) for class_label in unique(cryst.classes)] - formfactors = [FormFactor{Sunny.EMPTY_FF}(; atom) for atom in class_indices] - end - formfactors = upconvert_form_factors(formfactors) # Ensure formfactors have consistent type - return propagate_form_factors(sf, formfactors) -end - """ BinningParameters(binstart,binend,binwidth;covectors = I(4)) BinningParameters(binstart,binend;numbins,covectors = I(4)) @@ -203,7 +182,7 @@ end """ unit_resolution_binning_parameters(sf::StructureFactor) -Create [`BinningParameters`](@ref) which place one histogram bin centered at each possible Sunny scattering vector. +Create [`BinningParameters`](@ref) which place one histogram bin centered at each possible `(q,ω)` scattering vector of the crystal. This is the finest possible binning without creating bins with zero scattering vectors in them. This function can be used without reference to a [`StructureFactor`](@ref) using this alternate syntax to manually specify the bin centers for the energy axis and the lattice size: @@ -337,39 +316,29 @@ connected_path_bins(sf::StructureFactor, qs::Vector, density,args...;kwargs...) connected_path_bins(sw::SpinWaveTheory, ωvals, qs::Vector, density,args...;kwargs...) = connected_path_bins(sw.recipvecs_chem, ωvals, qs, density,args...;kwargs...) - """ - intensities_binned(sf::StructureFactor, params::BinningParameters, mode) + intensity, counts = intensities_binned(sf::StructureFactor, params::BinningParameters; formula, integrated_kernel) + +Given correlation data contained in a [`StructureFactor`](@ref) and [`BinningParameters`](@ref) describing the +shape of a histogram, compute the intensity and normalization for each histogram bin using a given [`intensity_fomula`](@ref), or `intensity_formula(sf,:perp)` by default. + +This is an alternative to [`intensities_interpolated`](@ref) which bins the scattering intensities into a histogram +instead of interpolating between them at specified `qs` values. See [`unit_resolution_binning_parameters`](@ref) +for a reasonable default choice of [`BinningParameters`](@ref) which roughly emulates [`intensities_interpolated`](@ref) with `interpolation = :round`. + +If a function `integrated_kernel(Δω)` is passed, it will be used as the CDF of a kernel function for energy broadening. +For example, +`integrated_kernel = Δω -> atan(Δω/η)/pi` (c.f. [`integrated_lorentzian`](@ref) implements Lorentzian broadening with parameter `η`. +Currently, energy broadening is only supported if the [`BinningParameters`](@ref) are such that the first three axes are purely spatial and the last (energy) axis is `[0,0,0,1]`. """ -# Given correlation data contained in `sf` and BinningParameters describing the -# shape of a histogram, compute the intensity and normalization for each histogram bin. -# -# This is an alternative to `intensities` which bins the scattering intensities into a histogram -# instead of interpolating between them at specified `qs` values. See `unit_resolution_binning_parameters` -# for a reasonable default choice of `BinningParameters` which roughly emulates `intensities` with `NoInterp`. -# -# If a function `integrated_kernel(Δω)` is passed, it will be used as the CDF of a kernel function for energy broadening. -# For example, -# `integrated_kernel = Δω -> atan(Δω/η)/pi` implements `lorentzian` broadening with parameter `η`. -# Currently, energy broadening is only supported if the `BinningParameters` are such that the first three axes are purely spatial and the last (energy) axis is `[0,0,0,1]`. -function intensities_binned(sf::StructureFactor, params::BinningParameters, mode; - integrated_kernel = nothing, - kT = nothing, - formfactors = nothing, +function intensities_binned(sf::StructureFactor, params::BinningParameters; + integrated_kernel = nothing,formula = intensity_formula(sf,:perp) ) (; binwidth, binstart, binend, covectors, numbins) = params output_intensities = zeros(Float64,numbins...) output_counts = zeros(Float64,numbins...) ωvals = ωs(sf) recip_vecs = 2π*inv(sf.crystal.latvecs)' - ffdata = prepare_form_factors(sf, formfactors) - contractor = if mode == :perp - DipoleFactor(sf) - elseif typeof(mode) <: Tuple{Int, Int} - Element(sf, mode) - else - Trace(sf) - end # Find an axis-aligned bounding box containing the histogram lower_aabb_q, upper_aabb_q = binning_parameters_aabb(params) @@ -399,8 +368,7 @@ function intensities_binned(sf::StructureFactor, params::BinningParameters, mode if all(xyztBin .<= numbins) && all(xyztBin .>= 1) ci = CartesianIndex(xyztBin.data) k = recip_vecs * q - NCorr, NAtoms = size(sf.data)[1:2] - intensity = calc_intensity(sf,k,base_cell,ω,iω,contractor, kT, ffdata, Val(NCorr), Val(NAtoms)) + intensity = formula.calc_intensity(sf,k,base_cell,iω) output_intensities[ci] += intensity output_counts[ci] += 1 end @@ -416,8 +384,7 @@ function intensities_binned(sf::StructureFactor, params::BinningParameters, mode # Calculate source scattering vector intensity only once ci = CartesianIndex(xyztBin.data) k = recip_vecs * q - NCorr, NAtoms = size(sf.data)[1:2] - intensity = calc_intensity(sf,k,base_cell,ω,iω,contractor, kT, ffdata, Val(NCorr), Val(NAtoms)) + intensity = formula.calc_intensity(sf,k,base_cell,iω) # Broaden from the source scattering vector (k,ω) to # each target bin ci_other for iωother = 1:numbins[4] @@ -544,10 +511,172 @@ function pruned_stencil_info(sf::StructureFactor, qs, interp::InterpolationSchem return (; qs_all, ks_all, idcs_all, counts) end +abstract type IntensityFormula end + +struct ClassicalIntensityFormula{T} <: IntensityFormula + kT :: Float64 + formfactors + string_formula :: String + calc_intensity :: Function +end + +function Base.show(io::IO, formula::ClassicalIntensityFormula{T}) where T + print(io,"ClassicalIntensityFormula{$T}") +end + +function Base.show(io::IO, ::MIME"text/plain", formula::ClassicalIntensityFormula{T}) where T + printstyled(io, "Classical Scattering Intensity Formula\n";bold=true, color=:underline) + + formula_lines = split(formula.string_formula,'\n') + + intensity_equals = " Intensity[ix_q,ix_ω] = " + println(io,"At discrete scattering modes S = S[ix_q,ix_ω], use:") + println(io) + println(io,intensity_equals,formula_lines[1]) + for i = 2:length(formula_lines) + precursor = repeat(' ', textwidth(intensity_equals)) + println(io,precursor,formula_lines[i]) + end + println(io) + + if isnothing(formula.formfactors) + printstyled(io, "No form factors specified\n";color=:yellow) + else + printstyled(io, "Form factors included in S ✓\n";color=:green) + end + if formula.kT == Inf + printstyled(io, "No temperature correction";color=:yellow) + print(io, " (kT = ∞)\n") + else + printstyled(io, "Temperature corrected (kT = $(formula.kT)) ✓\n";color = :green) + end + if T != Float64 + println(io,"Intensity :: $(T)") + end +end + +""" + formula = intensity_formula(sf::StructureFactor; kwargs...) + formula.calc_intensity(sf,q,ix_q,ix_ω) + +Establish a formula for computing the intensity of the discrete scattering modes `(q,ω)` using the correlation data ``𝒮^{αβ}(q,ω)`` stored in the [`StructureFactor`](@ref). +The `formula` returned from `intensity_formula` can be passed to [`intensities_interpolated`](@ref) or [`intensities_binned`](@ref). + +Sunny has several built-in formulas that can be selected by setting `contraction_mode` to one of these values: + +- `:perp` (default), which contracts ``𝒮^{αβ}(q,ω)`` with the dipole factor ``δ_{αβ} - q_{α}q_{β}``, returning the unpolarized intensity. +- `:trace`, which yields ``\\operatorname{tr} 𝒮(q,ω) = ∑_α 𝒮^{αα}(q,ω)`` +- `:full`, which will return all elements ``𝒮^{αβ}(𝐪,ω)`` without contraction. + +Additionally, there are keyword arguments providing temperature and form factor corrections: + +- `kT`: If a temperature is provided, the intensities will be rescaled by a + temperature- and ω-dependent classical-to-quantum factor. `kT` should be + specified when making comparisons with spin wave calculations or + experimental data. If `kT` is not specified, infinite temperature (no correction) is assumed. +- `formfactors`: To apply form factor corrections, provide this keyword with a + vector of `FormFactor`s, one for each unique site in the unit cell. The form factors + will be symmetry propagated to all equivalent sites. + +Alternatively, a custom formula can be specifed by providing a function `intensity = f(q,ω,correlations)` and specifying which correlations it requires: + + intensity_formula(f,sf::StructureFactor, required_correlations; kwargs...) + +The function is intended to be specified using `do` notation. For example, this custom formula sums the off-diagonal correlations: + + required = [(:Sx,:Sy),(:Sy,:Sz),(:Sx,:Sz)] + intensity_formula(sf,required,return_type = ComplexF64) do k, ω, off_diagonal_correlations + sum(off_diagonal_correlations) + end + +If your custom formula returns a type other than `Float64`, use the `return_type` keyword argument to flag this. +""" +function intensity_formula(sf::StructureFactor, elem::Tuple{Symbol,Symbol}; kwargs...) + string_formula = "S{$(elem[1]),$(elem[2])}[ix_q,ix_ω]" + intensity_formula(sf,Element(sf, elem); string_formula, kwargs...) +end +#intensity_formula(sf::StructureFactor, elem::Vector{Tuple{Symbol,Symbol}}; kwargs...) = intensity_formula(sf,Element(sf, elem); kwargs...) +intensity_formula(sf::StructureFactor; kwargs...) = intensity_formula(sf, :perp; kwargs...) +function intensity_formula(sf::StructureFactor, mode::Symbol; kwargs...) + if mode == :trace + contractor = Trace(sf) + string_formula = "Tr S" + elseif mode == :perp + contractor = DipoleFactor(sf) + string_formula = "∑_ij (I - Q⊗Q){i,j} S{i,j}\n\n(i,j = Sx,Sy,Sz)" + elseif mode == :full + contractor = FullTensor(sf) + string_formula = "S{α,β}" + end + intensity_formula(sf,contractor;string_formula,kwargs...) +end + +function intensity_formula(sf::StructureFactor, contractor::Contraction; kwargs...) + return_type = contraction_return_type(contractor) + intensity_formula(sf,required_correlations(contractor); return_type = return_type,kwargs...) do k,ω,correlations + intensity = contract(correlations, k, contractor) + end +end + +function intensity_formula(f::Function,sf::StructureFactor,required_correlations; kwargs...) + # SQTODO: This corr_ix may contain repeated correlations if the user does a silly + # thing like [(:Sx,:Sy),(:Sy,:Sx)], and this can technically be optimized so it's + # not computed twice + corr_ix = lookup_correlations(sf,required_correlations) + intensity_formula(f,sf,corr_ix;kwargs...) +end + +function intensity_formula(f::Function,sf::StructureFactor,corr_ix::AbstractVector{Int64}; kT = Inf, formfactors = nothing, return_type = Float64, string_formula = "f(Q,ω,S{α,β}[ix_q,ix_ω])") + # If temperature given, ensure it's greater than 0.0 + if iszero(kT) + error("`kT` must be greater than zero.") + end + + ffdata = prepare_form_factors(sf, formfactors) + NAtoms = size(sf.data)[2] + NCorr = length(corr_ix) + + ωs_sf = ωs(sf;negative_energies=true) + formula = function (sf::StructureFactor,k::Vec3,ix_q::CartesianIndex{3},ix_ω::Int64) + correlations = phase_averaged_elements(view(sf.data,corr_ix,:,:,ix_q,ix_ω), k, sf, ffdata, Val(NCorr), Val(NAtoms)) + + ω = ωs_sf[ix_ω] + intensity = f(k,ω,correlations) * classical_to_quantum(ω, kT) + + # Having this line saves the return_type in the function closure + # so that it can be read by intensities later + intensity :: return_type + end + ClassicalIntensityFormula{return_type}(kT,formfactors,string_formula,formula) +end + +function classical_to_quantum(ω, kT::Float64) + if kT == Inf + return 1.0 + end + if ω > 0 + ω/(kT*(1 - exp(-ω/kT))) + elseif iszero(ω) + 1.0 + else + -ω*exp(ω/kT)/(kT*(1 - exp(ω/kT))) + end +end + +function prepare_form_factors(sf, formfactors) + if isnothing(formfactors) + cryst = isnothing(sf.origin_crystal) ? sf.crystal : sf.origin_crystal + class_indices = [findfirst(==(class_label), cryst.classes) for class_label in unique(cryst.classes)] + formfactors = [FormFactor{Sunny.EMPTY_FF}(; atom) for atom in class_indices] + end + formfactors = upconvert_form_factors(formfactors) # Ensure formfactors have consistent type + return propagate_form_factors(sf, formfactors) +end + + """ - intensities(sf::StructureFactor, qs, mode; interpolation = nothing, - kT = nothing, formfactors = nothing, negative_energies = false) + intensities_interpolated(sf::StructureFactor, qs; interpolation = nothing, formula = intensity_formula(sf,:perp), negative_energies = false) The basic function for retrieving ``𝒮(𝐪,ω)`` information from a `StructureFactor`. Maps an array of wave vectors `qs` to an array of structure @@ -556,35 +685,20 @@ associated with the energy index can be retrieved by calling [`ωs`](@ref). The three coordinates of each wave vector are measured in reciprocal lattice units, i.e., multiples of the reciprocal lattice vectors. -- `mode`: Should be one of `:trace`, `:perp`, or `:full`. Determines an optional - contraction on the indices ``α`` and ``β`` of ``𝒮^{αβ}(q,ω)``. Setting - `trace` yields ``∑_α 𝒮^{αα}(q,ω)``. Setting `perp` will contract - ``𝒮^{αβ}(q,ω)`` with the dipole factor ``δ_{αβ} - q_{α}q_{β}``, returning - the unpolarized intensity. Setting `full` will return all elements - ``𝒮^{αβ}(𝐪,ω)`` without contraction. - `interpolation`: Since ``𝒮(𝐪, ω)`` is calculated on a finite lattice, data is only available at discrete wave vectors. By default, Sunny will round a requested `q` to the nearest available wave vector. Linear interpolation can be applied by setting `interpolation=:linear`. -- `kT`: If a temperature is provided, the intensities will be rescaled by a - temperature- and ω-dependent classical-to-quantum factor. `kT` should be - specified when making comparisons with spin wave calculations or - experimental data. -- `formfactors`: To apply form factor corrections, provide this keyword with a - vector of `FormFactor`s, one for each unique site in the unit cell. Sunny - will symmetry propagate the results to all equivalent sites. - `negative_energies`: If set to `true`, Sunny will return the periodic extension of the energy axis. Most users will not want this. """ -function intensities(sf::StructureFactor, qs, mode; - interpolation = :none, - kT = nothing, - formfactors = nothing, +function intensities_interpolated(sf::StructureFactor, qs; + formula = intensity_formula(sf,:perp) :: ClassicalIntensityFormula, + interpolation = :round, negative_energies = false, static_warn = true ) qs = Vec3.(qs) - NCorr, NAtoms = size(sf.data)[1:2] # If working on reshaped system, assume qs given as coordinates in terms of # reciprocal vectors of original crystal and convert them to qs in terms of @@ -597,57 +711,38 @@ function intensities(sf::StructureFactor, qs, mode; # Make sure it's a dynamical structure factor if static_warn && size(sf.data, 7) == 1 - error("`intensities` given a StructureFactor with no dynamical information. Call `instant_intensities` to retrieve instantaneous (static) structure factor data.") - end - - # If temperature given, ensure it's greater than 0.0 - if !isnothing(kT) && iszero(kT) - error("`kT` must be greater than zero.") + error("`intensities_interpolated` given a StructureFactor with no dynamical information. Call `instant_intensities_interpolated` to retrieve instantaneous (static) structure factor data.") end # Set up interpolation scheme - interp = if interpolation == :none + interp = if interpolation == :round NoInterp() elseif interpolation == :linear LinearInterp() end - # Set up element contraction - contractor = if mode == :trace - Trace(sf) - elseif mode == :perp - DipoleFactor(sf) - elseif mode == :full - FullTensor(sf) - elseif typeof(mode) <: Tuple{Int, Int} - Element(sf, mode) - end - - # Propagate form factor information (if any) - ffdata = prepare_form_factors(sf, formfactors) - # Precompute index information and preallocate ωvals = ωs(sf; negative_energies) nω = length(ωvals) stencil_info = pruned_stencil_info(sf, qs, interp) - intensities = zeros(contractor, size(qs)..., nω) + intensities = zeros(formula.calc_intensity.return_type, size(qs)..., nω) # Call type stable version of the function - intensities!(intensities, sf, qs, ωvals, interp, contractor, kT, ffdata, stencil_info, Val(NCorr), Val(NAtoms)) + intensities!(intensities, sf, qs, ωvals, interp, formula, stencil_info, formula.calc_intensity.return_type) return intensities end # Actual intensity calculation -function intensities!(intensities, sf::StructureFactor, q_targets::Array, ωvals, interp::InterpolationScheme, contraction::Contraction{T}, temp, ffdata, stencil_info, ::Val{NCorr}, ::Val{NAtoms}) where {T, NCorr, NAtoms} +function intensities!(intensities, sf::StructureFactor, q_targets::Array, ωvals, interp::InterpolationScheme{NInterp}, formula, stencil_info, T) where {NInterp} li_intensities = LinearIndices(intensities) ci_qs = CartesianIndices(q_targets) (; qs_all, ks_all, idcs_all, counts) = stencil_info for (iω, ω) in enumerate(ωvals) iq = 0 for (qs, ks, idcs, numrepeats) in zip(qs_all, ks_all, idcs_all, counts) - local_intensities = stencil_intensities(sf, ks, idcs, ω, iω, interp, contraction, temp, ffdata, Val(NCorr), Val(NAtoms)) + local_intensities = SVector{NInterp, T}(formula.calc_intensity(sf, ks[n], idcs[n], iω) for n in 1:NInterp) for _ in 1:numrepeats iq += 1 idx = li_intensities[CartesianIndex(ci_qs[iq], iω)] @@ -659,17 +754,17 @@ function intensities!(intensities, sf::StructureFactor, q_targets::Array, ωvals end """ - instant_intensities(sf::StructureFactor, qs, mode; kwargs...) + instant_intensities_interpolated(sf::StructureFactor, qs; kwargs...) Return ``𝒮(𝐪)`` intensities at wave vectors `qs`. The functionality is very -similar to [`intensities`](@ref), except the returned array has dimensions +similar to [`intensities_interpolated`](@ref), except the returned array has dimensions identical to `qs`. If called on a `StructureFactor` with dynamical information, i.e., ``𝒮(𝐪,ω)``, the ``ω`` information is integrated out. """ -function instant_intensities(sf::StructureFactor, qs, mode; kwargs...) +function instant_intensities_interpolated(sf::StructureFactor, qs; kwargs...) datadims = size(qs) ndims = length(datadims) - vals = intensities(sf, qs, mode; static_warn=false, kwargs...) + vals = intensities_interpolated(sf, qs; static_warn=false, kwargs...) static_vals = sum(vals, dims=(ndims+1,)) return reshape(static_vals, datadims) end @@ -729,7 +824,7 @@ integrated_lorentzian(η) = x -> atan(x/η)/π Performs a real-space convolution along the energy axis of an array of intensities. Assumes the format of the intensities array corresponds to what -would be returned by [`intensities`](@ref). `kernel` must be a function that +would be returned by [`intensities_interpolated`](@ref). `kernel` must be a function that takes two numbers: `kernel(ω, ω₀)`, where `ω` is a frequency, and `ω₀` is the center frequency of the kernel. Sunny provides [`lorentzian`](@ref) for the most common use case: diff --git a/src/StructureFactors/ElementContraction.jl b/src/StructureFactors/ElementContraction.jl index 9fca52d52..5890661b1 100644 --- a/src/StructureFactors/ElementContraction.jl +++ b/src/StructureFactors/ElementContraction.jl @@ -7,7 +7,9 @@ struct Trace{N} <: Contraction{Float64} indices :: SVector{N, Int64} end -struct DipoleFactor <: Contraction{Float64} end +struct DipoleFactor <: Contraction{Float64} + indices :: SVector{6,Int64} +end struct Element <: Contraction{ComplexF64} index :: Int64 @@ -19,41 +21,60 @@ struct FullTensor{NCorr} <: Contraction{SVector{NCorr, ComplexF64}} end ################################################################################ # Constructors ################################################################################ +Trace(swt::SpinWaveTheory) = Trace([1,5,9]) + function Trace(sf::StructureFactor{N}) where {N} # Collect all indices for matrix elements 𝒮^αβ where α=β indices = Int64[] - for (ci, idx) in sf.idxinfo - α, β = ci.I - if α == β - push!(indices, idx) + for (ki,i) = sf.observable_ixs + autocorrelation_index = CartesianIndex(i,i) + if haskey(sf.correlations,autocorrelation_index) + push!(indices,sf.correlations[autocorrelation_index]) + else + problematic_correlation = ki + error("Can't calculate trace because auto-correlation of the $problematic_correlation observable was not computed.") end end - # Check that there are the correct number of such elements - if N == 0 || sf.dipole_corrs - if length(indices) != 3 - error("Not all diagonal elements of the structure factor have been computed. Can't calculate trace.") - end - else - if length(indices) != N*N-1 - error("Not all diagonal elements of the structure factor have been computed. Can't calculate trace.") - end + + # SQ N.B.: This error doesn't make much sense, does it? + # So what if they used a different number from the default number of observables? + # Case in point: If you are doing dipole correlations in SU(N) mode, you're not taking + # the full trace, and this will error out. + + #= + total_autocorrelations = N == 0 ? 3 : N*N-1 + if length(indices) != total_autocorrelations + error("Unexpected number of observables were encounted. Expected $total_autocorrelations but actually have $(length(sf.observables)): $(keys(sf.observable_ixs))") end + =# + indices = sort(indices) - return Trace(SVector{length(indices), Int64}(indices)) + Trace(SVector{length(indices), Int64}(indices)) end -function DipoleFactor(sf::StructureFactor{N}) where {N} - if collect(keys(sf.idxinfo))[1:6] == [(1, 1), (1, 2), (2, 2), (1, 3), (2, 3), (3, 3)] .|> CartesianIndex - return DipoleFactor() +DipoleFactor(swt::SpinWaveTheory) = DipoleFactor([1,4,5,7,8,9]) + +function DipoleFactor(sf::StructureFactor{N}; spin_components = [:Sx,:Sy,:Sz]) where {N} + # Ensure that the observables themselves are present + for si in spin_components + if !haskey(sf.observable_ixs,si) + error("Observable $(si) missing, but required for dipole correction factor") + end end - error("Need to be in structure factor dipole mode to calculate depolarization correction.") + + # Ensure that the required correlations are also present + sx,sy,sz = spin_components + dipole_correlations = [(sx,sx),(sx,sy),(sy,sy),(sx,sz),(sy,sz),(sz,sz)] + indices = lookup_correlations(sf,dipole_correlations; err_msg = αβ -> "Missing correlation $(αβ), which is required to compute the depolarization correction.") + DipoleFactor(indices) end -function Element(sf::StructureFactor, pair) - index = sf.idxinfo[CartesianIndex(pair)] - return Element(index) +function Element(sf::StructureFactor, pair::Tuple{Symbol,Symbol}) + Element(only(lookup_correlations(sf,[pair]; err_msg = pair -> "Missing correlation $(pair), which was requested."))) end +FullTensor(swt::SpinWaveTheory) = FullTensor{9}() + function FullTensor(sf::StructureFactor{N}) where {N} FullTensor{size(sf.data, 1)}() end @@ -69,43 +90,42 @@ end ################################################################################ # Contraction methods ################################################################################ -function contract(elems, _, traceinfo::Trace) - intensity = 0.0 - for i in traceinfo.indices - # Diagonal elements should be real only. Finite imaginary component is - # usually on order 1e-17 and is due to roundoff in phase_averaged_elements. - intensity += real(elems[i]) - end - return intensity -end -function contract(elems, k::Vec3, ::DipoleFactor) + +# Diagonal elements should be real only. Finite imaginary component is +# usually on order 1e-17 and is due to roundoff in phase_averaged_elements. +contract(diagonal_elements, _, ::Trace) = sum(real(diagonal_elements)) + +function contract(dipole_elements, k::Vec3, dipoleinfo::DipoleFactor) dip_factor = polarization_matrix(k) # Note, can just take the real part since: # (1) diagonal elements are real by construction, and # (2) pairs of off diagonal contributions have the form x*conj(y) + conj(x)*y = 2real(x*conj(y)). - # Note also that if in dipole mode, which is guarenteed if the code makes it here, - # the order of indices is also guaranteed. - return dip_factor[1,1]*real(elems[1]) + - 2dip_factor[1,2]*real(elems[2]) + - dip_factor[2,2]*real(elems[4]) + # Note order - 2dip_factor[1,3]*real(elems[3]) + - 2dip_factor[2,3]*real(elems[5]) + - dip_factor[3,3]*real(elems[6]) + return dip_factor[1,1]*real(dipole_elements[1]) + + 2dip_factor[1,2]*real(dipole_elements[2]) + + dip_factor[2,2]*real(dipole_elements[3]) + + 2dip_factor[1,3]*real(dipole_elements[4]) + + 2dip_factor[2,3]*real(dipole_elements[5]) + + dip_factor[3,3]*real(dipole_elements[6]) end -function contract(elems, _, elem::Element) - return elems[elem.index] -end +contract(specific_element, _, ::Element) = only(specific_element) -function contract(elems, _, ::FullTensor) - return elems -end +contract(all_elems, _, ::FullTensor) = all_elems + +################################################################################ +# Contraction utils +################################################################################ +required_correlations(traceinfo::Trace) = traceinfo.indices +required_correlations(dipoleinfo::DipoleFactor) = dipoleinfo.indices +required_correlations(eleminfo::Element) = [eleminfo.index] +required_correlations(::FullTensor{NCorr}) where NCorr = 1:NCorr ################################################################################ # Contraction utils ################################################################################ -Base.zeros(::Contraction{T}, dims...) where T = zeros(T, dims...) \ No newline at end of file +Base.zeros(::Contraction{T}, dims...) where T = zeros(T, dims...) +contraction_return_type(::Contraction{T}) where T = T diff --git a/src/StructureFactors/Interpolation.jl b/src/StructureFactors/Interpolation.jl index 23b3163f6..e2ab4101f 100644 --- a/src/StructureFactors/Interpolation.jl +++ b/src/StructureFactors/Interpolation.jl @@ -6,10 +6,6 @@ struct LinearInterp <: InterpolationScheme{8} end ddtodo: Explanation of interpolation "API" =# -function stencil_intensities(sf::StructureFactor, ks, idcs, ω, iω, ::InterpolationScheme{NumInterp}, contraction::Contraction{T}, temp, ffdata, ::Val{NCorr}, ::Val{NAtoms}) where {NumInterp, T, NCorr, NAtoms} - return SVector{NumInterp, T}(calc_intensity(sf, ks[n], idcs[n], ω, iω, contraction, temp, ffdata, Val(NCorr), Val(NAtoms)) for n in 1:NumInterp) -end - function interpolated_intensity(::StructureFactor, _, _, stencil_intensities, ::NoInterp) return only(stencil_intensities) end @@ -72,4 +68,4 @@ function stencil_points(sf::StructureFactor, q, ::LinearInterp) map(i -> mod(m[i], Ls[i])+1, (1, 2, 3)) |> CartesianIndex{3} end return ms, ims -end \ No newline at end of file +end diff --git a/src/StructureFactors/PowderAveraging.jl b/src/StructureFactors/PowderAveraging.jl index 9df99e33a..f0154aa29 100644 --- a/src/StructureFactors/PowderAveraging.jl +++ b/src/StructureFactors/PowderAveraging.jl @@ -28,7 +28,7 @@ function spherical_shell(sf::StructureFactor, radius, density) end end -function powder_average(sf::StructureFactor, q_ias, mode, density; kwargs...) +function powder_average(sf::StructureFactor, q_ias, density; kwargs...) A = inv(inv(sf.crystal.latvecs)') # Transformation to convert from inverse angstroms to RLUs nω = length(ωs(sf)) output = zeros(Float64, length(q_ias), nω) # generalize this so matches contract @@ -38,7 +38,7 @@ function powder_average(sf::StructureFactor, q_ias, mode, density; kwargs...) numpoints = round(Int, area*density) fibpoints = numpoints == 0 ? [Vec3(0,0,0)] : r .* spherical_points_fibonacci(numpoints) qs = map(v->A*v, fibpoints) - vals = intensities(sf, qs, mode; kwargs...) + vals = intensities(sf, qs; kwargs...) vals = sum(vals, dims=1) / size(vals, 1) output[i,:] .= vals[1,:] end @@ -46,20 +46,23 @@ function powder_average(sf::StructureFactor, q_ias, mode, density; kwargs...) return output end -# Similar to `powder_average`, but the data is binned instead of interpolated. -# Also similar to `intensities_binned`, but the histogram x-axis is `|k|` in absolute units, which -# is a nonlinear function of `kx`,`ky`,`kz`. The y-axis is energy. -# -# Binning parameters are specified as tuples `(start,end,bin_width)`, -# e.g. `radial_binning_parameters = (0,6π,6π/55)`. -# -# Energy broadening is support in the same as `intensities_binned`. -function powder_averaged_bins(sf::StructureFactor, radial_binning_parameters, mode; +""" + powder_averaged_bins(sf::StructureFactor, radial_binning_parameters; formula + ω_binning_parameters, integrated_kernel = nothing, bzsize = nothing) + +Similar to [`powder_average`](@ref), but the data is binned instead of interpolated. +Also similar to [`intensities_binned`](@ref), but the histogram x-axis is `|k|` in absolute units, which is a nonlinear function of `kx`,`ky`,`kz`. The y-axis is energy. + +Binning parameters are specified as tuples `(start,end,bin_width)`, +e.g. `radial_binning_parameters = (0,6π,6π/55)`. + +Energy broadening is supported in the same way as `intensities_binned`. +""" +function powder_averaged_bins(sf::StructureFactor, radial_binning_parameters; ω_binning_parameters=unit_resolution_binning_parameters(ωs(sf)), - integrated_kernel=nothing, - bzsize=nothing, - kT=nothing, - formfactors=nothing, + integrated_kernel = nothing, + formula = intensity_formula(sf,:perp) :: ClassicalIntensityFormula, + bzsize=nothing ) ωstart,ωend,ωbinwidth = ω_binning_parameters rstart,rend,rbinwidth = radial_binning_parameters @@ -71,14 +74,6 @@ function powder_averaged_bins(sf::StructureFactor, radial_binning_parameters, mo output_counts = zeros(Float64,r_bin_count,ω_bin_count) ωvals = ωs(sf) recip_vecs = 2π*inv(sf.crystal.latvecs)' - ffdata = prepare_form_factors(sf, formfactors) - contractor = if mode == :perp - DipoleFactor(sf) - elseif typeof(mode) <: Tuple{Int, Int} - Element(sf, mode) - else - Trace(sf) - end # Loop over every scattering vector Ls = sf.latsize @@ -106,16 +101,14 @@ function powder_averaged_bins(sf::StructureFactor, radial_binning_parameters, mo # Check if the ω falls within the histogram ωbin = 1 .+ floor.(Int64,(ω .- ωstart) ./ ωbinwidth) if ωbin <= ω_bin_count && ωbin >= 1 - NCorr, NAtoms = size(sf.data)[1:2] - intensity = calc_intensity(sf,k,base_cell,ω,iω,contractor, kT, ffdata, Val(NCorr), Val(NAtoms)) + intensity = formula.calc_intensity(sf,k,base_cell,iω) output_intensities[rbin,ωbin] += intensity output_counts[rbin,ωbin] += 1 end else # `Energy broadening into bins' logic # Calculate source scattering vector intensity only once - NCorr, NAtoms = size(sf.data)[1:2] - intensity = calc_intensity(sf,k,base_cell,ω,iω,contractor, kT, ffdata, Val(NCorr), Val(NAtoms)) + intensity = formula.calc_intensity(sf,k,base_cell,iω) # Broaden from the source scattering vector (k,ω) to # each target bin (rbin,ωbin_other) for ωbin_other = 1:ω_bin_count diff --git a/src/StructureFactors/SampleGeneration.jl b/src/StructureFactors/SampleGeneration.jl index ed5448c03..340e496cd 100644 --- a/src/StructureFactors/SampleGeneration.jl +++ b/src/StructureFactors/SampleGeneration.jl @@ -1,66 +1,47 @@ -function observable_expectations!(buf, sys::System{N}, ops′::Array{ComplexF64, 3}) where N - Zs = sys.coherents - num_ops = size(ops′, 3) - ops = reinterpret(SMatrix{N, N, ComplexF64, N*N}, reshape(ops′, N*N, num_ops)) - - for (i, op) in enumerate(ops) - for site in all_sites(sys) - buf[i,site] = dot(Zs[site], op, Zs[site]) +function observable_values!(buf, sys::System{N}, ops; apply_g = true) where N + if N == 0 + for site in all_sites(sys), (i, op) in enumerate(ops) + dipole = sys.dipoles[site] + if apply_g + dipole = sys.gs[site] * dipole + end + buf[i,site] = op * dipole end - end - - return nothing -end - -function expectation_trajectory(sys, integrator, nsnaps, ops; kwargs...) - num_ops = size(ops, 3) - - traj_buf = zeros(ComplexF64, num_ops, size(sys.coherents)..., nsnaps) - expectation_trajectory!(traj_buf, sys, integrator, nsnaps, ops; kwargs...) - - return traj_buf -end - -function expectation_trajectory!(buf, sys, integrator, nsnaps, ops; measperiod = 1) - @assert size(ops, 3) == size(buf, 1) - - observable_expectations!(@view(buf[:,:,:,:,:,1]), sys, ops) - for n in 2:nsnaps - for _ in 1:measperiod - step!(sys, integrator) + else + Zs = sys.coherents + #num_ops = size(ops′, 3) + #ops = reinterpret(SMatrix{N, N, ComplexF64, N*N}, reshape(ops′, N*N, num_ops)) + + # SQTODO: This allocates :( + for (i, op) in enumerate(ops) + matrix_operator = convert(Matrix{ComplexF64},op) + for site in all_sites(sys) + buf[i,site] = dot(Zs[site], matrix_operator, Zs[site]) + end end - observable_expectations!(@view(buf[:,:,:,:,:,n]), sys, ops) end return nothing end +function trajectory(sys::System{N}, integrator, nsnaps, ops; kwargs...) where N + num_ops = length(ops) -function compute_mag!(M, sys::System, apply_g = true) - for site in all_sites(sys) - if apply_g - M[:, site] .= sys.gs[site] * sys.dipoles[site] - else - M[:, site] .= sys.dipoles[site] - end - end -end + traj_buf = zeros(N == 0 ? Float64 : ComplexF64, num_ops, sys.latsize..., natoms(sys.crystal), nsnaps) + trajectory!(traj_buf, sys, integrator, nsnaps, ops; kwargs...) -function dipole_trajectory(sys, integrator, nsnaps; kwargs...) - traj_buf = zeros(ComplexF64, 3, size(sys.dipoles)..., nsnaps) - dipole_trajectory!(traj_buf, sys, integrator, nsnaps; kwargs...) return traj_buf end -function dipole_trajectory!(buf, sys, integrator, nsnaps; measperiod = 1, apply_g = true) - @assert size(buf, 1) == 3 +function trajectory!(buf, sys, integrator, nsnaps, ops; measperiod = 1, apply_g = true) + @assert length(ops) == size(buf, 1) - compute_mag!(@view(buf[:,:,:,:,:,1]), sys, apply_g) + observable_values!(@view(buf[:,:,:,:,:,1]), sys, ops, apply_g = apply_g) for n in 2:nsnaps for _ in 1:measperiod step!(sys, integrator) end - compute_mag!(@view(buf[:,:,:,:,:,n]), sys, apply_g) + observable_values!(@view(buf[:,:,:,:,:,n]), sys, ops, apply_g = apply_g) end return nothing @@ -72,11 +53,7 @@ function new_sample!(sf::StructureFactor, sys::System; processtraj! = no_process @assert size(sys.dipoles) == size(samplebuf)[2:5] "`System` size not compatible with given `StructureFactor`" - if size(sf.observables)[1] == 0 - dipole_trajectory!(samplebuf, sys, integrator, nsnaps; measperiod, apply_g) - else - expectation_trajectory!(samplebuf, sys, integrator, nsnaps, sf.observables; measperiod) - end + trajectory!(samplebuf, sys, integrator, nsnaps, sf.observables; measperiod = measperiod, apply_g = apply_g) processtraj!(sf) @@ -103,15 +80,14 @@ function no_processing(::StructureFactor) end function accum_sample!(sf::StructureFactor) - (; data, idxinfo, samplebuf, copybuf, nsamples, fft!) = sf - latsize = size(samplebuf)[2:4] - natoms, nω = size(samplebuf)[5:6] + (; data, correlations, samplebuf, copybuf, nsamples, fft!) = sf + natoms = size(samplebuf)[5] fft! * samplebuf # Apply pre-planned and pre-normalized FFT nsamples[1] += 1 # Transfer to final memory layout while accumulating new samplebuf - for j in 1:natoms, i in 1:natoms, (ci, c) in idxinfo + for j in 1:natoms, i in 1:natoms, (ci, c) in correlations α, β = ci.I @. copybuf = @views samplebuf[α,:,:,:,i,:] * conj(samplebuf[β,:,:,:,j,:]) - data[c,i,j,:,:,:,:] @views data[c,i,j,:,:,:,:] .+= copybuf .* (1/nsamples[1]) diff --git a/src/StructureFactors/StructureFactors.jl b/src/StructureFactors/StructureFactors.jl index 5141ef553..d2e3e848f 100644 --- a/src/StructureFactors/StructureFactors.jl +++ b/src/StructureFactors/StructureFactors.jl @@ -6,9 +6,9 @@ struct StructureFactor{N} Δω :: Float64 # Energy step size # Correlation info (αβ indices of 𝒮^{αβ}(q,ω)) - dipole_corrs :: Bool # Whether using all correlations from dipoles - observables :: Array{ComplexF64, 3} # Operators corresponding to observables - idxinfo :: SortedDict{CartesianIndex{2}, Int64} # (α, β) to save from 𝒮^{αβ}(q, ω) + observables :: Vector{LinearMap} # Operators corresponding to observables + observable_ixs :: Dict{Symbol,Int64} # User-defined observable names + correlations :: SortedDict{CartesianIndex{2}, Int64} # (α, β) to save from 𝒮^{αβ}(q, ω) # Specs for sample generation and accumulation samplebuf :: Array{ComplexF64, 6} # New sample buffer @@ -21,9 +21,47 @@ struct StructureFactor{N} processtraj! :: Function # Function to perform post-processing on sample trajectories end -function Base.show(io::IO, ::MIME"text/plain", sf::StructureFactor) - modename = sf.dipole_corrs ? "Dipole correlations" : "Custom correlations" - printstyled(io, "StructureFactor [$modename]\n"; bold=true, color=:underline) +function Base.show(io::IO, sf::StructureFactor{N}) where N + modename = N == 0 ? "Dipole" : "SU($(N))" + print(io,"StructureFactor{$modename}") + print(io,all_observable_names(sf)) +end + +function Base.show(io::IO, ::MIME"text/plain", sf::StructureFactor{N}) where N + printstyled(io, "StructureFactor";bold=true, color=:underline) + modename = N == 0 ? "Dipole" : "SU($(N))" + print(io," ($(Base.format_bytes(Base.summarysize(sf))))\n") + print(io,"[") + if size(sf.data)[7] == 1 + printstyled(io,"S(q)";bold=true) + else + printstyled(io,"S(q,ω)";bold=true) + print(io," | nω = $(size(sf.data)[7])") + end + print(io," | $(sf.nsamples[1]) sample") + (sf.nsamples[1] > 1) && print(io,"s") + print(io,"]\n") + println(io,"Lattice: $(sf.latsize)×$(natoms(sf.crystal))") + print(io,"$(size(sf.data)[1]) correlations in $modename mode:\n") + + # Reverse the dictionary + observable_names = Dict(value => key for (key, value) in sf.observable_ixs) + + for i = 1:length(sf.observables) + print(io,i == 1 ? "╔ " : i == length(sf.observables) ? "╚ " : "║ ") + for j = 1:length(sf.observables) + if i > j + print(io,"⋅ ") + elseif haskey(sf.correlations,CartesianIndex(i,j)) + print(io,"⬤ ") + else + print(io,"• ") + end + end + print(io,observable_names[i]) + println(io) + end + printstyled(io,"") end Base.getproperty(sf::StructureFactor, sym::Symbol) = sym == :latsize ? size(sf.samplebuf)[2:4] : getfield(sf,sym) @@ -42,6 +80,27 @@ function merge!(sf::StructureFactor, others...) end end +# Finds the linear index according to sf.correlations of each correlation in corrs, where +# corrs is of the form [(:A,:B),(:B,:C),...] where :A,:B,:C are observable names. +function lookup_correlations(sf::StructureFactor,corrs; err_msg = αβ -> "Missing correlation $(αβ)") + indices = Vector{Int64}(undef,length(corrs)) + for (i,(α,β)) in enumerate(corrs) + αi = sf.observable_ixs[α] + βi = sf.observable_ixs[β] + # Make sure we're looking up the correlation with its properly sorted name + αi,βi = minmax(αi,βi) + idx = CartesianIndex(αi,βi) + + # Get index or fail with an error + indices[i] = get!(() -> error(err_msg(αβ)),sf.correlations,idx) + end + indices +end + +function all_observable_names(sf::StructureFactor) + observable_names = Dict(value => key for (key, value) in sf.observable_ixs) + [observable_names[i] for i in 1:length(observable_names)] +end """ StructureFactor @@ -55,35 +114,85 @@ function StructureFactor(sys::System{N}; Δt, nω, measperiod, process_trajectory = :none) where {N} # Set up correlation functions (which matrix elements αβ to save from 𝒮^{αβ}) - default_observables = false - default_correlations = false if isnothing(observables) - observables = zeros(ComplexF64, 0, 0, 3) # observables are empty in this case - default_observables = true + # Default observables are spin x,y,z + # projections (SU(N) mode) or components (dipole mode) + observable_ixs = Dict(:Sx => 1,:Sy => 2,:Sz => 3) + if N == 0 + dipole_component(i) = FunctionMap{Float64}(s -> s[i],1,3) + observables = dipole_component.([1,2,3]) + else + # SQTODO: Make this use the more optimized expected_spin function + # Doing this will also, by necessity, allow users to make the same + # type of optimization for their vector-valued observables. + observables = LinearMap{ComplexF64}.(spin_matrices(N)) + end else - (N == 0) && error("Structure Factor Error: Cannot provide matrices for observables when using dipolar `System`") + # If it was given as a list, preserve the user's preferred + # ordering of observables + if observables isa AbstractVector + # If they are pairs (:A => [...]), use the names + # and otherwise use alphabetical names + if !isempty(observables) && observables[1] isa Pair + observables = OrderedDict(observables) + else + dict = OrderedDict{Symbol,LinearMap}() + for i = 1:length(observables) + dict[Symbol('A' + i - 1)] = observables[i] + end + observables = dict + end + end + + # If observables were provided as (:name => matrix) pairs, + # reformat them to (:name => idx) and matrices[idx] + observable_ixs = Dict{Symbol,Int64}() + matrices = Vector{LinearMap}(undef,length(observables)) + for (i,name) in enumerate(keys(observables)) + next_available_ix = length(observable_ixs) + 1 + if haskey(observable_ixs,name) + error("Repeated observable name $name not allowed.") + end + observable_ixs[name] = next_available_ix + + # Convert dense matrices to LinearMap + if observables[name] isa Matrix + matrices[i] = LinearMap(observables[name]) + else + matrices[i] = observables[name] + end + end + observables = matrices end - nops = size(observables, 3) + + # By default, include all correlations if isnothing(correlations) correlations = [] - for i in 1:nops, j in i:nops - push!(correlations, (i, j)) + for oi in keys(observable_ixs), oj in keys(observable_ixs) + push!(correlations, (oi, oj)) end - default_correlations = true + elseif correlations isa AbstractVector{Tuple{Int64,Int64}} + # If the user used numeric indices to describe the correlations, + # we need to convert it to the names, so need to temporarily reverse + # the dictionary. + observable_names = Dict(value => key for (key, value) in observable_ixs) + correlations = [(observable_names[i],observable_names[j]) for (i,j) in correlations] end - dipole_corrs = default_observables && default_correlations - - # Construct look-up table for matrix elements - count = 1 - pairs = [] - for αβ in correlations - α, β = αβ - α, β = α < β ? (α, β) : (β, α) # Because SF is symmetric, only save diagonal and upper triangular - push!(pairs, (α, β) => count) - count += 1 + + # Construct look-up table for correlation matrix elements + idxinfo = SortedDict{CartesianIndex{2},Int64}() # CartesianIndex's sort to fastest order + for (α,β) in correlations + αi = observable_ixs[α] + βi = observable_ixs[β] + # Because correlation matrix is symmetric, only save diagonal and upper triangular + # by ensuring that all pairs are in sorted order + αi,βi = minmax(αi,βi) + idx = CartesianIndex(αi,βi) + + # Add this correlation to the list if it's not already listed + get!(() -> length(idxinfo) + 1,idxinfo,idx) end - pairs = map(i -> CartesianIndex(i.first) => i.second, pairs) # Convert to CartesianIndices - idxinfo = SortedDict{CartesianIndex{2}, Int64}(pairs) # CartesianIndices sort to fastest order + correlations = idxinfo # Set up trajectory processing function (e.g., symmetrize) processtraj! = if process_trajectory == :none @@ -99,7 +208,7 @@ function StructureFactor(sys::System{N}; Δt, nω, measperiod, # Preallocation na = natoms(sys.crystal) ncorr = length(correlations) - samplebuf = zeros(ComplexF64, nops, sys.latsize..., na, nω) + samplebuf = zeros(ComplexF64, length(observables), sys.latsize..., na, nω) copybuf = zeros(ComplexF64, sys.latsize..., nω) data = zeros(ComplexF64, ncorr, na, na, sys.latsize..., nω) @@ -114,8 +223,8 @@ function StructureFactor(sys::System{N}; Δt, nω, measperiod, origin_crystal = !isnothing(sys.origin) ? sys.origin.crystal : nothing # Make Structure factor and add an initial sample - sf = StructureFactor{N}(data, sys.crystal, origin_crystal, Δω, dipole_corrs, - observables, idxinfo, samplebuf, fft!, copybuf, measperiod, apply_g, integrator, + sf = StructureFactor{N}(data, sys.crystal, origin_crystal, Δω, + observables, observable_ixs, correlations, samplebuf, fft!, copybuf, measperiod, apply_g, integrator, nsamples, processtraj!) add_sample!(sf, sys; processtraj!) @@ -130,8 +239,8 @@ end Creates a `StructureFactor` for calculating and storing ``𝒮(𝐪,ω)`` data. This information will be obtained by running dynamical spin simulations on equilibrium snapshots, and measuring pair-correlations. The ``𝒮(𝐪,ω)`` data -can be retrieved by calling [`intensities`](@ref). Alternatively, -[`instant_intensities`](@ref) will integrate out ``ω`` to obtain ``𝒮(𝐪)``, +can be retrieved by calling [`intensities_interpolated`](@ref). Alternatively, +[`instant_intensities_interpolated`](@ref) will integrate out ``ω`` to obtain ``𝒮(𝐪)``, optionally applying classical-to-quantum correction factors. Prior to calling `DynamicStructureFactor`, ensure that `sys` represents a good @@ -153,18 +262,17 @@ Additional keyword options are the following: `:symmetrize`. The latter will symmetrize the trajectory in time, which can be useful for removing Fourier artifacts that arise when calculating the correlations. -- `observables`: Enables an advanced feature for SU(_N_) mode, allowing the user - to specify custom observables other than the three components of the dipole. - To use this features, `observables` must be given an `N×N×numops` array, - where the final index is used to retrieve each `N×N` operator. +- `observables`: Allows the user to specify custom observables. The `observables` + must be given as a list of complex `N×N` matrices or `LinearMap`s. It's + recommended to name each observable, for example: + `observables = [:A => a_observable_matrix, :B => b_map, ...]`. + By default, Sunny uses the 3 components of the dipole, `:Sx`, `:Sy` and `:Sz`. - `correlations`: Specify which correlation functions are calculated, i.e. which matrix elements ``αβ`` of ``𝒮^{αβ}(q,ω)`` are calculated and stored. Specified with a vector of tuples. By default Sunny records all auto- and - cross-correlations generated by the x, y, and z dipolar components (1, 2, - and 3 respectively). To retain only the xx and xy correlations, one would - set `correlations=[(1,1), (1,2)]`. If custom observables (`observables`) are - given, the indices are ordered in the same manner as the final index of - `ops`. + cross-correlations generated by all `observables`. + To retain only the xx and xy correlations, one would set + `correlations=[(:Sx,:Sx), (:Sx,:Sy)]` or `correlations=[(1,1),(1,2)]`. """ function DynamicStructureFactor(sys::System; Δt, nω, ωmax, kwargs...) nω = Int64(nω) @@ -183,7 +291,7 @@ Creates a `StructureFactor` object for calculating and storing instantaneous structure factor intensities ``𝒮(𝐪)``. This data will be calculated from the spin-spin correlations of equilibrium snapshots, absent any dynamical information. ``𝒮(𝐪)`` data can be retrieved by calling -[`instant_intensities`](@ref). +[`instant_intensities_interpolated`](@ref). _Important note_: When dealing with continuous (non-Ising) spins, consider creating a full [`DynamicStructureFactor`](@ref) object instead of an @@ -191,7 +299,7 @@ creating a full [`DynamicStructureFactor`](@ref) object instead of an which ``𝒮(𝐪)`` can be obtained by integrating out ``ω``. During this integration step, Sunny can incorporate temperature- and ``ω``-dependent classical-to-quantum correction factors to produce more accurate ``𝒮(𝐪)`` -estimates. See [`instant_intensities`](@ref) for more information. +estimates. See [`instant_intensities_interpolated`](@ref) for more information. Prior to calling `InstantStructureFactor`, ensure that `sys` represents a good equilibrium sample. Additional sample data may be accumulated by calling @@ -204,18 +312,17 @@ The following optional keywords are available: `:symmetrize`. The latter will symmetrize the trajectory in time, which can be useful for removing Fourier artifacts that arise when calculating the correlations. -- `observables`: Enables an advanced feature for SU(_N_) mode, allowing the user - to specify custom observables other than the three components of the dipole. - To use this features, `observables` must be given an `N×N×numops` array, - where the final index is used to retrieve each `N×N` operator. +- `observables`: Allows the user to specify custom observables. The `observables` + must be given as a list of complex `N×N` matrices or `LinearMap`s. It's + recommended to name each observable, for example: + `observables = [:A => a_observable_matrix, :B => b_map, ...]`. + By default, Sunny uses the 3 components of the dipole, `:Sx`, `:Sy` and `:Sz`. - `correlations`: Specify which correlation functions are calculated, i.e. which matrix elements ``αβ`` of ``𝒮^{αβ}(q,ω)`` are calculated and stored. Specified with a vector of tuples. By default Sunny records all auto- and - cross-correlations generated by the x, y, and z dipolar components (1, 2, - and 3 respectively). To retain only the xx and xy correlations, one would - set `correlations=[(1,1), (1,2)]`. If custom observables (`observables`) are - given, the indices are ordered in the same manner as the final index of - `observables`. + cross-correlations generated by all `observables`. + To retain only the xx and xy correlations, one would set + `correlations=[(:Sx,:Sx), (:Sx,:Sy)]` or `correlations=[(1,1),(1,2)]`. """ function InstantStructureFactor(sys::System; kwargs...) StructureFactor(sys; Δt=0.1, nω=1, measperiod=1, kwargs...) diff --git a/src/Sunny.jl b/src/Sunny.jl index 8950d861f..e10696613 100644 --- a/src/Sunny.jl +++ b/src/Sunny.jl @@ -1,6 +1,7 @@ module Sunny using LinearAlgebra +using LinearMaps import StaticArrays: SVector, SMatrix, SArray, MVector, MMatrix, SA, @SVector import Requires: @require import OffsetArrays: OffsetArray, OffsetMatrix, Origin @@ -10,7 +11,7 @@ import ProgressMeter: Progress, next! import Printf: @printf, @sprintf import Random: Random, randn! import DynamicPolynomials as DP -import DataStructures: SortedDict +import DataStructures: SortedDict, OrderedDict import Optim # Specific to Symmetry/ @@ -94,12 +95,13 @@ include("StructureFactors/BasisReduction.jl") include("StructureFactors/Interpolation.jl") include("StructureFactors/PowderAveraging.jl") include("StructureFactors/DataRetrieval.jl") +include("SpinWaveTheory/Intensities.jl") # Requires ElementContraction export DynamicStructureFactor, InstantStructureFactor, StructureFactor, FormFactor, - add_sample!, intensities, instant_intensities, broaden_energy, lorentzian, + add_sample!, intensities_interpolated, instant_intensities_interpolated, broaden_energy, lorentzian, connected_path, all_exact_wave_vectors, ωs, spherical_shell, merge!, BinningParameters, integrate_axes!, unit_resolution_binning_parameters, rlu_to_absolute_units!, intensities_binned, one_dimensional_cut_binning_parameters, axes_bincenters, - connected_path_bins, powder_averaged_bins + connected_path_bins, powder_averaged_bins, intensity_formula, integrated_lorentzian, count_bins include("SunnyGfx/SunnyGfx.jl") include("SunnyGfx/CrystalViewer.jl") diff --git a/src/System/System.jl b/src/System/System.jl index 384e01d08..2c2829ba2 100644 --- a/src/System/System.jl +++ b/src/System/System.jl @@ -75,6 +75,22 @@ function System(crystal::Crystal, latsize::NTuple{3,Int}, infos::Vector{SpinInfo return ret end +function Base.show(io::IO, sys::System{N}) where N + modename = if sys.mode==:SUN + "SU($N)" + elseif sys.mode==:dipole + "Dipole" + elseif sys.mode==:large_S + "Large-S" + else + error("Unreachable") + end + print(io,"System{$modename}[$(sys.latsize)×$(natoms(sys.crystal))]") + if !isnothing(sys.origin) + print(io,"[Reshape = $(cell_dimensions(sys))]") + end +end + function Base.show(io::IO, ::MIME"text/plain", sys::System{N}) where N modename = if sys.mode==:SUN "SU($N)" @@ -86,7 +102,7 @@ function Base.show(io::IO, ::MIME"text/plain", sys::System{N}) where N error("Unreachable") end printstyled(io, "System [$modename]\n"; bold=true, color=:underline) - println(io, "Cell size $(natoms(sys.crystal)), Lattice size $(sys.latsize)") + println(io, "Lattice: $(sys.latsize)×$(natoms(sys.crystal))") if !isnothing(sys.origin) println(io, "Reshaped cell geometry $(cell_dimensions(sys))") end diff --git a/test/test_structure_factors.jl b/test/test_structure_factors.jl index 37f0b3aee..f4c737c3f 100644 --- a/test/test_structure_factors.jl +++ b/test/test_structure_factors.jl @@ -27,10 +27,10 @@ sys = simple_model_sf(; mode=:SUN) thermalize_simple_model!(sys; kT=0.1) S = Sunny.spin_matrices(2) - observables = cat(S...; dims=3) + observables = Dict(:Sx => S[1], :Sy => S[2], :Sz => S[3]) sf = DynamicStructureFactor(sys; nω=100, ωmax=10.0, Δt=0.1, apply_g=false, observables) qgrid = all_exact_wave_vectors(sf) - vals = intensities(sf, qgrid, :trace; negative_energies=true) + vals = intensities_interpolated(sf, qgrid; formula = intensity_formula(sf,:trace),negative_energies=true) total_intensity_trace = sum(vals) @test isapprox(total_intensity_trace / prod(sys.latsize), 1.0; atol=1e-12) @@ -40,39 +40,42 @@ thermalize_simple_model!(sys; kT=0.1) sf = DynamicStructureFactor(sys; Δt=0.1, nω=100, ωmax=10.0, apply_g=false) add_sample!(sf, sys) - vals = intensities(sf, qgrid, :trace; negative_energies=true) + formula = intensity_formula(sf,:trace) + vals = intensities_interpolated(sf, qgrid; formula, negative_energies=true) total_intensity_trace = sum(vals) @test isapprox(total_intensity_trace / prod(sys.latsize), 1.0; atol=1e-12) # Test perp reduces intensity - vals = intensities(sf, qgrid, :perp; negative_energies=true) + vals = intensities_interpolated(sf, qgrid; negative_energies=true) # Formula is :perp by default total_intensity_unpolarized = sum(vals) @test total_intensity_unpolarized < total_intensity_trace # Test diagonal elements are approximately real (at one wave vector) - for α in 1:3 - intensities_symmetric = intensities(sf, [(0.25, 0.5, 0)], (α,α)) - @test sum(imag(intensities_symmetric)) < 1e-15 + diag_elems = [(α,α) for α in keys(sf.observable_ixs)] + formula_imaginary_parts = intensity_formula(sf,diag_elems) do k,ω,corr + sum(abs.(imag.(corr))) end + intensities_symmetric = intensities_interpolated(sf, [(0.25, 0.5, 0)]; formula = formula_imaginary_parts) + @test sum(imag(intensities_symmetric)) < 1e-15 # Test form factor correction works and is doing something. ddtodo: add example with sublattice formfactors = [FormFactor(1, "Fe2")] - vals = intensities(sf, qgrid, :trace; formfactors, negative_energies=true) + vals = intensities_interpolated(sf, qgrid; formula = intensity_formula(sf,:trace;formfactors = formfactors), negative_energies=true) total_intensity_ff = sum(vals) @test total_intensity_ff != total_intensity_trace # Test path function and interpolation working (no correctness implied here) qs, _ = connected_path(sf, [(0, 0, 0), (0, 1, 0), (1, 1, 0)], 20) - vals = intensities(sf, qs, :trace; interpolation=:linear) + vals = intensities_interpolated(sf, qs; formula, interpolation=:linear) @test length(size(vals)) == 2 # Test static from dynamic intensities working - static_vals = instant_intensities(sf, qgrid, :trace; negative_energies=true) + static_vals = instant_intensities_interpolated(sf, qgrid; formula, negative_energies=true) total_intensity_static = sum(static_vals) @test isapprox(total_intensity_static, total_intensity_trace; atol=1e-12) # Order of summation can lead to very small discrepancies @@ -80,7 +83,7 @@ sys = simple_model_sf(; mode=:dipole) thermalize_simple_model!(sys; kT=0.1) sf = InstantStructureFactor(sys; apply_g=false) - true_static_vals = instant_intensities(sf, qgrid, :trace) + true_static_vals = instant_intensities_interpolated(sf, qgrid; formula = intensity_formula(sf,:trace)) true_static_total = sum(true_static_vals) @test isapprox(true_static_total / prod(sys.latsize), 1.0; atol=1e-12) @@ -114,7 +117,7 @@ end sf = DynamicStructureFactor(sys; nω=25, ωmax=5.5, Δt=2Δt_langevin) qs, _ = connected_path(sf, [(0.0, 0.0, 0.0), (0.5, 0.0, 0.0), (0.5, 0.5, 0.0), (0.0, 0.0, 0.0)], 10) - data = intensities(sf, qs, :trace; interpolation=:linear, kT) + data = intensities_interpolated(sf, qs; formula = intensity_formula(sf,:trace;kT = kT), interpolation=:linear) refdata = [4.254958223663417 -6.420612698722547e-15 9.517855303472442e-15 2.254439019145773e-15 -7.411655927494534e-15 -3.415475412410896e-15 7.193291213800562e-16 9.03543153937406e-16 -9.219208521726188e-16 3.3188559986441705e-16 -4.05636043338745e-15 0.0 3.844317340337979e-15 -2.6066681841577996e-15 -2.839447505217288e-15 2.056986106201179e-15 -1.475037653587726e-15 9.403365037570114e-16 -5.807960757836408e-16 -8.320134259982439e-16 0.0 4.839967298001723e-16 3.1943784166793436e-15 -1.6962933006317804e-15 -9.403365036110586e-16; 4.050322690710251 0.001503739876456909 0.0011921434008736649 0.0017936849615682003 0.0008070642362539781 0.0019841446952351166 0.003333290236633009 0.01578638047521035 0.2065064395838012 0.2139565139050153 0.008695481317449428 0.0032619849411028076 0.002433260632334384 0.0022639058709873396 0.0014892407581326237 0.0013016014987539154 0.0016154825324611825 0.0012544930922307295 0.0011056948071601774 0.002783872916314107 0.0009854416043813062 0.0017992116724471591 0.000683853705951668 0.002129274756212289 0.0011434575660303405; 3.8456871577570855 0.0030074797529202387 0.002384286801737812 0.003587369923134146 0.001614128472515368 0.003968289390473648 0.0066665804732652986 0.0315727609504198 0.4130128791676033 0.42791302781003027 0.01739096263490291 0.006523969882205615 0.004866521264664924 0.0045278117419772865 0.0029784815162680867 0.002603202997505774 0.00323096506492384 0.0025089861844605184 0.0022113896143209355 0.005567745832629046 0.0019708832087626123 0.0035984233448938343 0.0013677074119001417 0.0042585495124262745 0.0022869151320616217; 3.6410516248039193 0.004511219629383568 0.0035764302026019597 0.005381054884700092 0.002421192708776758 0.005952434085712181 0.009999870709897588 0.047359141425629234 0.6195193187514055 0.6418695417150454 0.0260864439523564 0.009785954823308423 0.007299781896995466 0.0067917176129672335 0.004467722274403551 0.003904804496257633 0.004846447597386498 0.003763479276690309 0.0033170844214816943 0.008351618748943985 0.0029563248131439187 0.00539763501734051 0.0020515611178486158 0.00638782426864026 0.003430372698092903; 3.4364160918507536 0.006014959505846898 0.004768573603466106 0.007174739846266038 0.0032282569450381477 0.007936578780950712 0.013333160946529877 0.06314552190083869 0.8260257583352075 0.8558260556200602 0.03478192526980988 0.01304793976441123 0.009733042529326004 0.009055623483957179 0.005956963032539013 0.005206405995009491 0.006461930129849155 0.005017972368920097 0.004422779228642452 0.011135491665258923 0.003941766417525225 0.007196846689787185 0.002735414823797089 0.008517099024854244 0.0045738302641241835; 3.231780558897588 0.007518699382310227 0.005960717004330253 0.008968424807831984 0.004035321181299538 0.009920723476189244 0.016666451183162168 0.07893190237604812 1.0325321979190099 1.069782569525075 0.043477406587263366 0.01630992470551404 0.012166303161656546 0.011319529354947126 0.007446203790674477 0.00650800749376135 0.008077412662311812 0.006272465461149887 0.0055284740358032105 0.013919364581573862 0.004927208021906531 0.00899605836223386 0.003419268529745563 0.010646373781068232 0.005717287830155465; 3.0271450259444226 0.009022439258773557 0.007152860405194402 0.01076210976939793 0.004842385417560928 0.011904868171427777 0.019999741419794458 0.09471828285125757 1.239038637502812 1.2837390834300904 0.05217288790471685 0.019571909646616845 0.014599563793987087 0.013583435225937074 0.00893544454880994 0.007809608992513209 0.009692895194774473 0.0075269585533796765 0.00663416884296397 0.016703237497888802 0.005912649626287837 0.010795270034680535 0.004103122235694037 0.012775648537282216 0.006860745396186747; 2.822509492991256 0.010526179135236886 0.008345003806058548 0.012555794730963877 0.005649449653822317 0.013889012866666309 0.023333031656426747 0.11050466332646702 1.445545077086614 1.4976955973351054 0.060868369222170336 0.022833894587719655 0.01703282442631763 0.01584734109692702 0.010424685306945404 0.009111210491265068 0.01130837772723713 0.008781451645609466 0.007739863650124727 0.019487110414203743 0.006898091230669145 0.012594481707127212 0.00478697594164251 0.014904923293496202 0.00800420296221803; 2.6178739600380907 0.012029919011700217 0.009537147206922695 0.014349479692529822 0.006456513890083707 0.01587315756190484 0.026666321893059036 0.12629104380167647 1.652051516670416 1.7116521112401202 0.06956385053962381 0.02609587952882246 0.019466085058648164 0.018111246967916966 0.011913926065080866 0.010412811990016925 0.012923860259699785 0.010035944737839253 0.008845558457285485 0.02227098333051868 0.00788353283505045 0.014393693379573885 0.005470829647590983 0.017034198049710184 0.009147660528249307; 2.4132384270849245 0.013533658888163548 0.010729290607786842 0.01614316465409577 0.007263578126345098 0.017857302257143375 0.02999961212969133 0.1420774242768859 1.8585579562542185 1.9256086251451354 0.07825933185707731 0.02935786446992527 0.021899345690978712 0.020375152838906913 0.013403166823216331 0.011714413488768786 0.014539342792162445 0.011290437830069045 0.009951253264446245 0.025054856246833623 0.008868974439431757 0.016192905052020565 0.006154683353539458 0.019163472805924175 0.01029111809428059; 2.2086028941317593 0.015037398764626876 0.01192143400865099 0.017936849615661714 0.008070642362606488 0.019841446952381905 0.03333290236632362 0.15786380475209535 2.0650643958380206 2.1395651390501498 0.08695481317453078 0.03261984941102808 0.024332606323309247 0.022639058709896857 0.014892407581351793 0.013016014987520643 0.016154825324625102 0.012544930922298834 0.011056948071607002 0.027838729163148557 0.009854416043813062 0.017992116724467235 0.0068385370594879315 0.02129274756213816 0.01143457566031187; 2.003967361178593 0.016541138641090206 0.013113577409515138 0.019730534577227662 0.008877706598867878 0.021825591647620442 0.03666619260295591 0.1736501852273048 2.271570835421823 2.3535216529551652 0.09565029449198428 0.03588183435213089 0.02676586695563979 0.02490296458088681 0.016381648339487256 0.014317616486272505 0.01777030785708776 0.013799424014528626 0.012162642878767761 0.0306226020794635 0.010839857648194371 0.019791328396913912 0.007522390765436406 0.023422022318352144 0.012578033226343153; 1.7993318282254276 0.018044878517553535 0.014305720810379285 0.021524219538793604 0.009684770835129267 0.02380973634285897 0.039999482839588194 0.18943656570251424 2.478077275005625 2.5674781668601803 0.10434577580943775 0.03914381929323369 0.02919912758797033 0.027166870451876754 0.01787088909762272 0.01561921798502436 0.01938579038955042 0.015053917106758413 0.01326833768592852 0.03340647499577844 0.011825299252575675 0.02159054006936059 0.008206244471384879 0.025551297074566128 0.013721490792374435; 1.5946962952722614 0.019548618394016865 0.015497864211243433 0.023317904500359556 0.010491835071390659 0.025793881038097505 0.043332773076220486 0.20522294617772371 2.684583714589427 2.7814346807651953 0.11304125712689125 0.04240580423433651 0.03163238822030088 0.0294307763228667 0.019360129855758183 0.01692081948377622 0.021001272922013075 0.016308410198988207 0.01437403249308928 0.036190347912093385 0.012810740856956984 0.023389751741807265 0.008890098177333354 0.02768057183078012 0.014864948358405716; 1.3900607623190957 0.021052358270480194 0.01669000761210758 0.0251115894619255 0.011298899307652047 0.027778025733336032 0.04666606331285277 0.22100932665293316 2.891090154173229 2.9953911946702103 0.12173673844434474 0.04566778917543931 0.034065648852631414 0.03169468219385665 0.020849370613893647 0.01822242098252808 0.022616755454475734 0.017562903291217992 0.015479727300250035 0.03897422082840832 0.01379618246133829 0.02518896341425394 0.009573951883281827 0.0298098465869941 0.016008405924437; 1.18542522936593 0.022556098146943523 0.017882151012971727 0.026905274423491443 0.012105963543913438 0.029762170428574562 0.04999935354948506 0.23679570712814257 3.097596593757031 3.209347708575225 0.1304322197617982 0.04892977411654211 0.03649890948496195 0.03395858806484659 0.02233861137202911 0.019524022481279938 0.02423223798693839 0.01881739638344778 0.016585422107410795 0.04175809374472325 0.014781624065719595 0.02698817508670061 0.0102578055892303 0.03193912134320808 0.017151863490468277; 0.9807896964127646 0.024059838023406853 0.01907429441383587 0.02869895938505739 0.012913027780174825 0.031746315123813096 0.05333264378611735 0.25258208760335205 3.304103033340833 3.42330422248024 0.13912770107925168 0.05219175905764492 0.03893217011729248 0.036222493935836535 0.02382785213016457 0.020825623980031793 0.025847720519401044 0.020071889475677566 0.01769111691457155 0.04454196666103819 0.0157670656701009 0.02878738675914729 0.010941659295178773 0.03406839609942207 0.018295321056499555; 0.776154163459599 0.025563577899870182 0.020266437814700018 0.03049264434662333 0.013720092016436215 0.03373045981905163 0.05666593402274964 0.26836846807856146 3.510609472924635 3.6372607363852545 0.14782318239670517 0.05545374399874772 0.041365430749623025 0.03848639980682648 0.025317092888300034 0.022127225478783652 0.027463203051863703 0.021326382567907354 0.018796811721732306 0.047325839577353125 0.016752507274482204 0.03058659843159396 0.011625513001127246 0.036197670855636056 0.019438778622530836; 0.5715186305064324 0.027067317776333515 0.021458581215564168 0.03228632930818928 0.014527156252697609 0.03571460451429016 0.059999224259381936 0.28415484855377093 3.717115912508438 3.8512172502902704 0.15651866371415868 0.05871572893985054 0.043798691381953574 0.040750305677816435 0.0268063336464355 0.023428826977535518 0.029078685584326365 0.02258087566013715 0.01990250652889307 0.05010971249366808 0.017737948878863513 0.032385810104040645 0.012309366707075723 0.038326945611850044 0.02058223618856212; 0.36688309755326676 0.028571057652796845 0.022650724616428315 0.03408001426975523 0.015334220488958996 0.03769874920952869 0.06333251449601422 0.29994122902898035 3.9236223520922398 4.065173764195285 0.16521414503161214 0.06197771388095335 0.046231952014284115 0.04301421154880638 0.028295574404570964 0.024730428476287373 0.030694168116789024 0.02383536875236694 0.021008201336053824 0.05289358540998301 0.01872339048324482 0.034185021776487315 0.012993220413024196 0.040456220368064025 0.021725693754593402; 0.16224756460010117 0.03007479752926017 0.02384286801729246 0.03587369923132117 0.016141284725220386 0.03968289390476722 0.06666580473264652 0.3157276095041898 4.130128791676042 4.2791302781002996 0.17390962634906562 0.06523969882205616 0.04866521264661465 0.04527811741979632 0.029784815162706427 0.026032029975039228 0.032309650649251676 0.025089861844596727 0.022113896143214587 0.055677458326297946 0.019708832087626124 0.035984233448933985 0.013677074118972668 0.04258549512427801 0.02286915132062468; 0.02547366951406332 0.030608915692135863 0.024271147530223282 0.03662608658488363 0.01676090815809476 0.04064434266479797 0.06849494958147812 0.32145111202501353 4.199360924209134 4.3506321884138615 0.17752645496153693 0.06722822459047063 0.051291481078470945 0.05132773303506229 0.3088833337771953 0.061862859955905555 0.03518424091438777 0.026375405633927505 0.024374059691704853 0.05735620206994079 0.02029650166875649 0.03713960926899464 0.014297916726271169 0.043563776598251065 0.02400154960357742; 0.024423050162280906 0.02920379042783628 0.023171699267287443 0.03529587872243911 0.017005649984195115 0.03956039955441315 0.06731580365470832 0.30704885863706577 3.9940444426408024 4.137224891544513 0.17098597816404387 0.0666698320135085 0.054303765109378746 0.06494876813888033 1.1432004081043936 0.1667521469010012 0.040577046644870746 0.027723050817460265 0.028943160722854154 0.05682468746823945 0.020088627203384956 0.03700731338428324 0.014792737136269728 0.04224007150774224 0.02511182932037308; 0.023372430810498498 0.0277986651635367 0.02207225100435161 0.03396567085999459 0.017250391810295473 0.03847645644402834 0.06613665772793854 0.29264660524911806 3.7887279610724724 3.9238175946751657 0.1644455013665509 0.06611143943654638 0.05731604914028654 0.07856980324269833 1.9775174824315893 0.27164143384609646 0.04596985237535371 0.02907069600099302 0.033512261754003445 0.056293172866538115 0.01988075273801343 0.03687501749957184 0.015287557546268288 0.040916366417233424 0.02622210903716874; 0.022321811458716086 0.026393539899237126 0.020972802741415775 0.03263546299755008 0.017495133636395824 0.03739251333364352 0.06495751180116877 0.27824435186117036 3.583411479504143 3.7104102978058187 0.15790502456905786 0.06555304685958423 0.060328333171194334 0.09219083834651633 2.811834556758785 0.3765307207911917 0.05136265810583667 0.030418341184525778 0.03808136278515273 0.055761658264836776 0.019672878272641903 0.03674272161486043 0.015782377956266847 0.039592661326724614 0.027332388753964397; 0.02127119210693368 0.024988414634937547 0.019873354478479946 0.03130525513510556 0.01773987546249618 0.0363085702232587 0.06377836587439897 0.26384209847322265 3.3780949979358135 3.497003000936471 0.15136454777156488 0.0649946542826221 0.06334061720210213 0.10581187345033433 3.646151631085981 0.481420007736287 0.05675546383631963 0.03176598636805854 0.04265046381630202 0.055230143663135445 0.019465003807270377 0.036610425730149034 0.016277198366265407 0.0382689562362158 0.028442668470760057; 0.02022057275515127 0.023583289370637965 0.01877390621554411 0.029975047272661048 0.017984617288596533 0.03522462711287389 0.06259921994762921 0.24943984508527486 3.172778516367483 3.283595704067124 0.14482407097407182 0.06443626170565997 0.06635290123300994 0.11943290855415237 4.4804687054131795 0.5863092946813826 0.06214826956680262 0.033113631551591305 0.04721956484745132 0.054698629061434106 0.019257129341898847 0.036478129845437635 0.01677201877626397 0.03694525114570699 0.02955294818755572; 0.019169953403368863 0.022178164106338386 0.017674457952608275 0.028644839410216533 0.018229359114696888 0.034140684002489076 0.06142007402085942 0.23503759169732716 2.9674620347991536 3.0701884071977767 0.13828359417657882 0.06387786912869783 0.06936518526391774 0.13305394365797038 5.314785779740375 0.691198581626478 0.06754107529728559 0.03446127673512406 0.051788665878600613 0.05416711445973277 0.01904925487652732 0.036345833960726236 0.017266839186262525 0.03562154605519817 0.030663227904351376; 0.01811933405158645 0.020773038842038807 0.016575009689672442 0.02731463154777202 0.018474100940797242 0.03305674089210426 0.06024092809408964 0.22063533830937945 2.7621455532308237 2.8567811103284297 0.1317431173790858 0.0633194765517357 0.07237746929482552 0.14667497876178837 6.1491028540675705 0.7960878685715732 0.07293388102776854 0.03580892191865682 0.0563577669097499 0.05363559985803143 0.018841380411155795 0.03621353807601484 0.01776165959626108 0.03429784096468936 0.031773507621147036; 0.01706871469980404 0.019367913577739232 0.015475561426736608 0.025984423685327506 0.018718842766897593 0.03197279778171945 0.05906178216731986 0.20623308492143172 2.556829071662494 2.6433738134590823 0.1252026405815928 0.06276108397477356 0.07538975332573333 0.16029601386560638 6.983419928394766 0.9009771555166685 0.0783266867582515 0.03715656710218958 0.06092686794089919 0.05310408525633009 0.018633505945784265 0.03608124219130343 0.018256480006259644 0.03297413587418054 0.03288378733794269; 0.016018095348021632 0.01796278831343965 0.014376113163800774 0.024654215822882988 0.01896358459299795 0.03088885467133463 0.05788263624055008 0.191830831533484 2.351512590094164 2.4299665165897353 0.11866216378409977 0.062202691397811434 0.07840203735664114 0.17391704896942442 7.817737002721964 1.005866442461764 0.08371949248873449 0.038504212285722345 0.0654959689720485 0.05257257065462877 0.018425631480412742 0.03594894630659203 0.018751300416258204 0.031650430783671725 0.033994067054738356; 0.014967475996239217 0.01655766304914007 0.013276664900864937 0.02332400796043847 0.019208326419098302 0.029804911560949814 0.05670349031378029 0.17742857814553625 2.146196108525834 2.2165592197203874 0.11212168698660675 0.0616442988208493 0.08141432138754894 0.18753808407324243 8.652054077049161 1.1107557294068595 0.08911229821921746 0.0398518574692551 0.07006507000319778 0.052041056052927415 0.018217757015041212 0.03581665042188063 0.019246120826256763 0.030326725693162908 0.03510434677153401; 0.014436171477102593 0.01682110165656617 0.01384019039687941 0.02274678870906014 0.018926108021738303 0.029215191636604843 0.05459648496127452 0.16980955833438244 2.043332872550263 2.111326979767518 0.11430303417450743 0.12084061140332855 0.08845069039822534 0.1909162010813912 8.428096226730217 1.4943423937097053 0.11891818265150902 0.05149989182766467 0.07213667095360889 0.053481956722232236 0.02042586397586009 0.036505637914788676 0.019789529645368776 0.03159154617870763 0.03754997948150723; 0.013904866957965967 0.017084540263992268 0.014403715892893881 0.022169569457681816 0.018643889624378306 0.028625471712259876 0.052489479608768735 0.16219053852322862 1.9404696365746914 2.006094739814648 0.1164843813624081 0.18003692398580778 0.09548705940890173 0.19429431808953995 8.204138376411272 1.877929058012551 0.14872406708380057 0.06314792618607425 0.07420827190401999 0.05492285739153705 0.022633970936678967 0.03719462540769671 0.020332938464480785 0.032856366664252346 0.03999561219148044; 0.013373562438829341 0.017347978871418368 0.014967241388908354 0.021592350206303488 0.01836167122701831 0.028035751787914905 0.050382474256262966 0.1545715187120748 1.8376064005991202 1.9008624998617778 0.1186657285503088 0.23923323656828704 0.10252342841957814 0.19767243509768873 7.980180526092329 2.261515722315397 0.17852995151609213 0.07479596054448381 0.07627987285443111 0.05636375806084187 0.024842077897497844 0.03788361290060475 0.0208763472835928 0.03412118714979707 0.04244124490145365; 0.012842257919692714 0.01761141747884447 0.015530766884922827 0.021015130954925163 0.01807945282965831 0.027446031863569937 0.04827546890375718 0.146952498900921 1.7347431646235492 1.7956302599089078 0.12084707573820946 0.29842954915076625 0.10955979743025453 0.20105055210583747 7.756222675773384 2.6451023866182424 0.20833583594838367 0.0864439949028934 0.0783514738048422 0.05780465873014669 0.02705018485831672 0.03857260039351278 0.021419756102704805 0.03538600763534179 0.04488687761142686; 0.01231095340055609 0.01787485608627057 0.016094292380937302 0.020437911703546835 0.017797234432298314 0.02685631193922497 0.04616846355125141 0.13933347908976718 1.631879928647978 1.690398019956038 0.12302842292611015 0.3576258617332455 0.11659616644093093 0.20442866911398627 7.53226482545444 3.028689050921088 0.23814172038067527 0.09809202926130296 0.08042307475525332 0.05924555939945152 0.0292582918191356 0.039261587886420826 0.021963164921816818 0.036650828120886506 0.04733251032140008; 0.011779648881419465 0.01813829469369667 0.016657817876951773 0.01986069245216851 0.017515016034938317 0.026266592014880002 0.04406145819874564 0.1317144592786134 1.529016692672407 1.5851657800031684 0.12520977011401083 0.41682217431572477 0.12363253545160735 0.20780678612213505 7.308306975135497 3.4122757152239345 0.2679476048129668 0.10974006361971254 0.0824946757056644 0.06068646006875633 0.03146639877995448 0.03995057537932886 0.02250657374092883 0.037915648606431235 0.04977814303137329; 0.01124834436228284 0.018401733301122766 0.017221343372966248 0.019283473200790182 0.017232797637578318 0.025676872090535027 0.041954452846239866 0.12409543946745955 1.4261534566968355 1.4799335400502984 0.12739111730191152 0.47601848689820403 0.13066890446228374 0.21118490313028376 7.084349124816551 3.79586237952678 0.2977534892452584 0.12138809797812213 0.08456627665607552 0.06212736073806115 0.03367450574077335 0.0406395628722369 0.02304998256004084 0.03918046909197595 0.052223775741346505; 0.010717039843146214 0.018665171908548862 0.01778486886898072 0.018706253949411854 0.01695057924021832 0.02508715216619006 0.03984744749373409 0.11647641965630574 1.3232902207212645 1.3747013000974286 0.1295724644898122 0.5352147994806832 0.13770527347296013 0.21456302013843254 6.860391274497608 4.179449043829625 0.3275593736775499 0.1330361323365317 0.08663787760648661 0.06356826140736596 0.03588261270159223 0.041328550365144934 0.02359339137915285 0.040445289577520666 0.05466940845131972; 0.010185735324009588 0.018928610515974963 0.01834839436499519 0.01812903469803353 0.01666836084285832 0.024497432241845092 0.03774044214122831 0.10885739984515191 1.220426984745693 1.2694690601445584 0.13175381167771288 0.5944111120631626 0.14474164248363652 0.2179411371465813 6.636433424178663 4.563035708132472 0.35736525810984154 0.14468416669494127 0.08870947855689773 0.06500916207667079 0.03809071966241111 0.04201753785805298 0.024136800198264863 0.04171011006306539 0.057115041161292936; 0.009654430804872964 0.019192049123401063 0.018911919861009666 0.017551815446655204 0.016386142445498325 0.02390771231750012 0.035633436788722536 0.1012383800339981 1.117563748770122 1.1642368201916888 0.13393515886561355 0.6536074246456418 0.15177801149431291 0.22131925415473008 6.4124755738597194 4.946622372435317 0.3871711425421331 0.15633220105335083 0.09078107950730883 0.0664500627459756 0.040298826623229984 0.04270652535096102 0.024680209017376872 0.04297493054861011 0.059560673871266144; 0.009123126285736338 0.019455487730827163 0.019475445357024137 0.016974596195276873 0.016103924048138325 0.02331799239315515 0.03352643143621676 0.09361936022284428 1.0147005127945508 1.0590045802388188 0.13611650605351422 0.7128037372281211 0.15881438050498933 0.22469737116287886 6.188517723540775 5.330209036738163 0.4169770269744246 0.16798023541176044 0.09285268045771994 0.06789096341528042 0.04250693358404887 0.043395512843869055 0.025223617836488885 0.04423975103415483 0.06200630658123936; 0.008591821766599712 0.019718926338253263 0.020038970853038608 0.016397376943898548 0.01582170565077833 0.022728272468810182 0.03141942608371098 0.08600034041169047 0.9118372768189795 0.953772340285949 0.13829785324141491 0.7720000498106003 0.16585074951566575 0.22807548817102763 5.96455987322183 5.7137957010410085 0.44678291140671617 0.17962826977017 0.09492428140813104 0.06933186408458525 0.044715040544867746 0.04408450033677709 0.0257670266556009 0.045504571519699555 0.06445193929121257; 0.008060517247463086 0.019982364945679364 0.020602496349053083 0.015820157692520223 0.015539487253418333 0.022138552544465215 0.029312420731205213 0.07838132060053665 0.8089740408434083 0.848540100333079 0.1404792004293156 0.8311963623930796 0.17288711852634214 0.23145360517917637 5.740602022902886 6.097382365343855 0.47658879583900776 0.19127630412857957 0.09699588235854216 0.07077276475389008 0.04692314750568662 0.04477348782968513 0.026310435474712908 0.04676939200524427 0.06689757200118579; 0.007529212728326461 0.02024580355310546 0.021166021845067554 0.015242938441141895 0.015257268856058335 0.021548832620120244 0.027205415378699437 0.07076230078938284 0.7061108048678372 0.7433078603802092 0.14266054761721625 0.8903926749755587 0.1799234875370185 0.23483172218732512 5.516644172583941 6.4809690296467 0.5063946802712993 0.20292433848698915 0.09906748330895325 0.07221366542319488 0.0491312544665055 0.04546247532259317 0.026853844293824917 0.048034212490788986 0.069343204711159; 0.006997908209189836 0.02050924216053156 0.021729547341082026 0.014665719189763568 0.014975050458698336 0.020959112695775276 0.025098410026193667 0.06314328097822904 0.603247568892266 0.6380756204273396 0.14484189480511694 0.9495889875580379 0.18695985654769492 0.23820983919547392 5.292686322264998 6.864555693949546 0.5362005647035908 0.2145723728453987 0.10113908425936437 0.07365456609249971 0.05133936142732438 0.046151462815501205 0.02739725311293693 0.04929903297633371 0.0717888374211322; 0.006466603690053211 0.020772680767957657 0.0222930728370965 0.014088499938385244 0.01469283206133834 0.02036939277143031 0.02299140467368789 0.05552426116707522 0.5003843329166949 0.5328433804744697 0.14702324199301764 1.0087853001405171 0.19399622555837132 0.24158795620362267 5.068728471946054 7.248142358252391 0.5660064491358824 0.2262204072038083 0.10321068520977546 0.07509546676180452 0.05354746838814325 0.04684045030840924 0.027940661932048937 0.05056385346187843 0.07423447013110542; 0.005935299170916586 0.02103611937538376 0.022856598333110975 0.013511280687006917 0.014410613663978342 0.019779672847085338 0.020884399321182117 0.04790524135592142 0.3975210969411238 0.4276111405216 0.1492045891809183 1.0679816127229964 0.2010325945690477 0.24496607321177144 4.844770621627109 7.631729022555237 0.5958123335681739 0.23786844156221781 0.10528228616018658 0.07653636743110934 0.055755575348962125 0.04752943780131728 0.02848407075116095 0.051828673947423146 0.07668010284107864; 0.005403994651779959 0.021299557982809858 0.023420123829125446 0.012934061435628587 0.014128395266618342 0.019189952922740366 0.018777393968676337 0.040286221544767575 0.29465786096555224 0.3223789005687297 0.15138593636881897 1.127177925305476 0.20806896357972413 0.2483441902199202 4.620812771308164 8.015315686858084 0.6256182180004656 0.24951647592062745 0.10735388711059768 0.07797726810041417 0.05796368230978101 0.04821842529422532 0.029027479570272963 0.05309349443296787 0.07912573555105186; 0.0048726901326433345 0.02156299659023596 0.023983649325139918 0.012356842184250262 0.013846176869258346 0.0186002329983954 0.016670388616170564 0.03266720173361377 0.19179462498998115 0.21714666061585997 0.15356728355671967 1.186374237887955 0.21510533259040052 0.251722307228069 4.39685492098922 8.39890235116093 0.6554241024327572 0.26116451027903703 0.1094254880610088 0.07941816876971898 0.060171789270599886 0.048907412787133356 0.029570888389384972 0.05435831491851259 0.08157136826102507; 0.004341385613506709 0.021826435197662055 0.024547174821154392 0.011779622932871936 0.01356395847189835 0.01801051307405043 0.014563383263664789 0.02504818192245996 0.08893138901441007 0.11191442066299019 0.15574863074462034 1.2455705504704342 0.22214170160107694 0.25510042423621776 4.172897070670277 8.782489015463774 0.6852299868650487 0.2728125446374466 0.11149708901141989 0.08085906943902381 0.062379896231418756 0.0495964002800414 0.030114297208496982 0.055623135404057306 0.08401700097099828; 0.0040615115184178055 0.0218995017937112 0.0248178406041569 0.012447603611329593 0.014006174086251027 0.018088916017701247 0.013457842968559222 0.020585755770862534 0.020443597700697137 0.04191655754350475 0.15538767829138575 1.265703925957692 0.2256843047879328 0.25777441201324325 3.982241209196207 9.015864926829838 0.9209236273137634 0.31780578284792793 0.11806620553186306 0.08373787105232085 0.06626139719751045 0.05077077240674423 0.03136773374212183 0.057533342926944606 0.08670532348850002; 0.004284498271424343 0.021591824367006424 0.02450278696113547 0.015605984149459226 0.01589725772402906 0.019503564697343635 0.014355232788254065 0.022436516938377876 0.02070669571070052 0.04238744809078793 0.14994212655588057 1.2077114272545073 0.22223937632714771 0.25904014132802244 3.8581893254118884 8.948819332322339 1.568392779795327 0.42948942876255314 0.13363035319237038 0.0894924745536024 0.07348968617414778 0.052915913801036635 0.03404122570477236 0.060734324524517075 0.08987902562105884; 0.004507485024430881 0.02128414694030165 0.02418773331811404 0.018764364687588853 0.01778834136180709 0.020918213376986017 0.015252622607948907 0.024287278105893213 0.020969793720703907 0.042858338638071106 0.14449657482037542 1.1497189285513227 0.2187944478663626 0.26030587064280164 3.7341374416275697 8.881773737814838 2.2158619322768884 0.541173074677178 0.14919450085287767 0.09524707805488393 0.08071797515078508 0.05506105519532904 0.03671471766742288 0.06393530612208953 0.09305272775361766; 0.004730471777437419 0.020976469513596874 0.023872679675092615 0.021922745225718476 0.019679424999585117 0.022332862056628402 0.016150012427643752 0.02613803927340855 0.02123289173070729 0.043329229185354284 0.13905102308487027 1.0917264298481382 0.21534951940557753 0.26157159995758084 3.6100855578432514 8.814728143307338 2.8633310847584497 0.6528567205918029 0.16475864851338495 0.10100168155616544 0.08794626412742237 0.05720619658962145 0.039388209630073404 0.067136287719662 0.09622642988617648; 0.004953458530443957 0.020668792086892103 0.023557626032071188 0.0250811257638481 0.02157050863736315 0.023747510736270787 0.017047402247338594 0.02798880044092389 0.021495989740710673 0.04380011973263747 0.13360547134936512 1.0337339311449538 0.21190459094479241 0.26283732927236 3.486033674058933 8.747682548799839 3.510800237240011 0.7645403665064278 0.18032279617389224 0.10675628505744697 0.09517455310405967 0.05935133798391386 0.042061701592723925 0.07033726931723445 0.0994001320187353; 0.005176445283450496 0.02036111466018733 0.02324257238904976 0.028239506301977732 0.023461592275141184 0.025162159415913175 0.01794479206703344 0.029839561608439234 0.02175908775071406 0.044271010279920656 0.12815991961385997 0.975741432441769 0.20845966248400732 0.2641030585871392 3.3619817902746147 8.68063695429234 4.158269389721575 0.876224012421053 0.19588694383439959 0.11251088855872853 0.10240284208069701 0.06149647937820628 0.04473519355537446 0.07353825091480692 0.10257383415129413; 0.0053994320364570345 0.020053437233482553 0.022927518746028334 0.03139788684010736 0.025352675912919214 0.02657680809555556 0.01884218188672828 0.03169032277595457 0.02202218576071744 0.044741900827203834 0.12271436787835481 0.9177489337385843 0.20501473402322226 0.26536878790191837 3.237929906490296 8.61359135978484 4.805738542203136 0.9879076583356778 0.21145109149490687 0.11826549206001005 0.10963113105733432 0.06364162077249869 0.047408685518024976 0.07673923251237938 0.10574753628385294; 0.005622418789463572 0.01974575980677778 0.022612465103006907 0.03455626737823698 0.027243759550697244 0.027991456775197945 0.019739571706423122 0.03354108394346991 0.022285283770720827 0.04521279137448701 0.11726881614284966 0.8597564350353999 0.20156980556243714 0.2666345172166975 3.1138780227059772 8.54654576527734 5.453207694684697 1.0995913042503027 0.22701523915541416 0.1240200955612916 0.11685942003397164 0.06578676216679108 0.050082177480675505 0.07994021410995185 0.10892123841641177; 0.00584540554247011 0.01943808238007301 0.02229741145998548 0.037714647916366605 0.02913484318847527 0.029406105454840327 0.020636961526117967 0.03539184511098525 0.02254838178072421 0.0456836819217702 0.11182326440734452 0.8017639363322153 0.19812487710165205 0.2679002465314767 2.989826138921659 8.479500170769839 6.100676847166259 1.2112749501649276 0.24257938681592148 0.1297746990625731 0.12408770901060892 0.0679319035610835 0.05275566944332602 0.0831411957075243 0.11209494054897058; 0.006068392295476649 0.019130404953368235 0.021982357816964053 0.04087302845449624 0.03102592682625331 0.030820754134482722 0.021534351345812813 0.037242606278500595 0.022811479790727597 0.046154572469053384 0.10637771267183936 0.7437714376290308 0.19467994864086696 0.2691659758462559 2.86577425513734 8.412454576262341 6.7481459996478215 1.322958596079553 0.25814353447642885 0.13552930256385465 0.13131599798724627 0.07007704495537591 0.05542916140597656 0.08634217730509679 0.11526864268152942; 0.006291379048483187 0.01882272752666346 0.021667304173942623 0.044031408992625864 0.03291701046403134 0.032235402814125104 0.022431741165507654 0.039093367446015936 0.02307457780073098 0.04662546301633656 0.1009321609363342 0.6857789389258461 0.19123502018008187 0.2704317051610351 2.7417223713530214 8.34540898175484 7.395615152129384 1.4346422419941778 0.27370768213693614 0.1412839060651362 0.13854428696388357 0.07222218634966832 0.05810265336862708 0.08954315890266924 0.11844234481408822; 0.005945081012206664 0.018728614232401696 0.021312636474666684 0.03969675887247996 0.03013447620092704 0.029931408080412503 0.021513683629800093 0.03686056769195569 0.022496886914371107 0.04497445040176341 0.10844225562728177 0.7683861149030353 0.19523977143692084 0.26229337853762197 2.7027644540659646 7.932870273176352 6.494573878963147 1.271925705507353 0.2515791870342675 0.13091647293808206 0.12877665842528058 0.06741193787578825 0.055905307232000365 0.08358857548193607 0.11235173262125246; 0.005616970312609012 0.018674835132401393 0.020956193159252665 0.03564332430908099 0.02746709515115821 0.027702435010922732 0.020654269060515474 0.03476257544674996 0.021886661394123187 0.04340914016455194 0.11658653432646833 0.8566088291489652 0.19993838789299778 0.2528813932325319 2.6424060977794643 7.497630336384534 5.651780983017238 1.1188362209713925 0.23014691851890634 0.1208606516865352 0.11921456336553957 0.06267549094089367 0.053623068977771984 0.07767744495037049 0.10650157434603617; 0.0053070469496902285 0.018661390226662543 0.02059797422770055 0.03187110530242895 0.024914867314724837 0.025548483605655772 0.019853497457653788 0.032799390710398724 0.021243901239987192 0.041929532304702075 0.12536499703389378 0.9504470816636346 0.20533086954831262 0.24219574924576462 2.560647302493519 7.039689171379381 4.867236464291655 0.975373788386296 0.20941087659085253 0.11111644231049544 0.10985800178466058 0.05801284554498452 0.0512559386059419 0.07180976730797248 0.1008918699884393; 0.005015310923450316 0.018688279515185154 0.020237979680010357 0.028380101852523824 0.02247779269162692 0.023469553864611623 0.019111368821215037 0.030971013482902 0.02056860645196315 0.040535626822213865 0.1347776437495583 1.049900872447045 0.21141721640286545 0.23023644657732034 2.4574880682081295 6.55904677816089 4.140940322786398 0.8415384077520633 0.189371061250106 0.10168384480996284 0.10070697368264354 0.053424001688060836 0.04880391611651015 0.06598554255474204 0.09552261954846188; 0.004741762233889273 0.01875550299796922 0.01987620951618208 0.025170313959365642 0.020155871281864468 0.0214656457877903 0.018427883151199226 0.029277443764259788 0.019860777030051048 0.0392274237170873 0.1448244744734617 1.1549702014991954 0.2181974284566562 0.21700348522719903 2.3329283949232966 6.055703156729068 3.472892558501471 0.7173300790686956 0.170027472496667 0.09256285918493745 0.09176147905948855 0.04890895937012262 0.046267001509476734 0.0602047706906792 0.0903938230261039; 0.0044864008810071 0.018863060675014743 0.019512663736215713 0.022241741622954392 0.01794910308543748 0.019536759375191795 0.017803040447606354 0.02771868155447208 0.019120412974250885 0.03800492298932238 0.15550548920560409 1.2656550688200858 0.22567150570968486 0.2024968651954006 2.1869682826390195 5.529658307083912 2.863093171436871 0.6027488023361921 0.15138011033053536 0.08375348543541924 0.08302151791519556 0.04446771859116987 0.04364519478484164 0.054467451715783946 0.08550548042136535; 0.004249226864803798 0.01901095254632172 0.019147342340111263 0.019594384843290073 0.015857488102345954 0.01768289462681611 0.017236840710436414 0.026294726853538883 0.018347514284562667 0.036868124638919086 0.16682068794598542 1.381955474409716 0.23383944816195149 0.18671658648192518 2.0196077313552983 4.980912229225421 2.311542161592598 0.4977945775545527 0.13342897475171112 0.0752557235614082 0.07448709024976458 0.04010027935120258 0.04093849594260487 0.048773585630056265 0.08085759173424628; 0.004030240185279365 0.019199178611890154 0.018780245327868726 0.017228243620372674 0.01388102633258988 0.015904051542663236 0.016729283939689414 0.025005579661460185 0.017542080960986386 0.03581702866587744 0.17877007069460574 1.5038714182680875 0.24270125581345603 0.16966264908677256 1.8308467410721319 4.4094649231535925 1.818239528968651 0.402467404723777 0.11617406576019418 0.0670695735629043 0.06615819606319553 0.035806641650220725 0.0381469049827664 0.04312317243349614 0.07645015696474661; 0.003829440842433803 0.019427738871720054 0.018411372699488112 0.015143317954202217 0.012019717776169279 0.014200230122733193 0.016280370135365357 0.02385123997823601 0.01670411300352205 0.034851635070197444 0.19135363745146505 1.6314029003951989 0.2522569286641986 0.151335053009943 1.6206853117895224 3.8153163888684336 1.383185273565033 0.3167672838438662 0.09961538335598473 0.0591950354399076 0.058034835355488555 0.03158680548822436 0.035270421905326275 0.03751621212610364 0.07228317611286642; 0.00364682883626711 0.0196966333258114 0.018040724454969403 0.013339607844778686 0.010273562433084136 0.012571430367025963 0.015890099297464228 0.022831707803866324 0.015833610412169654 0.03397194385187907 0.2045713882165632 1.7645499207910493 0.26250646671417893 0.13173379825143633 1.3891234435074682 3.19846662636994 1.0063793953817421 0.24069421491481943 0.08375292753908264 0.05163210919241808 0.05011700812664356 0.027440770865213446 0.032309046710284464 0.03195270470787871 0.06835664917860562; 0.003482404166779289 0.020005861974164207 0.017668300594312618 0.011817113292102093 0.008642560303334457 0.011017652275541556 0.015558471425986046 0.021946983138351164 0.014930573186929204 0.03317795501092236 0.21842332298990041 1.903312479455641 0.2734498699633974 0.11085888481125263 1.1361611362259707 2.5589156356581126 0.6878218944187784 0.17424819793663693 0.06858669830948796 0.04438079482043574 0.04240471437666059 0.023368537781188 0.029262779397640992 0.026432650178821376 0.06467057616196431; 0.003336166833970336 0.02035542481677846 0.017294101117517738 0.010575834296172424 0.007126711386920239 0.00953889584827996 0.015285486520930792 0.0211970659816905 0.01399500132780069 0.032469668547327274 0.2329094417714765 2.047690576388972 0.28508713841185357 0.08871031268939183 0.8617983899450286 1.8966634167329506 0.42751277067614196 0.11742923290931864 0.05411669566720066 0.03744109232396056 0.03489795410553959 0.019370106236148012 0.02613161996739583 0.020956048538931606 0.06122495706294241; 0.0032081168378402546 0.02074532185365419 0.01691812602458478 0.009615770856989688 0.005726015683841483 0.008135161085241191 0.015071144582298482 0.02058195633388435 0.013026894834784123 0.03184708446109385 0.2480297445612916 2.1976842115910444 0.2974182720595479 0.065288081885854 0.5660352046646426 1.2117099695944544 0.22545202415383264 0.07023731983286455 0.04034291961222077 0.03081300170299257 0.027596727313280618 0.015445476230093489 0.022915568419549 0.015522899788209431 0.05801979188153998; 0.0030982541783890418 0.02117555308479136 0.016540375315513736 0.008936922974553878 0.004440473194098185 0.006806447986425233 0.014915445610089104 0.020101654194932702 0.012026253707879492 0.03131020275222205 0.2637842313593456 2.353293385061856 0.31044327090647994 0.0405921924006391 0.24887158038481239 0.5040552942426242 0.08163965485185048 0.03267245870727467 0.027265370144548252 0.024496522957531738 0.020501033999883626 0.011594647763024421 0.019614624754100488 0.01013320392665483 0.05505508061775698; 0.003986060199204506 0.021667187117796857 0.01636513474953283 0.008961312317088707 0.003984787990792592 0.006593319207341276 0.01568255954966047 0.026585774437496455 0.09323269334405329 0.09894405893660566 0.2720159537942247 2.405810542249199 0.3132574140559028 0.023802033162204958 0.025510611235400365 0.019857044300259814 0.018342042545192346 0.013150324561259558 0.018986928009989176 0.02074879622283559 0.016003939766600226 0.009483012447212029 0.01737059247264471 0.007269691016504726 0.05240627180170216; 0.013372182982866916 0.022155568535481448 0.01656164409942635 0.00986570251646603 0.004944644709510053 0.008389942979407048 0.01805866197172783 0.04581474041317424 0.3275943486229148 0.29335635548569283 0.2640765795721831 2.2464238510404266 0.29485531428654443 0.024856788637303444 0.02538637723863029 0.02035090249159238 0.018934962928222495 0.013624389939802934 0.019109317598796076 0.02154713402505953 0.016320151162677467 0.010655363193894584 0.017354130581092475 0.009278526594886734 0.05001626806215685; 0.032558243211696786 0.022555698110233387 0.01670380885351584 0.010683421218062494 0.0058261254852875945 0.0100419967795584 0.020226415798254266 0.06360472661611265 0.544986248491847 0.4736383548454688 0.25603431347488276 2.0924962817597463 0.2769960802296381 0.025777578738165316 0.025208183080951168 0.020760576121720744 0.019440354410629333 0.014032203543619663 0.019176199241579848 0.022237342289137325 0.016574357459609775 0.01172061742801308 0.017296059828413073 0.01112596827477222 0.047667899388239905; 0.06154424088569435 0.02286757584205268 0.0167916290118013 0.011414468421878102 0.006629230318125225 0.011549480607795348 0.022185821029239796 0.0799557330463119 0.7454083929508517 0.639790057015935 0.24788915550232368 1.9440278344071564 0.25967971188518363 0.02656440346479059 0.02497602876236301 0.02108606519064492 0.019858216992412844 0.014373765372709759 0.019187572938340496 0.022819421015068994 0.01676655865739716 0.012678775149567532 0.017196380214606517 0.012812016056161203 0.045361165779951286; 0.10033017600485925 0.023091201730939315 0.016825104574282733 0.012058844127912846 0.0073539592080229315 0.012912394464117873 0.023936877664684393 0.09486775970377176 0.9288607819999266 0.7918114619970896 0.23964110565450605 1.8010185089826598 0.2429062092531812 0.02721726281717924 0.024689914282865807 0.021327369698364906 0.020188550673573034 0.014649075427073206 0.019143438689078005 0.023293370202854526 0.01689675475603962 0.013529836358557923 0.017055091739672804 0.014336669939053664 0.04309606723729104; 0.1489160485691916 0.023226575776893302 0.01680423554096014 0.012616548336166732 0.008000312154980722 0.014130738348525986 0.025479585704588074 0.10834080658849235 1.095343415639073 0.9297025697889336 0.23129016393142968 1.663468305486254 0.22667557233363075 0.027736156795331285 0.024349839642459563 0.021484489644880696 0.020431355454109904 0.014858133706710012 0.019043796493792398 0.02365918985249392 0.01696494575553715 0.014273801054984263 0.01687219440361193 0.01569992992344961 0.04087260376025915; 0.20730185857869146 0.023273697979914635 0.016729021911833515 0.013087581046639757 0.008568289158998594 0.015204512261019681 0.026813945148950833 0.12037487370047362 1.2448562938682906 1.053463380391467 0.2228363303330947 1.5313772239179404 0.21098780112653223 0.028121085399246716 0.023955804841144273 0.021557425030192295 0.020586631334023452 0.015000940211620177 0.018888646352483662 0.023916879963987184 0.016971131655889753 0.014910669238846548 0.016647688206423896 0.01690179600934904 0.0386907753488556; 0.27548760603335876 0.02323256834000332 0.016599463686902863 0.013471942259331925 0.00905789022007655 0.016133716201598966 0.02793995599777268 0.13096996103971564 1.3773994166875794 1.1630938938046897 0.21427960485950098 1.404745264277718 0.1958428956318856 0.02837204862892554 0.023507809878919945 0.0215461758542997 0.02065437831331368 0.015077494941803701 0.018677988265151796 0.024066440537334314 0.01691531245709743 0.015440440910144782 0.016381573148108702 0.017942268196751963 0.03655058200308041; 0.3534732909331936 0.023103186857159355 0.016415560866168185 0.013769631974243236 0.00946911533821459 0.016918350170263836 0.028857618251053603 0.14012606860621835 1.4929727840969398 1.258594110028602 0.20561998751064864 1.2835724265655875 0.18124085584969096 0.028489046484367758 0.023005854755786582 0.02145074211720292 0.02063459639198059 0.015087797897260585 0.0184118222317968 0.024107871572535307 0.016797488159160184 0.01586311606887896 0.01607384922866635 0.01882134648565837 0.03445202372293357; 0.4412589132781957 0.022885553531382732 0.016177313449629478 0.013980650191373683 0.009801964513412713 0.017558414167014293 0.029566931908793604 0.14784319639998178 1.5915763960963716 1.3399640290632033 0.19685747828653763 1.1678587107815483 0.16718168177994827 0.028472078965573355 0.022449939471744174 0.021271123818901942 0.020527285570024176 0.015031849077990826 0.018090148252418677 0.024041173069590166 0.016617658762078012 0.01617869471504909 0.015724516448096842 0.019539030876068265 0.03239510050841508; 0.5388444730683659 0.02257966836267346 0.015884721437286742 0.014104996910723276 0.010056437745670919 0.01805390819185034 0.030067896970992693 0.15412134442100595 1.6732102526858752 1.4072036509084944 0.1879920771871679 1.0576041169256003 0.1536653734226575 0.028321146072542356 0.021840064026792713 0.021007320959396773 0.020332445847444442 0.014909648483994422 0.017712966327017427 0.02386634502849889 0.016375824265850908 0.016387176848655163 0.015333574806400168 0.020095321367981644 0.030379812359524937; 0.6462299703037032 0.022185531351031542 0.015537784829139983 0.014142672132292006 0.010232535034989207 0.01840483224477197 0.03036051343765086 0.15896051266929084 1.7378743538654502 1.4603129755644748 0.1790237842125395 0.9528086449977444 0.14069193077781872 0.028036247805274743 0.021176228420932218 0.020659333538687413 0.020050077224241385 0.014721196115271385 0.01728027645559305 0.023583387449261486 0.01607198467047888 0.016488562469697182 0.014901024303576338 0.020490217961398507 0.028406159276263145; 0.7634154049842079 0.02170314249645697 0.015136503625189195 0.014093675856079875 0.010330256381367578 0.018611186325779183 0.030444781308768105 0.16236070114483642 1.7855686996350961 1.4992920030311443 0.1699525993626525 0.8534722949979804 0.12826135384543189 0.02761738416377052 0.02045843265416268 0.02022716155677386 0.01968017970041501 0.014466491971821701 0.016792078638145545 0.02319230033187794 0.01570613997596193 0.016482851578175153 0.01442686493962535 0.02072372065631886 0.026474141258629708; 0.8904007771098802 0.021132501798949745 0.014680877825434376 0.01395800808208689 0.010349601784806034 0.01867297043487199 0.03032070058434444 0.16432190984764272 1.8162932899948134 1.5241407333085033 0.16077852263750675 0.7595950669263075 0.11637364262549699 0.027064555148029683 0.019686676726484102 0.019710805013656112 0.019222753275965315 0.014145536053645377 0.016248372874674916 0.02269308367634826 0.01527829018230005 0.016370044174089063 0.013911096714547204 0.020795829452742706 0.024583758306624626; 1.0271860866807196 0.02047360925850987 0.014170907429875533 0.013735668810313045 0.01029057124530457 0.01859018457205038 0.029988271264379848 0.1648441387777098 1.8300481249446021 1.5348591663965518 0.1515015540371024 0.6711769607827264 0.10502879711801405 0.02637776075805224 0.018860960637896485 0.019110263909334176 0.0186777979508923 0.01375832836074241 0.015649159165181153 0.022085737482672452 0.014788435289493248 0.016150140257438925 0.013353719628341896 0.020706544350670028 0.022735010420247895; 1.1737713336967275 0.019726464875137347 0.013606592438512661 0.01342665804075834 0.010153164762863192 0.018362828737314357 0.02944749334887434 0.16392738793503753 1.8268332044844626 1.5314473022952897 0.14212169356143933 0.5882179765672367 0.09422681732298305 0.02555700099383819 0.017981284388399823 0.01842553824380805 0.01804531372519596 0.013304868893112803 0.014994437509664268 0.02137026175085051 0.014236575297541517 0.01582313982822474 0.012754733681009432 0.020455865350100848 0.020927897599499513; 1.3301565181579027 0.018891068648832166 0.012987932851345755 0.013030975773422774 0.009937382337481895 0.017990902930663917 0.028698366837827907 0.16157165731962597 1.806648528614394 1.5139051410047166 0.13263894121051756 0.5107181142798384 0.08396770324040395 0.024602275855387518 0.017047647977994115 0.017656628017077722 0.017325300598876294 0.012785157650756551 0.01428420790812425 0.02054665648088242 0.013622710206444854 0.015389042886446491 0.012114138872549803 0.020043792451035138 0.019162419844379473; 1.4963416400642449 0.017967420579594336 0.012314928668374827 0.012548622008306351 0.00964322396916068 0.017474407152099067 0.02774089173124056 0.15777694693147512 1.7694940973343969 1.482232682524833 0.12305329698433715 0.43867737392053197 0.07425145487027687 0.023513585342700247 0.01606005140667937 0.016803533229143207 0.016517758571933314 0.01219919463367366 0.013518470360561105 0.01961492167276821 0.012946840016203269 0.014847849432104196 0.011431935202963019 0.019470325653472925 0.017438577154887794; 1.6723266994157542 0.016955520667423854 0.011587579889599872 0.011979596745409066 0.00927068965789955 0.0168133414016198 0.026575068029112296 0.152543256770585 1.7153699106444709 1.4364299268556389 0.11336476088289807 0.3720957554893173 0.06507807221260173 0.022290929455776362 0.015018494674455582 0.015866253880004504 0.015622687644367016 0.01154697984186413 0.012697224866974834 0.01857505732650786 0.012208964726816759 0.014199559465197843 0.010708122672249078 0.018735464957414197 0.015756369531024467; 1.8581116962124316 0.015855368912320728 0.010805886515020888 0.011323899984730925 0.008819779403698503 0.016007705679226124 0.025200895731443115 0.14587058683695564 1.6442759685446169 1.3764968739971342 0.10357333290620035 0.3109732589861941 0.05644755526737853 0.02093430819461587 0.013922977781322754 0.014844789969661605 0.014640087816177393 0.010828513275327956 0.011820471427365434 0.017427063442101383 0.011409084338285322 0.013444172985727441 0.009942701280407975 0.01783921036285896 0.014115796972789495; 2.0536966304542763 0.01466696531428494 0.009969848544637876 0.010581531726271926 0.008290493206557538 0.015057499984918031 0.023618374838233012 0.13775893713058696 1.556212271034834 1.3024335239493188 0.09367901305424393 0.2553098844111626 0.048359904034607305 0.019443721559218766 0.012773500727280883 0.013739141498114517 0.013569959087364454 0.010043794934065143 0.01088821004173291 0.016170940019548767 0.010547198850608959 0.012581689993692986 0.009135671027439715 0.016781561869807206 0.012516859480182869; 2.259081502141289 0.0133903098733165 0.009079465978450833 0.00975249197003206 0.007682831066476654 0.01396272431869552 0.021827505349481982 0.12820830765147892 1.451178818115122 1.2142398767121922 0.08368180132702879 0.20510563176422236 0.04081511851428798 0.017819169549585047 0.011570063512329965 0.01254930846536323 0.012412301457928184 0.009192824818075682 0.00990044071007725 0.014806687058850007 0.009623308263787666 0.011612110489094472 0.00828703191334429 0.015562519478258931 0.010959557053204588; 2.4742663112734693 0.012025402589415417 0.008134738816459767 0.008836780716011343 0.006996792983455857 0.012723378680558603 0.019828287265190044 0.11721869839963167 1.3291756097854819 1.1119159322857555 0.07358169772455506 0.16036050104537403 0.03381319870642065 0.01606065216571473 0.010312666136470014 0.011275290871407757 0.011167114927868602 0.008275602927359585 0.008857163432398469 0.013334304560005126 0.008637412577821451 0.01053543447193191 0.007396783938121711 0.014182083188214152 0.009443889691854668; 2.6992510578508155 0.01057224346258168 0.007135667058664673 0.007834397964209766 0.006232378957495141 0.011339463070507271 0.01762072058535718 0.10479010937504513 1.190202646045913 0.9954616906700081 0.06337870224682263 0.12107449225461722 0.02735414461100525 0.014168169407607796 0.009001308599701018 0.00991708871624809 0.009834399497185694 0.007292129261916846 0.00775837820869656 0.011753792523014103 0.007589511792710309 0.009351661942205294 0.006464927101771972 0.012640252999672856 0.007969857396133099; 2.934035741873331 0.009030832492815293 0.006082250705065552 0.00674534371462733 0.005389588988594509 0.009810977488541525 0.015204805309983405 0.09092254057771933 1.0342599268964159 0.8648771518649502 0.05307281489383153 0.08724760539195207 0.021437956228041816 0.012141721275264254 0.007635990902022981 0.008474701999884234 0.008414155165879471 0.0062424038217474675 0.006604085038971525 0.010065150947876954 0.006479605908454241 0.008060792899914626 0.005491461404295074 0.010937028912635052 0.006537460166039883; 3.178620363341012 0.007401169680116254 0.004974489755662401 0.005569617967264034 0.00446842307675396 0.008137921934661365 0.012580541439068706 0.07561599200765423 0.8613474523369898 0.7201623158705817 0.042664035665581766 0.05887984045737848 0.016064633557530322 0.009981307768684104 0.006216713043435901 0.006948130722316183 0.006906381933949924 0.005126426606851446 0.005394283923223359 0.008268379834593666 0.005307694925053246 0.006662827345059906 0.004476386845691018 0.009072410927100729 0.005146698001575016; 3.4330049222538626 0.005683255024484565 0.003812384210455224 0.00430722072211988 0.0034688812219734943 0.006320296408866792 0.009747928972613094 0.058870463664849854 0.6714652223676354 0.5613171826869026 0.032152364562073335 0.03597119745089648 0.011234176599470784 0.007686928887867345 0.004743475023939781 0.005337374883543943 0.005311079801397059 0.003944197617228784 0.004128974861452068 0.006363479183164246 0.004073778842507326 0.005157765277641132 0.003419703425959802 0.007046399043069897 0.0037975709027385013; 3.69718941861188 0.0038770885259202153 0.0025959340694440124 0.0029581519791948593 0.002390963424253105 0.004358100911157798 0.006706967910616544 0.040685955549306096 0.4646132369883512 0.3883417523139119 0.02153780158330617 0.018521676372506 0.006946585353863173 0.005258584632813963 0.0032162768435346108 0.0036424344835675002 0.0036282487682208625 0.0026957168528794727 0.0028081578536576415 0.00435044899358868 0.0027778576608164726 0.003545606697658297 0.0023214111451014216 0.004858993260542538 0.0024900788695303316; 3.9711738524150646 0.001982670184423223 0.0013251393326287788 0.001522411738488986 0.0012346696835928055 0.0022513354415343977 0.0034576582530790904 0.02106246766102314 0.24079149619913948 0.20123602475161143 0.010820346729280395 0.006531277222207204 0.003201859820707534 0.002696275003523983 0.0016351185022204064 0.001863309522386875 0.0018578888344213539 0.001380984313803527 0.0014318328998400942 0.0022292892658669914 0.0014199313799806994 0.0018263516051114168 0.0011815100031158875 0.002510193579518678 0.00122422190195052; 4.254958223663417 -6.420612698722547e-15 9.517855303472442e-15 2.254439019145773e-15 -7.411655927494534e-15 -3.415475412410896e-15 7.193291213800562e-16 9.03543153937406e-16 -9.219208521726188e-16 3.3188559986441705e-16 -4.05636043338745e-15 0.0 3.844317340337979e-15 -2.6066681841577996e-15 -2.839447505217288e-15 2.056986106201179e-15 -1.475037653587726e-15 9.403365037570114e-16 -5.807960757836408e-16 -8.320134259982439e-16 0.0 4.839967298001723e-16 3.1943784166793436e-15 -1.6962933006317804e-15 -9.403365036110586e-16]