Proposed modifications to the interface

[Lazersmoke]

Extremely open to changes on all of this! This is just to get the ball rolling so we can start landing breaking changes :)

The main changes are:

Changes to the intensities interface:

• intensities renamed to intensities_interpolated, and option interpolation = :none renamed to interpolation = :round
• Options kT, formfactors, and contractor are removed from intensities and intensities_binned
• New function intensity_formula(sf, :perp; kT, formfactors) (documentation reproduced below)
• Support for custom intensity formulae via do-notation
• Same changes applied to powder averaging functions (they now accept intensity formulas)

Supporting changes to the internals:

• dipole_trajectory and expectation_trajectory have been combined
• Everything is now a custom observable (subject to benchmarking), including the dipole component observables via LinearMaps (new dependency). [N.B. the custom dipole operators is a Martin requested feature in disguise]
• Observables are now named by Symbols in StructureFactor in order to facilitate users referring to them (e.g. when specifying correlations to save) without having to use indices all the time. Providing observable names is optional (default alphabetical) and indexing by number is still supported.
• show(::StructureFactor) facelift:

StructureFactor (37.348 MiB)
[S(q,ω) | nω = 239 | 1 sample]
6 correlations on (16, 16, 4) lattice:
╔ ⬤ ⬤ ⬤ Sx
║ ⋅ ⬤ ⬤ Sy
╚ ⋅ ⋅ ⬤ Sz

Sparsity pattern looks better in reality:

Example creating a SU(2) system with only Sx and Sz

S = Sunny.spin_matrices(2)
observables = [:Sx => S[1], ], :Sz => S[3]]
sf = StructureFactor(...; observables)

Documentation for intensity_formula

(reproduced here for clarity)

formula = intensity_formula(sf::StructureFactor, contraction_mode; kwargs...)
formula(sf,k,cell,iω)

Establish a formula for computing 𝒮(𝐪,ω) from the correlation data stored in the StructureFactor. The formula returned from intensity_formula can be passed to intensities_interpolated or intensities_binned.

Sunny has several built-in formulas that can be selected by setting contraction_mode to one of these values:

• :trace, which yields tr⁡𝒮(q,ω)=∑α𝒮αα(q,ω)
• :perp, which contracts 𝒮αβ(q,ω) with the dipole factor δαβ−qαqβ​, returning the unpolarized intensity.
• :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 FormFactors, 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 specified 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 q, ω, 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.

Modified keyword arguments to StructureFactor

• observables: Allows the user to specify custom observables. The observables must be given as a list of complex N×N matrices or LinearMaps. 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 Sαβ(q,ω) are calculated and stored. Specified with a vector of tuples. By default Sunny records all auto- and 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)].

[ddahlbom]

Thanks for this! Most of this looks great. I'll to go over this with a bit more care tomorrow and leave comments in the code where relevant.

Some initial thoughts:

1. The updated SF display looks great.
2. I like the optional use of symbols for specifying correlations.
3. Didn't know about LinearMaps. This could actually be useful for a number of other things that are going on (e.g., for KPM, Lanczos, etc.)
4. I like the approach to intensity_formula from a programming standpoint. I think maybe we'll want a different name. I was slightly put off by the notion that a user would have to configure a function to do a standard calculation, but it looks like you have the basic use cases covered such that the user may not even realize they're passing a function to intensities, just some configuration thing. In general, we've been thinking of the typical user as someone who knows some Matlab or Python, is not primarily a programmer, and probably is just learning Julia to get their calculation done. But we do want everything configurable for the advanced user, and this is a nice way to do it.
5. It would be good to hear more about your experiences benchmarking some of these changes.

[Lazersmoke]

Re: 4.

FYI: it defaults to formula = intensity_formula(sf,:perp), so for the "standard" calculation nothing needs to be specified.

Agreed that the user shouldn't need to think about the fact that they're passing a function to intensities, unless it's their custom function. "Formula" is meant to imply mathematical formula: specifically a formula for the scattering cross section in terms of the structure factor, e.g.:

Many such formulae exist, none more correct than the other so I wouldn't call it "corrections"

It may be better to introduce a struct for this instead of a function. intensity_formula would become a constructor IntensityFormula (or CrossSectionFormula mayhaps?). Benefits include:

