diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index ea8771af8..560d3ff93 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -208,7 +208,7 @@ struct Context {
         }
         {
             const RefinementStage s4("refit.ml");
-            EstimateLikelihoodParameters(refl, fmodel, rows);
+            EstimateLikelihoodParameters(refl, fmodel, rows, nthreads);
         }
         engine->SetRowModel(rows);
         return ComputeRFactors(refl, fmodel);
diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp
index ccbe8e3ab..61120304d 100644
--- a/image_analysis/structure_refinement/RefinementXray.cpp
+++ b/image_analysis/structure_refinement/RefinementXray.cpp
@@ -229,7 +229,7 @@ RefinementReflections PrepareRefinementReflections(const std::vector<MergedRefle
 }
 
 void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::vector<std::complex<float>> &fmodel,
-                                  XrayRowModel &rows) {
+                                  XrayRowModel &rows, size_t nthreads) {
     const size_t n = refl.Size();
     std::vector<double> fc(n);
     for (size_t i = 0; i < n; i++)
@@ -248,7 +248,8 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
     const size_t bins = std::clamp<size_t>(free.size() / FREE_PER_BIN, 1, MAX_BINS);
 
     std::vector<double> centre(bins), bin_d(bins), bin_beta(bins);
-    for (size_t b = 0; b < bins; b++) {
+    // each bin is fitted on its own
+    ParallelFor(static_cast<int>(bins), nthreads, [&](int b) {
         const std::vector<size_t> idx(free.begin() + b * free.size() / bins, free.begin() + (b + 1) * free.size() / bins);
         double s2 = 0, sfo = 0, sfc2 = 0, sfofc = 0;
         for (size_t i : idx) {
@@ -271,7 +272,7 @@ void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::
         }
         bin_d[b] = d;
         bin_beta[b] = beta;
-    }
+    });
     rows.d.resize(n);
     rows.var.resize(n);
     for (size_t i = 0; i < n; i++) {
diff --git a/image_analysis/structure_refinement/RefinementXray.h b/image_analysis/structure_refinement/RefinementXray.h
index a0a459f21..fc37a9e85 100644
--- a/image_analysis/structure_refinement/RefinementXray.h
+++ b/image_analysis/structure_refinement/RefinementXray.h
@@ -58,7 +58,7 @@ struct XrayRowModel {
 // epsilon beta + sigma^2. `fmodel` is k_total (F + fmask), without D.
 // Following Lunin & Skovoroda (1995) Acta Cryst. A51, 880-887, and Afonine et al. (2012) Acta Cryst. D68, 352-367
 void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::vector<std::complex<float>> &fmodel,
-                                  XrayRowModel &rows);
+                                  XrayRowModel &rows, size_t nthreads = 1);
 
 // R-work and R-free of k_total (F + fmask) against the observations.
 struct RFactors {