diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index 070c38d85..93d5383c1 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -729,7 +729,7 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
         auto restraints_now = [&]() {
             const RefinementStage stage("repulsions");
             RestraintSet set = topo.geometry;
-            set.Append(BuildRepulsions(topo, ctx.Model(), cell, sg, REPULSION_MARGIN));
+            set.Append(BuildRepulsions(topo, ctx.Model(), cell, sg, REPULSION_MARGIN, ctx.nthreads));
             set.Finalize(ctx.atoms.size());
             return set;
         };
@@ -737,7 +737,7 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
         {
             const RestraintSet set = restraints_now();
             ctx.engine->SetRestraints(0, set);
-            RestraintSet repulsions = BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0);
+            RestraintSet repulsions = BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0, ctx.nthreads);
             result.geometry_start = MeasureGeometry(topo.geometry, repulsions, ctx.XyzB());
         }
         double w0 = 0;
@@ -783,7 +783,7 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
             t.r_free = rf.r_free;
             {
                 const RefinementStage s1("measure");
-                t.geometry = MeasureGeometry(topo.geometry, BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0), ctx.XyzB());
+                t.geometry = MeasureGeometry(topo.geometry, BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0, ctx.nthreads), ctx.XyzB());
             }
             t.within_bounds = t.geometry.rmsd_bonds <= MAX_RMSD_BONDS && t.geometry.rmsd_angles <= MAX_RMSD_ANGLES;
             logger.Info("Refinement, X-ray weight {:.3g}, B {}: R-work {:.4f}, R-free {:.4f}, rmsd bonds {:.4f} A, "
diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp
index d1c8e698a..1f94562e3 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.cpp
+++ b/image_analysis/structure_refinement/RefinementRestraints.cpp
@@ -13,6 +13,7 @@
 #include "gemmi/riding_h.hpp"   // place_hydrogens_on_all_atoms
 #include "gemmi/topo.hpp"
 
+#include "../../common/ParallelFor.h"
 #include "MonomerLibrary.h"
 #include "RefinementMath.h"
 
@@ -457,78 +458,88 @@ RefinementTopology BuildRefinementTopology(gemmi::Structure &st, bool hydrogens)
     return out;
 }
 
