Support Asymmetric Correlations

src/Intensities/ElementContraction.jl

@@ -19,6 +19,9 @@ struct FullTensor{NCorr,NSquare,NObs,NObs2} <: Contraction{SMatrix{NObs, NObs, C
     indices :: SVector{NSquare, Int64}
 end
 
+struct AllAvailable{NCorr} <: Contraction{SVector{NCorr, ComplexF64}}
+end
+
 
 ################################################################################
 # Constructors
@@ -52,7 +55,7 @@ function Trace(obs::ObservableInfo)
     Trace(SVector{length(indices), Int64}(indices))
 end
 
-
+# Warning: This dipole factor assumes Sαβ symmetric (i.e. it assumes equilibrium)
 function DipoleFactor(obs::ObservableInfo; spin_components = [:Sx,:Sy,:Sz])
     # Ensure that the observables themselves are present
     for si in spin_components
@@ -119,9 +122,11 @@ end
 contract(specific_element, _, ::Element) = only(specific_element)
 
 function contract(all_elems, _, full::FullTensor{NCorr,NSquare,NObs,NObs2}) where {NCorr, NSquare,NObs,NObs2}
-    # 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))
+    reshape(all_elems[full.indices],NObs,NObs)
+end
+
+function contract(all_elems, _, ::AllAvailable{NCorr}) where NCorr
+    all_elems
 end
 
 ################################################################################
@@ -131,6 +136,7 @@ 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
+required_correlations(::AllAvailable{NCorr}) where NCorr = 1:NCorr
 
 
 ################################################################################
@@ -144,10 +150,14 @@ function contractor_from_mode(source, mode::Symbol)
         string_formula = "Tr S"
     elseif mode == :perp
         contractor = DipoleFactor(source.observables)
-        string_formula = "∑_ij (I - Q⊗Q){i,j} S{i,j}\n\n(i,j = Sx,Sy,Sz)"
+        string_formula = "∑_ij (I - Q⊗Q){i,j} S{i,j}\n\n(i,j = Sx,Sy,Sz) and assuming S = S†"
     elseif mode == :full
         contractor = FullTensor(source.observables)
         string_formula = "S{α,β}"
+    elseif mode == :all_available
+        corrs = keys(source.observables.correlations)
+        contractor = AllAvailable{length(corrs)}()
+        string_formula = "[" * join(map(x -> "S{$(x.I[1]),$(x.I[2])}",corrs),", ") * "]"
     end
     return contractor, string_formula
 end
@@ -160,6 +170,7 @@ Sunny has several built-in formulas that can be selected by setting `contraction
 - `:trace` (default), which yields ``\\operatorname{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.
+- `:all_available`, which will return each computed correlation as determined by `sc.observables` or `swt.observables`.
 """
 function intensity_formula(swt::SpinWaveTheory, mode::Symbol; kwargs...)
     contractor, string_formula = contractor_from_mode(swt, mode)

src/Operators/Observables.jl

@@ -13,9 +13,7 @@ function Base.show(io::IO, ::MIME"text/plain", obs::ObservableInfo)
     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))
+            if haskey(obs.correlations,CartesianIndex(i,j))
                 print(io,"⬤ ")
             else
                 print(io,"• ")
@@ -101,9 +99,6 @@ function parse_observables(N; observables, correlations)
     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
@@ -121,12 +116,10 @@ function lookup_correlations(obs::ObservableInfo,corrs; err_msg = αβ -> "Missi
     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)
+        indices[i] = get!(() -> error(err_msg((α,β))),obs.correlations,idx)
     end
     indices
 end

src/SampledCorrelations/CorrelationSampling.jl

@@ -84,7 +84,11 @@ function accum_sample!(sc::SampledCorrelations)
     (; data, variance, observables, samplebuf, nsamples, fft!) = sc
     natoms = size(samplebuf)[5]
 
+    sourcebuf = copy(samplebuf)
+    sourcebuf[:,:,:,:,:,(1+size(samplebuf,6)÷2):end] .= 0
+
     fft! * samplebuf # Apply pre-planned and pre-normalized FFT
+    fft! * sourcebuf
     count = nsamples[1] += 1
 
     # Note that iterating over the `correlations` (a SortedDict) causes
@@ -94,14 +98,19 @@ function accum_sample!(sc::SampledCorrelations)
     for j in 1:natoms, i in 1:natoms, (ci, c) in observables.correlations  
         α, β = ci.I
 
+        # Convention is: S{α,β} means <α(t)> <β(0)>, i.e. α is delayed
         sample_α = @view samplebuf[α,:,:,:,i,:]
-        sample_β = @view samplebuf[β,:,:,:,j,:]
+        sample_β = @view sourcebuf[β,:,:,:,j,:]
+
+        correlation_plus_trash = FFTW.ifft(sample_α .* conj.(sample_β),4)
+        correlation = FFTW.fft(correlation_plus_trash[:,:,:,1:size(data,7)],4)
+
         databuf  = @view data[c,i,j,:,:,:,:]
 
         if isnothing(variance)
             for k in eachindex(databuf)
                 # Store the diff for one complex number on the stack.
-                diff = sample_α[k] * conj(sample_β[k]) - databuf[k]
+                diff = correlation[k] - databuf[k]
 
                 # Accumulate into running average
                 databuf[k] += diff * (1/count)
@@ -113,7 +122,7 @@ function accum_sample!(sc::SampledCorrelations)
                 μ_old = databuf[k]
 
                 # Update running mean.
-                matrixelem = sample_α[k] * conj(sample_β[k])
+                matrixelem = correlation[k]
                 databuf[k] += (matrixelem - databuf[k]) * (1/count)
                 μ = databuf[k]
 

src/SampledCorrelations/DataRetrieval.jl

@@ -115,7 +115,7 @@ The function is intended to be specified using `do` notation. For example, this
 
 If your custom formula returns a type other than `Float64`, use the `return_type` keyword argument to flag this.
 """
-function intensity_formula(f::Function,sc::SampledCorrelations,required_correlations; kwargs...)
+function intensity_formula(f::Function,sc,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

src/SampledCorrelations/SampledCorrelations.jl

@@ -160,12 +160,12 @@ function dynamical_correlations(sys::System{N}; Δt, nω, ωmax,
 
     # Preallocation
     na = natoms(sys.crystal)
-    samplebuf = zeros(ComplexF64, num_observables(observables), sys.latsize..., na, nω) 
+    samplebuf = zeros(ComplexF64, num_observables(observables), sys.latsize..., na, 2nω)
     data = zeros(ComplexF64, num_correlations(observables), na, na, sys.latsize..., nω)
     variance = calculate_errors ? zeros(Float64, size(data)...) : nothing
 
     # Normalize FFT according to physical convention
-    normalizationFactor = 1/(nω * √(prod(sys.latsize)))
+    normalizationFactor = 1/(nω * √(prod(sys.latsize))) # Note nω not 2nω due to zeroing
     fft! = normalizationFactor * FFTW.plan_fft!(samplebuf, (2,3,4,6))
 
     # Other initialization

test/test_correlation_sampling.jl

@@ -1,3 +1,40 @@
+@testitem "Asymmetric Correlations" begin
+
+    latsize = (1,1,1)
+    cryst = Sunny.fcc_primitive_crystal()
+    sys = System(cryst, latsize, [SpinInfo(1; S = 1/2, g=2)], :SUN; seed = 0)
+    sc = dynamical_correlations(sys; Δt = 0.1, ωmax = 10.0, nω=100, observables = [:A => I(2), :B => I(2)])
+    ts = range(0,1,length = size(sc.samplebuf,6))
+    #As = cos.(10π .* ts)
+    #Bs = 2 .* sin.(10π .* ts)
+    As = exp.(-(ts .- 0.15).^2 ./ (2 * 0.05^2))
+    Bs = exp.(-(ts .- 0.35).^2 ./ (2 * 0.1^2))
+    sc.samplebuf[1,1,1,1,1,:] .= As
+    sc.samplebuf[2,1,1,1,1,:] .= Bs
+    Sunny.accum_sample!(sc)
+    real_data = real(ifft(sc.data,7))
+    real_data .*= 199*199
+
+    # Reference calculation
+    q11 = zeros(199)
+    q12 = zeros(199)
+    q21 = zeros(199)
+    q22 = zeros(199)
+    for t = 0:198
+        for tau = 0:198
+            q11[1+t] += As[1+(t+tau)] * As[1+(tau)]
+            q12[1+t] += As[1+(t+tau)] * Bs[1+(tau)]
+            q21[1+t] += Bs[1+(t+tau)] * As[1+(tau)]
+            q22[1+t] += Bs[1+(t+tau)] * Bs[1+(tau)]
+        end
+    end
+
+    @test isapprox(q11,real_data[sc.observables.correlations[CartesianIndex(1,1)],1,1,1,1,1,:];atol = 1e-8)
+    @test isapprox(q12,real_data[sc.observables.correlations[CartesianIndex(1,2)],1,1,1,1,1,:];atol = 1e-8)
+    @test isapprox(q21,real_data[sc.observables.correlations[CartesianIndex(2,1)],1,1,1,1,1,:];atol = 1e-8)
+    @test isapprox(q22,real_data[sc.observables.correlations[CartesianIndex(2,2)],1,1,1,1,1,:];atol = 1e-8)
+end
+
 @testitem "Correlation sampling" begin
 
     function simple_model_sc(; mode, seed=111)