diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp
index 6470a8a0..8994d5d7 100644
--- a/image_analysis/scale_merge/RotationScaleMerge.cpp
+++ b/image_analysis/scale_merge/RotationScaleMerge.cpp
@@ -3182,19 +3182,29 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
     const bool pass_filters = (for_search && search_min_zeta > 0.0) || min_cc_for_image > 0.0;
     if (pass_filters) {
         corr_before_pass_filters.resize(partials.size());
-        for (size_t i = 0; i < partials.size(); ++i)
-            corr_before_pass_filters[i] = partials[i].corr;
+        ParallelChunks(static_cast<int>(partials.size()),
+                       ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+            for (int i = lo; i < hi; ++i) corr_before_pass_filters[i] = partials[i].corr;
+        });
     }
     bool corr_filtered = false;
 
     // --- 2a. On the de-novo search pass only, drop observations whose Lorentz geometry is poor. ---
     if (for_search && search_min_zeta > 0.0) {
-        int64_t n_dropped = 0;
-        for (auto &o : partials)
-            if (!(std::isfinite(o.zeta) && o.zeta >= search_min_zeta)) {
-                if (std::isfinite(o.corr) && o.corr > 0.0f) ++n_dropped;
-                o.corr = 0.0f;
+        std::atomic<int64_t> dropped = 0;
+        ParallelChunks(static_cast<int>(partials.size()),
+                       ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+            int64_t local = 0;
+            for (int i = lo; i < hi; ++i) {
+                Obs &o = partials[i];
+                if (!(std::isfinite(o.zeta) && o.zeta >= search_min_zeta)) {
+                    if (std::isfinite(o.corr) && o.corr > 0.0f) ++local;
+                    o.corr = 0.0f;
+                }
             }
+            dropped += local;
+        });
+        const int64_t n_dropped = dropped.load();
         corr_filtered = corr_filtered || n_dropped > 0;
         if (n_dropped > 0)
             logger.Info("Space-group search: ignoring {} observations with |zeta| < {:.2f} "
@@ -3212,8 +3222,11 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
                 ++n_rejected;
             }
         if (n_rejected > 0) {
-            for (auto &o : partials)
-                if (reject[o.frame]) o.corr = 0.0f;
+            ParallelChunks(static_cast<int>(partials.size()),
+                           ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+                for (int i = lo; i < hi; ++i)
+                    if (reject[partials[i].frame]) partials[i].corr = 0.0f;
+            });
             corr_filtered = true;
             logger.Info("Rejected {} of {} frames correlating below {:.2f} with the merged reference",
                         n_rejected, n_frames, min_cc_for_image);
@@ -3226,7 +3239,10 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
     // left it doing nothing whenever the CPU combine was in use (the diagnostic dump).
     if (gpu_active_ && corr_filtered) {
         std::vector<float> corr(partials.size());
-        for (size_t i = 0; i < partials.size(); ++i) corr[i] = partials[i].corr;
+        ParallelChunks(static_cast<int>(partials.size()),
+                       ThreadsForWork(partials.size(), nthreads), [&](int lo, int hi) {
+            for (int i = lo; i < hi; ++i) corr[i] = partials[i].corr;
+        });
         gpu_->SetCorr(corr.data());
     }
 #endif
@@ -3280,7 +3296,8 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
         gpu_->GetFullsPxPy(fpx.data(), fpy.data());
         gpu_->GetFullsVariance(fvb.data(), fvi.data());
         if (scaled_fulls_on_gpu) gpu_->GetFullsCorr(fcorr.data());
-        for (int i = 0; i < nf; ++i) {
+        ParallelChunks(nf, ThreadsForWork(static_cast<size_t>(nf), nthreads), [&](int lo, int hi) {
+        for (int i = lo; i < hi; ++i) {
             Obs &o = fulls[i];
             o.h = fh[i]; o.k = fk[i]; o.l = fl[i];
             o.I = fI[i]; o.sigma = fsig[i]; o.d = fd[i];
@@ -3289,6 +3306,7 @@ RotationScaleMerge::Result RotationScaleMerge::Run(bool for_search) {
             o.image_number = fimg[i]; o.frame = fframe[i]; o.px = fpx[i]; o.py = fpy[i];
             o.on_ice = fon[i]; o.group = fgroup[i];
         }
+        });
         logger.Info("3D combine{} (GPU): {} fulls", scaled_fulls_on_gpu ? " + scale-fulls" : "", nf);
         combined_on_gpu = true;
     }