Lswt general observables

examples/fei2_classical.jl

@@ -70,6 +70,7 @@ set_exchange!(sys, [J′2apm 0.0 0.0; 0.0 J′2apm 0.0; 0.0 0.0 J′2azz], Bond(
 D = 2.165#hide
 S = spin_operators(sys, 1)#hide
 set_onsite_coupling!(sys, -D*S[3]^2, 1)#hide
+sys
 
 # ## Finding a ground state
 

examples/fei2_tutorial.jl

@@ -162,7 +162,8 @@ set_exchange!(sys, [J′2apm 0.0    0.0;
 
 D = 2.165
 S = spin_operators(sys, 1)
-set_onsite_coupling!(sys, -D*S[3]^2, 1)
+set_onsite_coupling!(sys, -D*S[3]^2, 1);
+sys#hide
 
 # Any anisotropy operator can be converted to a linear combination of Stevens
 # operators with [`print_stevens_expansion`](@ref).
@@ -233,6 +234,9 @@ suggest_magnetic_supercell([[0, -1/4, 1/4]], sys.latsize)
 # ground-state, which makes optimization much easier.
 
 sys_min = reshape_supercell(sys, [1 0 0; 0 1 -2; 0 1 2])
+
+#
+
 randomize_spins!(sys_min)
 minimize_energy!(sys_min)
 plot_spins(sys_min; ghost_radius=3)

examples/out_of_equilibrium.jl

@@ -64,6 +64,7 @@ field = set_external_field!(sys, [0.0 0.0 h]);
 D = 19.0
 Sz = Sunny.spin_operators(sys, 1)[3]
 set_onsite_coupling!(sys, D*Sz^2, 1);
+sys#hide
 
 # Initialize system to an infinite temperature (fully randomized) initial
 # condition.
@@ -131,4 +132,4 @@ plot_triangular_plaquettes(sun_berry_curvature, frames; resolution=(1800,600),
 # that the process has generated a number of well-formed skyrmions of both
 # positive (red) and negative (blue) charge in addition to a number of other
 # metastable spin configurations. A full-sized version of this figure is
-# available in [Dahlbom et al.](https://doi.org/10.1103/PhysRevB.106.235154).
+# available in [Dahlbom et al.](https://doi.org/10.1103/PhysRevB.106.235154).

src/Intensities/ElementContraction.jl

@@ -15,21 +15,21 @@ struct Element <: Contraction{ComplexF64}
     index :: Int64
 end
 
-struct FullTensor{NCorr} <: Contraction{SVector{NCorr, ComplexF64}} end
+struct FullTensor{NCorr,NSquare,NObs} <: Contraction{SMatrix{NObs, NObs, ComplexF64}}
+    indices :: SVector{NSquare, Int64}
+end
 
 
 ################################################################################
 # Constructors
 ################################################################################
-Trace(swt::SpinWaveTheory) = Trace(@SVector[1,5,9])
-
-function Trace(sc::SampledCorrelations{N}) where {N}
+function Trace(obs::ObservableInfo)
     # Collect all indices for matrix elements 𝒮^αβ where α=β
     indices = Int64[]
-    for (ki,i) = sc.observable_ixs
+    for (ki,i) = obs.observable_ixs
         autocorrelation_index = CartesianIndex(i,i)
-        if haskey(sc.correlations,autocorrelation_index)
-            push!(indices,sc.correlations[autocorrelation_index])
+        if haskey(obs.correlations,autocorrelation_index)
+            push!(indices,obs.correlations[autocorrelation_index])
         else
             problematic_correlation = ki
             error("Can't calculate trace because auto-correlation of the $problematic_correlation observable was not computed.")
@@ -52,31 +52,33 @@ function Trace(sc::SampledCorrelations{N}) where {N}
     Trace(SVector{length(indices), Int64}(indices))
 end
 
-DipoleFactor(swt::SpinWaveTheory) = DipoleFactor([1,4,5,7,8,9])
-
-function DipoleFactor(sc::SampledCorrelations{N}; spin_components = [:Sx,:Sy,:Sz]) where {N}
+function DipoleFactor(obs::ObservableInfo; spin_components = [:Sx,:Sy,:Sz])
     # Ensure that the observables themselves are present
     for si in spin_components
-        if !haskey(sc.observable_ixs,si)
+        if !haskey(obs.observable_ixs,si)
             error("Observable $(si) missing, but required for dipole correction factor")
         end
     end
 
     # 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(sc,dipole_correlations; err_msg = αβ -> "Missing correlation $(αβ), which is required to compute the depolarization correction.")
+    indices = lookup_correlations(obs,dipole_correlations; err_msg = αβ -> "Missing correlation $(αβ), which is required to compute the depolarization correction.")
     DipoleFactor(indices)
 end
 
-function Element(sc::SampledCorrelations, pair::Tuple{Symbol,Symbol})
-    Element(only(lookup_correlations(sc,[pair]; err_msg = pair -> "Missing correlation $(pair), which was requested.")))
+function Element(obs::ObservableInfo, pair::Tuple{Symbol,Symbol})
+    Element(only(lookup_correlations(obs,[pair]; err_msg = pair -> "Missing correlation $(pair), which was requested.")))
 end
 
-FullTensor(swt::SpinWaveTheory) = FullTensor{9}()
-
-function FullTensor(sc::SampledCorrelations{N}) where {N}
-    FullTensor{size(sc.data, 1)}()
+function FullTensor(obs::ObservableInfo)
+    n_obs = num_observables(obs)
+    tensor_elements = Matrix{Tuple{Symbol,Symbol}}(undef,n_obs,n_obs)
+    for (ki,i) = obs.observable_ixs, (kj,j) = obs.observable_ixs
+      tensor_elements[i,j] = (ki,kj) # Required to put matrix in correct order
+    end
+    indices = lookup_correlations(obs, collect(tensor_elements); err_msg = αβ -> "Missing correlation $(αβ). All correlations are required to return the full tensor.")
+    FullTensor{num_correlations(obs),length(indices),num_observables(obs)}(indices)
 end
 
 ################################################################################
@@ -113,7 +115,11 @@ end
 
 contract(specific_element, _, ::Element) = only(specific_element)
 
-contract(all_elems, _, ::FullTensor) = all_elems
+function contract(all_elems, _, full::FullTensor{NCorr,NSquare,NObs}) where {NCorr, NSquare,NObs}
+    # This Hermitian takes only the upper triangular part of its argument
+    # and ensures that Sαβ has exactly the correct symmetry
+    Hermitian(reshape(all_elems[full.indices],NObs,NObs))
+end
 
 ################################################################################
 # Contraction utils
@@ -131,13 +137,13 @@ Base.zeros(::Contraction{T}, dims...) where T = zeros(T, dims...)
 
 function contractor_from_mode(source, mode::Symbol)
     if mode == :trace
-        contractor = Trace(source)
+        contractor = Trace(source.observables)
         string_formula = "Tr S"
     elseif mode == :perp
-        contractor = DipoleFactor(source)
+        contractor = DipoleFactor(source.observables)
         string_formula = "∑_ij (I - Q⊗Q){i,j} S{i,j}\n\n(i,j = Sx,Sy,Sz)"
     elseif mode == :full
-        contractor = FullTensor(source)
+        contractor = FullTensor(source.observables)
         string_formula = "S{α,β}"
     end
     return contractor, string_formula
@@ -184,4 +190,4 @@ function error_formula(sc::SampledCorrelations, mode::Symbol; kwargs...)
         error("Error information not available for this `SampledCorrelation`.")
     end
     intensity_formula(sc, mode; calculate_errors=true, kwargs...)
-end
+end

src/Operators/Observables.jl

@@ -0,0 +1,141 @@
+
+struct ObservableInfo
+    # Correlation info (αβ indices of S^{αβ}(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 S^{αβ}(q, ω)
+end
+
+function Base.show(io::IO, ::MIME"text/plain", obs::ObservableInfo)
+    # Reverse the dictionary
+    observable_names = Dict(value => key for (key, value) in obs.observable_ixs)
+
+    for i = 1:length(obs.observables)
+        print(io,i == 1 ? "╔ " : i == length(obs.observables) ? "╚ " : "║ ")
+        for j = 1:length(obs.observables)
+            if i > j
+                print(io,"⋅ ")
+            elseif haskey(obs.correlations,CartesianIndex(i,j))
+                print(io,"⬤ ")
+            else
+                print(io,"• ")
+            end
+        end
+        print(io,observable_names[i])
+        println(io)
+    end
+    printstyled(io,"")
+end
+
+
+function parse_observables(N;observables = nothing,correlations = nothing)
+    # Set up correlation functions (which matrix elements αβ to save from 𝒮^{αβ})
+    if isnothing(observables)
+        # 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
+        # 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
+
+    # By default, include all correlations
+    if isnothing(correlations)
+        correlations = []
+        for oi in keys(observable_ixs), oj in keys(observable_ixs)
+            push!(correlations, (oi, oj))
+        end
+    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
+
+    # 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
+    ObservableInfo(observables, observable_ixs, idxinfo)
+end
+
+
+
+# Finds the linear index according to sc.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(obs::ObservableInfo,corrs; err_msg = αβ -> "Missing correlation $(αβ)")
+    indices = Vector{Int64}(undef,length(corrs))
+    for (i,(α,β)) in enumerate(corrs)
+        αi = obs.observable_ixs[α]
+        βi = obs.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(αβ)),obs.correlations,idx)
+    end
+    indices
+end
+
+function all_observable_names(obs::ObservableInfo)
+    observable_names = Dict(value => key for (key, value) in obs.observable_ixs)
+    [observable_names[i] for i in 1:length(observable_names)]
+end
+
+num_correlations(obs::ObservableInfo) = length(obs.correlations)
+num_observables(obs::ObservableInfo) = length(obs.observables)
+
+

src/Operators/Stevens.jl

@@ -206,6 +206,10 @@ print_stevens_expansion(S[1]^4 + S[2]^4 + S[3]^4)
 ```
 """
 function print_stevens_expansion(op::Matrix{ComplexF64})
+  println(show_stevens_expansion(op))
+end
+
+function show_stevens_expansion(op::Matrix{ComplexF64})
     op ≈ op' || error("Requires Hermitian operator")
     terms = String[]
 
@@ -226,5 +230,6 @@ function print_stevens_expansion(op::Matrix{ComplexF64})
     str = join(terms, " + ")
     # Remove redundant plus signs and print
     str = replace(str, "+ -" => "- ")
-    println(str)
 end
+
+

src/SampledCorrelations/CorrelationSampling.jl

@@ -49,12 +49,12 @@ function trajectory!(buf, sys, Δt, nsnaps, ops; measperiod = 1, apply_g = true)
 end
 
 function new_sample!(sc::SampledCorrelations, sys::System; processtraj! = no_processing)
-    (; Δt, samplebuf, measperiod, apply_g) = sc
+    (; Δt, samplebuf, measperiod, apply_g, observables) = sc
     nsnaps = size(samplebuf, 6)
 
     @assert size(sys.dipoles) == size(samplebuf)[2:5] "`System` size not compatible with given `SampledCorrelations`"
 
-    trajectory!(samplebuf, sys, Δt, nsnaps, sc.observables; measperiod = measperiod, apply_g = apply_g)
+    trajectory!(samplebuf, sys, Δt, nsnaps, observables.observables; measperiod = measperiod, apply_g = apply_g)
 
     processtraj!(sc)
 
@@ -81,7 +81,7 @@ function no_processing(::SampledCorrelations)
 end
 
 function accum_sample!(sc::SampledCorrelations)
-    (; data, variance, correlations, samplebuf, nsamples, fft!) = sc
+    (; data, variance, observables, samplebuf, nsamples, fft!) = sc
     natoms = size(samplebuf)[5]
 
     fft! * samplebuf # Apply pre-planned and pre-normalized FFT
@@ -91,7 +91,7 @@ function accum_sample!(sc::SampledCorrelations)
     # allocations here. The contents of the loop contains no allocations. There
     # does not seem to be a big performance penalty associated with these
     # allocations.
-    for j in 1:natoms, i in 1:natoms, (ci, c) in correlations  
+    for j in 1:natoms, i in 1:natoms, (ci, c) in observables.correlations  
         α, β = ci.I
 
         sample_α = @view samplebuf[α,:,:,:,i,:]

src/SampledCorrelations/DataRetrieval.jl

@@ -126,7 +126,7 @@ function intensity_formula(f::Function,sc::SampledCorrelations,required_correlat
     # 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(sc,required_correlations)
+    corr_ix = lookup_correlations(sc.observables,required_correlations)
     intensity_formula(f,sc,corr_ix;kwargs...)
 end