• It could be pretty-printed, perhaps even as a mathematical formula (albeit typeset in unicode lol)
• The return_type mechanism would be way less hacky and could be a type parameter

[kbarros]

You both have more experience than me, so I'll be happy to try out whatever you decide. Given how many moving parts there are, it might be good to push a prototype to main and collect feedback. What do you think -- can we label main as the 0.5.0 development branch and start making breaking changes on it?

[ddahlbom]

I'm definitely for merging as soon as possible. I think maybe we should tweak some of the naming, but maybe that's best done in conjunction with users.

Two things I think should be addressed beforehand.

1. I do think a better approach is to have a struct for the correction information. kT and formfactors will be handled in different ways with the SWT/KPM material, so the same function couldn't be reused in both contexts. I think specifying what is now called a contractor as a function makes a lot of sense, and that could be reused in exactly the same way in both parts of the code. So, maybe the solution is a struct with contractor a function kept in the struct along with kT and formfactors?

2. I think benchmarking should be done to the level that we're confident this approach doesn't present any fundamental obstacles to achieving similar or better performance to the existing situation. Part of the goal of the original intensities function was basically to create plots quickly, maybe interactively.

[Lazersmoke]

Re: Benchmarking

Firstly, the only possible slow-down is with the changes to the internals; intensities should see no performance impact.

My holistic benchmark for this is just to time how long it takes to perform one "forwards problem" for the FeI2 example. Without this PR, it takes 59-60 seconds, and after this PR it takes up to 60.6 seconds. So whatever performance regression there may be with this doesn't impact things enough to worry about in my opinion.

[Lazersmoke]

kT and formfactors will be handled in different ways with the SWT/KPM material, so the same function couldn't be reused in both contexts

Even worse, the shape of the calculation of intensities is not the same for classical (SF) vs quantum (SWT) calculations. I'm working on making ClassicalIntensityFormula and SpinWaveIntensityFormula to address this; we'll see if the result is likeable.

[Lazersmoke]

[Lazersmoke]

There's no docs for intensity_formula(::SpinWaveTheory) at the moment, but also it doesn't do much because there aren't form factors/kT/custom observables for SWT yet. The Contraction code is reused though!

If the docs build, this is OK to be a prototype on main @kbarros

[ddahlbom]

This looks great to me. I think it would be good to get into main once tests and docs build.

[ddahlbom]

Maybe I'm not following. If you're doing dipole correlations in SU(N) mode, it will take the first branch (since sf.dipole_corrs will be true). (The checking for 3 is not necessary, however, when sf.dipole_corrs is true.)

I guess there is ambiguity here in terms of what is meant by "trace." Basically, "trace" means Sxx+Syy+Szz unless the user has specified custom observables/correlations, in which case it means \sum Taa where a runs from 1 to N^2-1.

My feeling was that there should be some check here to make sure your not just summing some random collection of diagonal elements. But I guess we expect a user specifying arbitrary correlations to be sophisticated enough to know what they're doing.

[Lazersmoke]

If you're in SU(N) mode, you may still have a number different from N*N-1 number of correlations included in your trace, for example if you have M observables, you have M =/= N*N-1 terms.

The code now just directly checks that, for each observable in the structure factor, that observable's autocorrelation is included in the trace.

[ddahlbom]

Indeed, the point of the check was precisely to prevent the user from thinking they're getting a trace if M is not equal to the dimension of the Lie algebra su(N), namely N^2-1. This has physical significance. (The check is quick and dirty in the sense that it doesn't check is the N^2-1 operators indeed form an orthonormal basis for NxN Hermitian matrices.)

I'm not sure this should be called Trace if it's just a sum or some available set of autocorrelations, and I do think there should be a trace option. Maybe Trace should just mean Sxx+Syy+Szz, which is what most people will expect. If someone wants to do a trace of a full set of SU(N) observables, they'll just have to write their own formula.

[ddahlbom]

The code now just directly checks that, for each observable in the structure factor, that observable's autocorrelation is included in the trace.

I guess this is fine too, now that I read it a bit more carefully, since, in the default situation, the user will get Sxx+Syy+Szz. In any case, just explaining the reasoning behind the check.