diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index 184e671ed..7e8034a1e 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -1201,9 +1201,11 @@ std::vector<MissingReflection> MissingReflections(Context &ctx, const gemmi::SMa
     std::unique_ptr<RefinementEngine> engine = MakeRefinementEngineCPU(ctx.cell, *ctx.sg, ctx.d_min, missing, ctx.nthreads);
     engine->SetBlur(RefinementBlur(ctx.cell, *ctx.sg, ctx.d_min, ctx.b_floor));
     engine->SetAniso(&ctx.aniso);
+    // the atoms and the solvent side by side: neither reads the other
     std::vector<std::complex<float>> fcalc, fmask;
-    engine->Fcalc(ctx.Model(), fcalc);
+    std::future<void> atoms = std::async(std::launch::async, [&] { engine->Fcalc(ctx.Model(), fcalc); });
     ctx.mask->Fmask(ctx.Model(), missing.hkl, fmask);
+    atoms.get();
 
     std::vector<std::pair<double, double>> d_by_s;   // (1/d^2, D) of the observed reflections, by resolution
     for (size_t m = 0; m < ctx.refl.Size(); m++)
@@ -1298,7 +1300,11 @@ void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double
         fwt_data.v.push_back({ctx.refl.hkl[m], std::polar(static_cast<float>(fwt), static_cast<float>(phi))});
     }
     result.maps = maps;
-    // The difference map and the 2mFo-DFc map beside the anomalous one below: none reads another.
+    // The reflections the data do not have, the difference map and the 2mFo-DFc map beside the anomalous one below:
+    // none reads another.
+    std::future<std::vector<MissingReflection>> missing = std::async(std::launch::async, [&ctx, &b_aniso] {
+        return MissingReflections(ctx, b_aniso);
+    });
     const std::array<int, 3> size = gemmi::get_size_for_hkl(diff, {{0, 0, 0}}, 3.0);
     std::future<void> difference_map = std::async(std::launch::async, [&diff, &size, &result] {
         result.difference_map = std::make_shared<gemmi::Grid<float>>(MapFromFPhi(gemmi::get_f_phi_on_grid<float>(diff, size, true)));
@@ -1354,7 +1360,7 @@ void MakeMaps(Context &ctx, const std::vector<MergedReflection> &merged, double
     }
     // The filled reflections after the observed ones: FWT/PHWT and the model's FC/PHIC, every other column missing.
     const float nan = std::numeric_limits<float>::quiet_NaN();
-    for (const MissingReflection &r : MissingReflections(ctx, b_aniso)) {
+    for (const MissingReflection &r : missing.get()) {
         const double phi_deg = std::arg(r.f_model) * 180.0 / PI;
         const double k_aniso = std::exp(-0.25 * b_aniso.r_u_r(r.hkl));
         const float row[] = {static_cast<float>(r.hkl[0]), static_cast<float>(r.hkl[1]), static_cast<float>(r.hkl[2]), nan,