diff --git a/src/SampledCorrelations/DataRetrieval.jl b/src/SampledCorrelations/DataRetrieval.jl
index a5b5c01ff..5c301aeeb 100644
--- a/src/SampledCorrelations/DataRetrieval.jl
+++ b/src/SampledCorrelations/DataRetrieval.jl
@@ -106,7 +106,7 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
     NPos = Val{length(sc.crystal.positions)}()
 
     # Intensities calculation
-    intensities_rounded!(intensities, sc.data, sc.crystal, sc.positions, sc.measure.combiner, ffs, q_idx_info, ωidcs, NCorr, NPos)
+    intensities_aux!(intensities, sc.data, sc.crystal, sc.positions, sc.measure.combiner, ffs, q_idx_info, ωidcs, NCorr, NPos)
 
     # Convert to a q-space density in original (not reshaped) RLU.
     intensities .*= det(sc.crystal.recipvecs) / det(crystal.recipvecs)
@@ -136,20 +136,20 @@ function intensities(sc::SampledCorrelations, qpts; energies, kernel=nothing, kT
     end
 end
 
-function intensities_rounded!(intensities, data, crystal, positions, combiner, ff_atoms, q_idx_info, ωidcs, ::Val{NCorr}, ::Val{NPos}) where {NCorr, NPos}
+function intensities_aux!(intensities, data, crystal, positions, combiner, ff_atoms, q_idx_info, ωidcs, ::Val{NCorr}, ::Val{NPos}) where {NCorr, NPos}
     (; qabs, idcs, counts) = q_idx_info 
     (; recipvecs) = crystal 
     qidx = 1
     for (qabs, idx, count) in zip(qabs, idcs, counts)
-        prefactors = prefactors_for_phase_averaging(qabs, recipvecs, @view(positions[idx,:]), ff_atoms, Val(NCorr), Val(NPos))
+        prefactors = prefactors_for_phase_averaging(qabs, recipvecs, view(positions, idx, :), ff_atoms, Val{NCorr}(), Val{NPos}())
 
         # Perform phase-averaging over all omega
         for (n, iω) in enumerate(ωidcs)
-            elems = zero(MVector{NCorr,ComplexF64})
+            elems = zero(SVector{NCorr, ComplexF64})
             for j in 1:NPos, i in 1:NPos
-                elems .+= (prefactors[i] * conj(prefactors[j])) .* view(data, :, i, j, idx, iω)
+                elems += (prefactors[i] * conj(prefactors[j])) * SVector{NCorr}(view(data, :, i, j, idx, iω))
             end
-            val = combiner(qabs, SVector{NCorr, ComplexF64}(elems))
+            val = combiner(qabs, elems)
             intensities[n, qidx] = val
         end
 
diff --git a/src/SampledCorrelations/PhaseAveraging.jl b/src/SampledCorrelations/PhaseAveraging.jl
index 96f528ec1..fe1f352db 100644
--- a/src/SampledCorrelations/PhaseAveraging.jl
+++ b/src/SampledCorrelations/PhaseAveraging.jl
@@ -1,21 +1,3 @@
-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), Val{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])), Val{NAtoms}())
-
-    for j in 1:NAtoms, i in 1:NAtoms
-        elems .+= (prefactor[i] * conj(prefactor[j])) .* view(data, :, i, j)
-    end
-
-    return SVector{NCorr,ComplexF64}(elems)
-end
-
 function prefactors_for_phase_averaging(q_absolute::Vec3, recipvecs, positions, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms}
     # Form factor
     ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), Val{NAtoms}())
@@ -38,12 +20,12 @@ end
 # not used, but would be employed in a parallel intensities-like pipeline for
 # error propagation.
 function error_basis_reduction(data, q_absolute::Vec3, _::Crystal, ff_atoms, ::Val{NCorr}, ::Val{NAtoms}) where {NCorr, NAtoms}
-    elems = zero(MVector{NCorr, ComplexF64})
+    elems = zero(SVector{NCorr, ComplexF64})
     ffs = ntuple(i -> compute_form_factor(ff_atoms[i], q_absolute⋅q_absolute), NAtoms)
 
     for j in 1:NAtoms, i in 1:NAtoms
-        elems .+= ffs[i] .* ffs[j] .* view(data, :, i, j)
+        elems += (ffs[i] * ffs[j]) * SVector{NCorr}(view(data, :, i, j))
     end
 
-    return SVector{NCorr,Float64}(elems / (NAtoms*NAtoms))
+    return elems / (NAtoms*NAtoms)
 end