-RestraintSet BuildRepulsions(const RefinementTopology &topo, const gemmi::Model &model, const gemmi::UnitCell &cell_in,
-                             const gemmi::SpaceGroup &sg, double margin) {
-    RestraintSet out;
+RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell_in,
+                             const gemmi::SpaceGroup &sg, double margin, size_t nthreads) {
     gemmi::UnitCell cell = cell_in;
     cell.set_cell_images_from_spacegroup(&sg);
-    gemmi::Model copy = model;
     float max_radius = 0;
     for (float r : topo.vdw_radius)
         max_radius = std::max(max_radius, r);
     const double reach = 2 * max_radius + margin;
-    gemmi::NeighborSearch ns(copy, cell, std::max(reach, 5.0));
+    gemmi::NeighborSearch ns(model, cell, std::max(reach, 5.0));
     ns.populate(false);
 
-    std::vector<std::vector<int>> first(copy.chains.size());
+    std::vector<std::vector<int>> first(model.chains.size());
     std::vector<gemmi::Position> pos;
-    for (size_t c = 0; c < copy.chains.size(); c++)
-        for (const gemmi::Residue &r : copy.chains[c].residues) {
+    for (size_t c = 0; c < model.chains.size(); c++)
+        for (const gemmi::Residue &r : model.chains[c].residues) {
             first[c].push_back(static_cast<int>(pos.size()));
             for (const gemmi::Atom &a : r.atoms)
                 pos.push_back(a.pos);
         }
 
-    for (size_t i = 0; i < pos.size(); i++) {
-        const gemmi::Fractional fi = cell.fractionalize(pos[i]);
-        ns.for_each(pos[i], '\0', reach, [&](gemmi::NeighborSearch::Mark &m, double) {
-            const int j = first[m.chain_idx][m.residue_idx] + m.atom_idx;
-            if (topo.rider_parent[i] >= 0 || topo.rider_parent[j] >= 0 || topo.element[i] == gemmi::El::H ||
-                topo.element[j] == gemmi::El::H)
-                return;
-            if (!topo.refinable[i] && !topo.refinable[j])
-                return;
-            if (m.image_idx == 0 && (static_cast<int>(i) >= j ||
-                                     std::binary_search(topo.excluded[i].begin(), topo.excluded[i].end(), j)))
-                return;
-            if (m.image_idx != 0 && static_cast<int>(i) > j)
-                return;
-            if (topo.altloc[i] && topo.altloc[j] && topo.altloc[i] != topo.altloc[j])
-                return;
-            // The image of j nearest to i, as a Cartesian x -> R x + t.
-            const gemmi::FTransform ft = ns.get_image_transformation(m.image_idx);
-            const gemmi::Fractional fj = ft.apply(cell.fractionalize(pos[j]));
-            const gemmi::Vec3 shift(std::round(fi.x - fj.x), std::round(fi.y - fj.y), std::round(fi.z - fj.z));
-            if (m.image_idx == 0 && (shift.x != 0 || shift.y != 0 || shift.z != 0))
-                return;   // the same copy's own lattice translate: the identity term covers the contact
-            const gemmi::Mat33 rot = cell.orth.mat.multiply(ft.mat).multiply(cell.frac.mat);
-            const gemmi::Vec3 tran = cell.orth.mat.multiply(ft.vec + shift);
-            const gemmi::Position image(rot.multiply(pos[j]) + tran);
-            const double d = pos[i].dist(image);
+    // The atoms in blocks searched on the threads, each block's pairs listed in atom order; the blocks are
+    // joined in order below, so the list is the same for any number of threads.
+    constexpr int BLOCKS = 64;
+    std::vector<std::vector<RestraintRecord>> found(BLOCKS);
+    ParallelFor(BLOCKS, nthreads, [&](int block) {
+        const size_t lo = pos.size() * block / BLOCKS, hi = pos.size() * (block + 1) / BLOCKS;
+        for (size_t i = lo; i < hi; i++) {
+            if (topo.rider_parent[i] >= 0 || topo.element[i] == gemmi::El::H)
+                continue;
+            const gemmi::Fractional fi = cell.fractionalize(pos[i]);
+            ns.for_each(pos[i], '\0', reach, [&](gemmi::NeighborSearch::Mark &m, double) {
+                const int j = first[m.chain_idx][m.residue_idx] + m.atom_idx;
+                if (topo.rider_parent[j] >= 0 || topo.element[j] == gemmi::El::H)
+                    return;
+                if (!topo.refinable[i] && !topo.refinable[j])
+                    return;
+                if (m.image_idx == 0 && (static_cast<int>(i) >= j ||
+                                         std::binary_search(topo.excluded[i].begin(), topo.excluded[i].end(), j)))
+                    return;
+                if (m.image_idx != 0 && static_cast<int>(i) > j)
+                    return;
+                if (topo.altloc[i] && topo.altloc[j] && topo.altloc[i] != topo.altloc[j])
+                    return;
+                // The image of j nearest to i, as a Cartesian x -> R x + t.
+                const gemmi::FTransform ft = ns.get_image_transformation(m.image_idx);
+                const gemmi::Fractional fj = ft.apply(cell.fractionalize(pos[j]));
+                const gemmi::Vec3 shift(std::round(fi.x - fj.x), std::round(fi.y - fj.y), std::round(fi.z - fj.z));
+                if (m.image_idx == 0 && (shift.x != 0 || shift.y != 0 || shift.z != 0))
+                    return;   // the same copy's own lattice translate: the identity term covers the contact
+                const gemmi::Mat33 rot = cell.orth.mat.multiply(ft.mat).multiply(cell.frac.mat);
+                const gemmi::Vec3 tran = cell.orth.mat.multiply(ft.vec + shift);
+                const gemmi::Position image(rot.multiply(pos[j]) + tran);
+                const double d = pos[i].dist(image);
 
-            float d0 = topo.vdw_radius[i] + topo.vdw_radius[j];
-            if (topo.metal[i] || topo.metal[j]) {
-                const gemmi::El other = topo.element[topo.metal[i] ? j : static_cast<int>(i)];
-                if (other == gemmi::El::N || other == gemmi::El::O || other == gemmi::El::S)
-                    d0 = METAL_CONTACT;
-            } else if (topo.hbond[i] && topo.hbond[j])
-                d0 = HBOND_CONTACT;
-            if (!(d < d0 + margin))
-                return;
-            RestraintRecord r;
-            r.kind = static_cast<int>(RestraintKind::Repulsion);
-            r.atom[0] = static_cast<int>(i);
-            r.atom[1] = j;
-            r.ideal = d0;
-            r.sigma = REPULSION_SIGMA;
-            for (int a = 0; a < 3; a++) {
-                for (int b = 0; b < 3; b++)
-                    r.rot[3 * a + b] = static_cast<float>(rot[a][b]);
-                r.tran[a] = static_cast<float>(tran.at(a));
-            }
+                float d0 = topo.vdw_radius[i] + topo.vdw_radius[j];
+                if (topo.metal[i] || topo.metal[j]) {
+                    const gemmi::El other = topo.element[topo.metal[i] ? j : static_cast<int>(i)];
+                    if (other == gemmi::El::N || other == gemmi::El::O || other == gemmi::El::S)
+                        d0 = METAL_CONTACT;
+                } else if (topo.hbond[i] && topo.hbond[j])
+                    d0 = HBOND_CONTACT;
+                if (!(d < d0 + margin))
+                    return;
+                RestraintRecord r;
+                r.kind = static_cast<int>(RestraintKind::Repulsion);
+                r.atom[0] = static_cast<int>(i);
+                r.atom[1] = j;
+                r.ideal = d0;
+                r.sigma = REPULSION_SIGMA;
+                for (int a = 0; a < 3; a++) {
+                    for (int b = 0; b < 3; b++)
+                        r.rot[3 * a + b] = static_cast<float>(rot[a][b]);
+                    r.tran[a] = static_cast<float>(tran.at(a));
+                }
+                found[block].push_back(r);
+            });
+        }
+    });
+    RestraintSet out;
+    for (const std::vector<RestraintRecord> &records : found)
+        for (const RestraintRecord &r : records)
             out.Add(r, 2);
-        });
-    }
     out.Finalize(pos.size());
     return out;
 }
diff --git a/image_analysis/structure_refinement/RefinementRestraints.h b/image_analysis/structure_refinement/RefinementRestraints.h
index 64b2e2ce2..7617e9094 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.h
+++ b/image_analysis/structure_refinement/RefinementRestraints.h
@@ -84,9 +84,10 @@ void PlaceRidingHydrogens(RefinementTopology &topo);
 // The repulsions of the model as it is now: every pair of atoms, symmetry mates included, closer than its
 // contact distance plus `margin` and not excluded, where at least one of the two is refinable and the two are
 // not alternative conformations of each other. Contact distances: van der Waals radii summed, 2.6 A where
-// both can hydrogen-bond, 1.9 A between a metal and N, O or S.
-RestraintSet BuildRepulsions(const RefinementTopology &topo, const gemmi::Model &model, const gemmi::UnitCell &cell,
-                             const gemmi::SpaceGroup &sg, double margin);
+// both can hydrogen-bond, 1.9 A between a metal and N, O or S. `model` is not changed (gemmi's NeighborSearch
+// takes it non-const); the search runs on `nthreads` threads.
+RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model, const gemmi::UnitCell &cell,
+                             const gemmi::SpaceGroup &sg, double margin, size_t nthreads = 1);
 
 // The usual summary of a model's geometry.
 struct GeometryStatistics {