Optimize phase averaging

SunnySuite · Oct 31, 2023 · 3c43a43 · 3c43a43
1 parent dee8b4f
commit 3c43a43
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 6 deletions.
diff --git a/src/Intensities/Interpolation.jl b/src/Intensities/Interpolation.jl
@@ -148,13 +148,13 @@ function intensities_interpolated(sc::SampledCorrelations, qs, formula::Classica
     intensities = zeros(return_type, size(qs)..., nω)
 
     # Call type-stable version of the function
-    intensities_interpolated!(intensities, sc, qs, ωvals, interp, formula, stencil_info, return_type)
+    intensities_interpolated!(intensities, sc, qs, ωvals, interp, formula, stencil_info, Val(return_type))
 
     return intensities
 end
 
 
-function intensities_interpolated!(intensities, sc::SampledCorrelations, q_targets::Array, ωvals, interp::InterpolationScheme{NInterp}, formula, stencil_info, T) where {NInterp}
+function intensities_interpolated!(intensities, sc::SampledCorrelations, q_targets::Array, ωvals, interp::InterpolationScheme{NInterp}, formula, stencil_info, ::Val{T}) where {NInterp, T}
     li_intensities = LinearIndices(intensities)
     ci_targets = CartesianIndices(q_targets)
 

diff --git a/src/SampledCorrelations/BasisReduction.jl b/src/SampledCorrelations/BasisReduction.jl
@@ -1,13 +1,16 @@
 function phase_averaged_elements(data, q_absolute::Vec3, crystal::Crystal, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms}
     elems = zero(MVector{NCorr,ComplexF64})
+
+    # Form factor
     ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), NAtoms)
 
-    # Real space position of each atom within the unit cell
-    rs = ntuple(i -> crystal.latvecs * crystal.positions[i], NAtoms)
+    # Overall phase factor for each site
+    q = crystal.recipvecs \ q_absolute
+    r = crystal.positions
+    prefactor = ntuple(i -> ffs[i] * exp(- 2π*im * (q ⋅ r[i])), NAtoms)
 
     for j in 1:NAtoms, i in 1:NAtoms
-        phase = exp(im*(q_absolute ⋅ (rs[j] - rs[i])))
-        elems .+= phase .* ffs[i] .* ffs[j] .* view(data, :, i, j)
+        elems .+= (prefactor[i] * conj(prefactor[j])) .* view(data, :, i, j)
     end
 
     return SVector{NCorr,ComplexF64}(elems)