Merge branch 'af-mr' into rc175: experimental --af-model, single-copy molecular replacement deciding open symmetry hypotheses into p_mr.mtz

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
2026-10-11 14:38:08 +02:00
parent 62b7d8b2b9 60a8bf3981
commit 57fef1cf44
12 changed files with 1575 additions and 1 deletions
@@ -438,6 +438,35 @@ Problems Whose Variables Separate" (1973), SIAM J. Numer. Anal. 10, 413-432
 for solving separable nonlinear least squares problems" (1975), BIT 15, 49-57
 [doi:10.1007/BF01932995](https://doi.org/10.1007/BF01932995).

+**Predicted models** — `--af-model` (experimental) places an AlphaFold prediction by molecular
+replacement. J. Jumper, R. Evans, A. Pritzel et al., "Highly accurate protein structure prediction
+with AlphaFold" (2021), Nature 596, 583-589
+[doi:10.1038/s41586-021-03819-2](https://doi.org/10.1038/s41586-021-03819-2); the models are those of
+the AlphaFold Protein Structure Database: M. Varadi, S. Anyango, M. Deshpande et al., "AlphaFold
+Protein Structure Database: massively expanding the structural coverage of protein-sequence space
+with high-accuracy models" (2022), Nucleic Acids Res. 50, D439-D444
+[doi:10.1093/nar/gkab1061](https://doi.org/10.1093/nar/gkab1061). The model is prepared as
+phenix.process_predicted_model prepares one - residues below pLDDT 70 removed, B-factors from the
+coordinate error the pLDDT predicts: R. D. Oeffner, T. Croll, C. Millán et al., "Putting AlphaFold
+models to work with phenix.process_predicted_model and ISOLDE" (2022), Acta Cryst. D78, 1303-1314
+[doi:10.1107/S2059798322010026](https://doi.org/10.1107/S2059798322010026).
+
+**Molecular replacement** — the search scores each orientation of the model against the data by
+interpolating its molecular transform, computed once by FFT, as Navaza's AMoRe does: J. Navaza,
+"AMoRe: an automated package for molecular replacement" (1994), Acta Cryst. A50, 157-163
+[doi:10.1107/S0108767393007597](https://doi.org/10.1107/S0108767393007597). The translation
+function is the full-symmetry one, all pairs of symmetry copies on one Fourier synthesis:
+R. A. Crowther and D. M. Blow, "A method of positioning a known molecule in an unknown crystal
+structure" (1967), Acta Cryst. 23, 544-548
+[doi:10.1107/S0365110X67003172](https://doi.org/10.1107/S0365110X67003172). A placement is called
+a solution, as Phaser calls one, where its translation-function Z-score reaches 8: A. J. McCoy,
+R. W. Grosse-Kunstleve, P. D. Adams, M. D. Winn, L. C. Storoni and R. J. Read, "Phaser
+crystallographic software" (2007), J. Appl. Cryst. 40, 658-674
+[doi:10.1107/S0021889807021206](https://doi.org/10.1107/S0021889807021206). Whether one copy fits
+the cell is read off the Matthews coefficient: B. W. Matthews, "Solvent content of protein
+crystals" (1968), J. Mol. Biol. 33, 491-497
+[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2).
+
 ## Software and computing methods

 Decoding bitshuffle+LZ4 images on the GPU, rather than decompressing them on the host and uploading
@@ -109,7 +109,8 @@ The header block above section 1 records **how the result was produced**: `RUGNU
 Rates, per-image costs and progress remain on stdout only.

 Sections, in order: the `SUMMARY`, then `1. DATA SET AND GEOMETRY`, `2. CRYSTAL`,
-`3. MERGED DATA`, `4. DIAGNOSTICS`, and `5. MODEL VALIDATION` only with `--model`
+`3. MERGED DATA`, `4. DIAGNOSTICS`, `5. MODEL VALIDATION` only with `--model`, and the unnumbered
+`MODEL PRIOR (EXPERIMENTAL)` only with `--af-model` ([below](#model-prior---af-model-experimental))
 (see [The summary](#the-summary-and-what-the-run-decided) below). The numbering is contiguous, and
 the fifth section appearing renumbers nothing; a stage that did not run states that inside its
 section — `MERGE= NOT_PERFORMED` — rather than the section disappearing.
@@ -958,3 +959,34 @@ the scaling has mistreated.
 files are **byte for byte** what a run with no model would have written — same space group, same indexing, same `.mtz`, `.cif`, `.hkl` and
 `_unmerged.mtz`. A rejected model is therefore safe to try: it costs the null's compute and changes
 nothing else.
+
+## Model prior (`--af-model`, experimental)
+
+`--af-model <file>` takes a predicted model - an AlphaFold mmCIF, with the pLDDT in the B-factor
+column - and places it by molecular replacement: residues below pLDDT 70 removed, B-factors from the
+pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A), a translation function
+for the 50 best orientations in every space group the merge cannot tell apart, and the best three
+placements per group through the rigid body, bulk solvent and R-free (to `SCORE_D_MIN`, 3 A or the
+data's limit). One copy in the asymmetric unit only, for now.
+
+The table has one row per group and merge. A row is `PLACED` where its translation-function Z-score
+(`TF_Z`) reaches 8, at most 5% of its C-alpha atoms clash with symmetry mates (`CLASH`), one copy fits
+the cell (`MATTHEWS`: Vm below 1.5 A^3/Da; `MULTI_COPY`: Vm of 4 or more, which is reported but not
+supported yet) and its R-work stands at least 3 sd above the same model refitted from random
+orientations (`NULL_SIGMA`, or `NULL` where it does not). `MODEL_PRIOR_DECISION=` is then:
+
+- `AGREES` - the placed group with the lowest R-free is the answer's, ahead of every other group of
+  the merge that packs by 0.03 in R-free;
+- `ENANTIOMORPH` - it is the answer's enantiomorph, and the answer's files take that label (a label
+  only, as with `--model`; no reflection moves);
+- `SPACE_GROUP` - another group with the answer's rotations; written as `<prefix>_mr.mtz`/`.cif` from
+  the answer's merge, the reflections that group extinguishes left out;
+- `TWIN_SUBGROUP` - no group of the answer's merge can hold the model at all (`MATTHEWS`/`CLASH` in
+  every row) and a group of the twin hypothesis' subgroup can; written as `<prefix>_mr.mtz`/`.cif`
+  from the twin hypothesis' merge, with its law and fraction. The twin hypothesis' merge is searched
+  only in this case;
+- `UNDECIDED` or `NOT_PLACED` - nothing is written.
+
+With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `<prefix>_mr_model.pdb`/`.cif`.
+Apart from the enantiomorph label, the answer's own files are never changed by the model.
+
@@ -2,6 +2,7 @@ ADD_LIBRARY(JFJochStructureRefinement STATIC
        ModelFFT.cpp ModelFFT.h
        ModelGrid.cpp ModelGrid.h
        ModelScaling.cpp ModelScaling.h
+        MolecularReplacement.cpp MolecularReplacement.h
        ModelValidation.cpp ModelValidation.h
        RigidBodyRefine.cpp RigidBodyRefine.h
        $<$<BOOL:${JFJOCH_CUDA_AVAILABLE}>:RigidBodyGPU.cpp>
@@ -0,0 +1,994 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "MolecularReplacement.h"
+
+#include <algorithm>
+#include <array>
+#include <atomic>
+#include <chrono>
+#include <cmath>
+#include <future>
+#include <numeric>
+#include <random>
+#include <unordered_map>
+
+#include <gemmi/mmread_gz.hpp>   // read_structure_gz
+#include <gemmi/modify.hpp>      // remove_hydrogens
+#include <gemmi/dencalc.hpp>     // DensityCalculator
+#include <gemmi/it92.hpp>        // IT92 x-ray form factors
+#include <gemmi/scaling.hpp>     // Scaling (bulk solvent + anisotropic B)
+
+#include "../../common/JFJochMath.h"   // PI
+#include "../../common/Logger.h"
+#include "../../common/ParallelFor.h"
+#include "ModelFFT.h"
+#include "ModelGrid.h"
+#include "ModelScaling.h"
+#include "RigidBodyRefine.h"
+#ifdef JFJOCH_USE_CUDA
+#include "RigidBodyGPU.h"
+#endif
+
+namespace {
+
+using Table = gemmi::IT92<float>;
+
+// The search's resolution range: the low-resolution end is where the bulk solvent, which the search
+// does not model, stops dominating; the high-resolution end is where a predicted model still agrees
+// with the crystal's structure in its details.
+constexpr double SEARCH_D_MIN = 4.0;
+constexpr double SEARCH_D_MAX = 12.0;
+constexpr int SEARCH_BINS = 10;
+constexpr double ROTATION_STEP_DEG = 6.0;
+// Orientations the translation function is run for, per merge, and placements per space group taken
+// through the rigid body.
+constexpr int TRANSLATION_ROTATIONS = 50;
+constexpr int RESCORED_PER_GROUP = 3;
+// The resolution the placements are refined and scored to: whether a model sits in the right place
+// is a low-resolution question, the rigid body's ladder ends at 3.5 A anyway, and every R-factor of
+// the comparison is made to the same limit.
+constexpr double SCORE_D_MIN = 3.0;
+
+// The decision. A translation-function peak this many standard deviations above the function's mean
+// is a solution, where the model packs and one copy fits the cell; the McCoy et al. (2007) J. Appl.
+// Cryst. 40, 658-674 convention for a single-copy search, taken over to this Z-score's own null on
+// the Crowther & Blow function (measured: the highest peak of a wrong orientation reads 4-6).
+constexpr double TF_Z_PLACED = 8.0;
+constexpr double CLASH_MAX = 0.05;
+// Matthews (1968) J. Mol. Biol. 33, 491-497: protein crystals have Vm of about 1.7 to 3.5 A^3/Da.
+// Below VM_MIN one copy does not fit the cell; at or above VM_TWO_COPIES two copies would.
+constexpr double VM_MIN = 1.5;
+constexpr double VM_TWO_COPIES = 4.0;
+constexpr double DALTON_PER_RESIDUE = 110.0;
+// How far below every other candidate of its merge the chosen one's R-free must be.
+constexpr double R_FREE_MARGIN = 0.03;
+// The null (as in model validation): the same model refitted from random orientations about its
+// centroid, at a coarse limit, and how far above that the real placement's R-work has to sit.
+constexpr int NULL_REPLICATES = 9;
+constexpr unsigned NULL_SEED = 20261011;
+constexpr double NULL_D_MIN = 3.5;
+constexpr double NULL_SIGMA_MIN = 3.0;
+
+double seconds_since(std::chrono::steady_clock::time_point t0) {
+    return std::chrono::duration<double>(std::chrono::steady_clock::now() - t0).count();
+}
+
+// A rotation drawn uniformly from SO(3), through a uniform random unit quaternion.
+// Following Shoemake (1992) Graphics Gems III, 124-132
+gemmi::Mat33 random_rotation(std::mt19937 &rng) {
+    std::uniform_real_distribution<double> u(0.0, 1.0);
+    const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng);
+    const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1);
+    const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3);
+    return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w),     2 * (x * z + y * w),
+            2 * (x * y + z * w),     1 - 2 * (x * x + z * z), 2 * (y * z - x * w),
+            2 * (x * z - y * w),     2 * (y * z + x * w),     1 - 2 * (x * x + y * y)};
+}
+
+double trace(const gemmi::Mat33 &m) {
+    return m.a[0][0] + m.a[1][1] + m.a[2][2];
+}
+
+gemmi::Vec3 miller_vec(const std::array<int, 3> &h) {
+    return {static_cast<double>(h[0]), static_cast<double>(h[1]), static_cast<double>(h[2])};
+}
+
+// h R for a symmetry operator's rotation: the index the operator takes h to.
+gemmi::Miller rotated_index(const gemmi::Op &op, const gemmi::Miller &h) {
+    const auto k = op.apply_to_hkl_without_division(h);
+    return {{k[0] / gemmi::Op::DEN, k[1] / gemmi::Op::DEN, k[2] / gemmi::Op::DEN}};
+}
+
+long hkl_key(const gemmi::Miller &h) {
+    return (static_cast<long>(h[0]) + 512) * 1048576L + (static_cast<long>(h[1]) + 512) * 1024L + (h[2] + 512);
+}
+
+int good_even_size(int n) {
+    while (n % 2 != 0 || !gemmi::has_small_factorization(n))
+        ++n;
+    return n;
+}
+
+// --- placements scored against the data: R-work / R-free after the bulk solvent and the scale ---
+
+struct PlacementScore {
+    double r_work = 1, r_free = 1;
+};
+
+gemmi::AsuData<gemmi::ValueSigma<float>> working_amplitudes(const std::vector<MergedReflection> &merged,
+                                                            const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                                            double d_min) {
+    const gemmi::GroupOps gops = sg.operations();
+    const gemmi::ReciprocalAsu asu(&sg);
+    gemmi::AsuData<gemmi::ValueSigma<float>> fobs;
+    fobs.unit_cell_ = cell;
+    fobs.spacegroup_ = &sg;
+    for (const MergedReflection &r : merged) {
+        if (std::isnan(r.F) || r.rfree_flag || r.d < d_min)
+            continue;
+        gemmi::Miller h{{r.h, r.k, r.l}};
+        if (!asu.is_in(h))
+            h = asu.to_asu(h, gops).first;
+        fobs.v.push_back({h, {r.F, 1.0f}});
+    }
+    fobs.ensure_asu();
+    fobs.ensure_sorted();
+    return fobs;
+}
+
+// The model's Fcalc and bulk-solvent mask to d_min, scaled to the working amplitudes as model
+// validation scales them (ModelValidation.cpp, fit_model), and the R-factors over the reflections
+// to d_min.
+PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
+                               const std::vector<MergedReflection> &merged, double d_min, size_t nthreads) {
+    gemmi::DensityCalculator<Table, float> dc;
+    dc.d_min = d_min;
+    dc.rate = 1.5;
+    dc.set_grid_cell_and_spacegroup(st);
+    dc.set_refmac_compatible_blur(st.models[0]);
+    gemmi::Grid<float> mask;
+    mask.unit_cell = dc.grid.unit_cell;
+    mask.spacegroup = dc.grid.spacegroup;
+    mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
+    const std::vector<size_t> leaders = OrbitLeaders(mask, nthreads);
+    PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
+    gemmi::AsuData<std::complex<float>> fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
+    PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
+    gemmi::AsuData<std::complex<float>> fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
+
+    gemmi::AsuData<gemmi::ValueSigma<float>> fobs = working_amplitudes(merged, st.cell, sg, d_min);
+    gemmi::Scaling<float> scaling(st.cell, &sg);
+    scaling.use_solvent = true;
+    scaling.prepare_points(fmodel, fobs, &fmask);
+    FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
+    scaling.scale_data(fmodel, &fmask);
+
+    const gemmi::GroupOps gops = sg.operations();
+    const gemmi::ReciprocalAsu asu(&sg);
+    std::unordered_map<long, std::pair<double, bool>> obs;
+    for (const MergedReflection &r : merged) {
+        if (std::isnan(r.F) || r.d < d_min)
+            continue;
+        gemmi::Miller h{{r.h, r.k, r.l}};
+        if (!asu.is_in(h))
+            h = asu.to_asu(h, gops).first;
+        obs[hkl_key(h)] = {r.F, r.rfree_flag};
+    }
+    double num_w = 0, den_w = 0, num_f = 0, den_f = 0;
+    for (const auto &hv : fmodel.v) {
+        const auto it = obs.find(hkl_key(hv.hkl));
+        if (it == obs.end())
+            continue;
+        const double fo = it->second.first, fc = std::abs(hv.value);
+        if (it->second.second) {
+            num_f += std::fabs(fo - fc);
+            den_f += fo;
+        } else {
+            num_w += std::fabs(fo - fc);
+            den_w += fo;
+        }
+    }
+    PlacementScore out;
+    out.r_work = den_w > 0 ? num_w / den_w : 1;
+    out.r_free = den_f > 0 ? num_f / den_f : 1;
+    return out;
+}
+
+// The rigid-body engines for one model in one group, where there is a GPU.
+struct RigidBodyEngines {
+#ifdef JFJOCH_USE_CUDA
+    std::unique_ptr<RigidBodyGPUPool> pool;
+#endif
+    RigidBodyGPUPool *get() {
+#ifdef JFJOCH_USE_CUDA
+        return pool.get();
+#else
+        return nullptr;
+#endif
+    }
+};
+
+RigidBodyEngines make_engines(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
+                              const std::vector<MergedReflection> &merged, double d_min, size_t engines,
+                              Logger &logger) {
+    RigidBodyEngines out;
+#ifdef JFJOCH_USE_CUDA
+    const double finest = RigidBodyLadder(d_min).back();
+    size_t observations = 0;
+    for (const MergedReflection &r : merged)
+        if (!std::isnan(r.F) && r.d >= finest)
+            ++observations;
+    out.pool = RigidBodyGPUPool::Create(st.models[0], st.cell, sg, d_min, observations, engines, logger);
+#else
+    (void) st; (void) sg; (void) merged; (void) d_min; (void) engines; (void) logger;
+#endif
+    return out;
+}
+
+// The placement through the rigid body, and its R-factors where the rigid body left it. A search
+// placement is a few degrees and a grid step from its answer, well inside the rigid body's reach, so
+// the step is taken without the R-free check model validation makes before committing one.
+PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg,
+                                const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
+                                size_t nthreads) {
+    const std::vector<gemmi::Position> start = ModelPositions(st.models[0]);
+    Logger quiet = Logger::Buffered();
+    const RigidBodyRefineResult rb = RefineRigidBody(st.models[0], st.cell, sg,
+                                                     working_amplitudes(merged, st.cell, sg, d_min), d_min, quiet,
+                                                     nthreads, pool);
+    if (!rb.converged)
+        SetModelPositions(st.models[0], start);
+    return score_placement(st, sg, merged, d_min, nthreads);
+}
+
+// The placed model's R-work at NULL_D_MIN against the same model refitted from random orientations
+// about its centroid, in standard deviations of the latter: the null model validation builds for a
+// model that claims something (ModelValidation.cpp). An orientation the crystal cannot tell from the
+// placement's - within the rigid body's reach of one of its symmetry classes - is drawn again.
+double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg,
+                  const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
+                  size_t nthreads) {
+    const double null_d_min = std::max(d_min, NULL_D_MIN);
+    std::vector<MergedReflection> coarse;
+    for (const MergedReflection &r : merged)
+        if (r.d >= null_d_min)
+            coarse.push_back(r);
+    const double real = score_placement(placed, sg, coarse, null_d_min, nthreads).r_work;
+
+    const std::vector<gemmi::Mat33> pg = CartesianPointGroup(sg, placed.cell);
+    const double reach = RigidBodyReachDeg(placed.models[0], null_d_min);
+    std::vector<gemmi::Mat33> draws;
+    std::mt19937 rng(NULL_SEED);
+    while (static_cast<int>(draws.size()) < NULL_REPLICATES) {
+        const gemmi::Mat33 r = random_rotation(rng);
+        if (RotationDistanceDeg(r.multiply(orientation), orientation, pg) >= reach)
+            draws.push_back(r);
+    }
+    std::vector<double> r_work(NULL_REPLICATES);
+    std::atomic<int> next{0};
+    std::vector<std::future<void>> workers;
+    for (size_t t = 0; t < std::min<size_t>(std::max<size_t>(nthreads, 1), NULL_REPLICATES); t++)
+        workers.push_back(std::async(std::launch::async, [&] {
+            for (int i = next++; i < NULL_REPLICATES; i = next++) {
+                gemmi::Structure rep = placed;
+                std::vector<gemmi::Position> pos = ModelPositions(rep.models[0]);
+                gemmi::Vec3 c;
+                for (const gemmi::Position &p : pos)
+                    c += p;
+                c *= 1.0 / static_cast<double>(pos.size());
+                for (gemmi::Position &p : pos)
+                    p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + c);
+                SetModelPositions(rep.models[0], pos);
+                r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads).r_work;
+            }
+        }));
+    for (auto &w : workers)
+        w.get();
+    const double mean = std::accumulate(r_work.begin(), r_work.end(), 0.0) / NULL_REPLICATES;
+    double s2 = 0;
+    for (double x : r_work)
+        s2 += (x - mean) * (x - mean);
+    const double sd = std::sqrt(s2 / (NULL_REPLICATES - 1));
+    return sd > 0 ? (mean - real) / sd : 0.0;
+}
+
+bool same_rotations_and_centring(const gemmi::GroupOps &a, const gemmi::GroupOps &b) {
+    if (a.sym_ops.size() != b.sym_ops.size() || a.cen_ops.size() != b.cen_ops.size())
+        return false;
+    for (const gemmi::Op &op : a.sym_ops)
+        if (b.find_by_rotation(op.rot) == nullptr)
+            return false;
+    for (const gemmi::Op::Tran &c : a.cen_ops)
+        if (std::find(b.cen_ops.begin(), b.cen_ops.end(), c) == b.cen_ops.end())
+            return false;
+    return true;
+}
+
+bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
+    if (!a.is_enantiomorphic())
+        return false;
+    gemmi::GroupOps e = a.operations();
+    e.change_basis_forward(a.change_of_hand_op());
+    const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e);
+    return other != nullptr && other->number == b.number;
+}
+
+}  // namespace
+
+// --- the model ---
+
+SearchModel PrepareSearchModel(const std::string &path) {
+    SearchModel out;
+    try {
+        out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
+    } catch (const std::exception &e) {
+        out.failure_reason = std::string("cannot read the model: ") + e.what();
+        return out;
+    }
+    if (out.st.models.empty()) {
+        out.failure_reason = "the model has no atoms";
+        return out;
+    }
+    out.st.models.resize(1);
+    gemmi::remove_hydrogens(out.st);
+    gemmi::Model &model = out.st.models[0];
+    for (gemmi::Chain &ch : model.chains) {
+        for (const gemmi::Residue &r : ch.residues)
+            for (const gemmi::Atom &a : r.atoms)
+                if (a.b_iso > 100.0f) {
+                    out.failure_reason = "the B-factor column holds values above 100, so it is not a pLDDT: "
+                                         "not a predicted model";
+                    return out;
+                }
+        out.residues_read += static_cast<int>(ch.residues.size());
+        // The residues a prediction is confident of: pLDDT 70 and above, as Oeffner et al. (2022)
+        // Acta Cryst. D78, 1303-1314 trim an AlphaFold model for molecular replacement.
+        std::erase_if(ch.residues, [](const gemmi::Residue &r) { return r.atoms.empty() || r.atoms[0].b_iso < 70.0f; });
+        out.residues_kept += static_cast<int>(ch.residues.size());
+    }
+    std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); });
+    if (out.residues_kept == 0) {
+        out.failure_reason = "no residue of the model has a pLDDT of 70 or more";
+        return out;
+    }
+    // B from pLDDT through the coordinate error it predicts: rmsd = 1.5 exp(4 (0.7 - pLDDT/100)) A
+    // and B = 8 pi^2 rmsd^2 / 3. Following Oeffner et al. (2022) Acta Cryst. D78, 1303-1314
+    std::vector<gemmi::Position> pos;
+    for (gemmi::Chain &ch : model.chains)
+        for (gemmi::Residue &r : ch.residues)
+            for (gemmi::Atom &a : r.atoms) {
+                const double rmsd = 1.5 * std::exp(4.0 * (0.7 - a.b_iso / 100.0));
+                a.b_iso = static_cast<float>(8.0 * PI * PI * rmsd * rmsd / 3.0);
+                pos.push_back(a.pos);
+            }
+    gemmi::Vec3 c;
+    for (const gemmi::Position &p : pos)
+        c += p;
+    c *= 1.0 / static_cast<double>(pos.size());
+    for (gemmi::Position &p : pos) {
+        p = gemmi::Position(gemmi::Vec3(p) - c);
+        out.radius_A = std::max(out.radius_A, p.length());
+    }
+    SetModelPositions(model, pos);
+    out.mass_da = DALTON_PER_RESIDUE * out.residues_kept;
+    return out;
+}
+
+// --- rotations ---
+
+std::vector<gemmi::Mat33> UniformRotations(double step_deg) {
+    const double step = step_deg * PI / 180.0;
+    const double golden = PI * (3.0 - std::sqrt(5.0));
+    const int shells = std::max(1, static_cast<int>(std::lround(PI / step)));
+    std::vector<gemmi::Mat33> out{gemmi::Mat33()};
+    for (int i = 1; i <= shells; i++) {
+        const double kappa = PI * i / shells;
+        const double reach = step / (2.0 * std::sin(kappa / 2.0));   // the axis spacing this shell needs
+        int n = std::max(1, static_cast<int>(std::lround(4.0 * PI / (reach * reach))));
+        // At 180 deg an axis and its opposite are the same rotation: one hemisphere.
+        const bool half = i == shells;
+        if (half)
+            n = (n + 1) / 2;
+        for (int k = 0; k < n; k++) {
+            const double z = half ? 1.0 - (k + 0.5) / n : 1.0 - (2.0 * k + 1.0) / n;
+            const double r = std::sqrt(std::max(0.0, 1.0 - z * z));
+            const double x = r * std::cos(golden * k), y = r * std::sin(golden * k);
+            const double s = std::sin(kappa), c1 = 1.0 - std::cos(kappa);
+            out.push_back({1 + c1 * (x * x - 1), -s * z + c1 * x * y, s * y + c1 * x * z,
+                           s * z + c1 * x * y, 1 + c1 * (y * y - 1), -s * x + c1 * y * z,
+                           -s * y + c1 * x * z, s * x + c1 * y * z, 1 + c1 * (z * z - 1)});
+        }
+    }
+    return out;
+}
+
+std::vector<gemmi::Mat33> CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) {
+    std::vector<gemmi::Mat33> out;
+    for (const gemmi::Op &op : sg.operations().sym_ops) {
+        gemmi::Mat33 f;
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++)
+                f.a[i][j] = static_cast<double>(op.rot[i][j]) / gemmi::Op::DEN;
+        out.push_back(cell.orth.mat.multiply(f).multiply(cell.frac.mat));
+    }
+    return out;
+}
+
+std::vector<gemmi::Mat33> RotationAsu(const std::vector<gemmi::Mat33> &all,
+                                      const std::vector<gemmi::Mat33> &point_group) {
+    std::vector<gemmi::Mat33> out;
+    for (const gemmi::Mat33 &r : all) {
+        const double t = trace(r);
+        if (std::all_of(point_group.begin(), point_group.end(),
+                        [&](const gemmi::Mat33 &s) { return t >= trace(s.multiply(r)) - 1e-9; }))
+            out.push_back(r);
+    }
+    return out;
+}
+
+double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b,
+                           const std::vector<gemmi::Mat33> &point_group) {
+    double best = -1.0;
+    for (const gemmi::Mat33 &s : point_group)
+        best = std::max(best, trace(a.multiply(s.multiply(b).transpose())));
+    return std::acos(std::clamp((best - 1.0) / 2.0, -1.0, 1.0)) * 180.0 / PI;
+}
+
+// --- the molecular transform ---
+
+MolecularTransform::MolecularTransform(const SearchModel &model, double d_min, size_t nthreads) {
+    // On a P1 box six times the model's radius: the transform of a body of radius r turns its phase
+    // by pi r / box per grid step, and trilinear interpolation reads it to a few percent at pi/6.
+    // Interpolating the molecular transform of one model in place of a Fourier summation per
+    // orientation, following Navaza (1994) Acta Cryst. A50, 157-163
+    box_ = std::max(6.0 * model.radius_A, 120.0);
+    gemmi::Structure st = model.st;
+    st.cell = gemmi::UnitCell(box_, box_, box_, 90, 90, 90);
+    st.spacegroup_hm = "P 1";
+    st.setup_cell_images();
+    gemmi::DensityCalculator<Table, float> dc;
+    dc.d_min = d_min;
+    dc.rate = 1.5;
+    dc.set_grid_cell_and_spacegroup(st);
+    dc.set_refmac_compatible_blur(st.models[0]);
+    PutModelDensityOnGrid(dc, st.models[0], {}, nthreads);
+    grid_ = MapToFPhi(dc.grid);
+    blur_ = dc.blur;
+}
+
+std::complex<double> MolecularTransform::At(const gemmi::Vec3 &s) const {
+    const double u[3] = {s.x * box_, s.y * box_, s.z * box_};
+    const int i0[3] = {static_cast<int>(std::floor(u[0])), static_cast<int>(std::floor(u[1])),
+                       static_cast<int>(std::floor(u[2]))};
+    const double f[3] = {u[0] - i0[0], u[1] - i0[1], u[2] - i0[2]};
+    const auto get = [&](int i, int j, int l) {
+        bool conj = false;
+        if (l < 0) {   // the half-l grid holds l >= 0; F(-h) = conj F(h)
+            i = -i; j = -j; l = -l;
+            conj = true;
+        }
+        i %= grid_.nu; if (i < 0) i += grid_.nu;
+        j %= grid_.nv; if (j < 0) j += grid_.nv;
+        const std::complex<float> v = grid_.data[i + static_cast<size_t>(grid_.nu) * (j + static_cast<size_t>(grid_.nv) * l)];
+        return conj ? std::conj(std::complex<double>(v)) : std::complex<double>(v);
+    };
+    std::complex<double> out = 0;
+    for (int dx = 0; dx < 2; dx++)
+        for (int dy = 0; dy < 2; dy++)
+            for (int dz = 0; dz < 2; dz++) {
+                const double w = (dx ? f[0] : 1 - f[0]) * (dy ? f[1] : 1 - f[1]) * (dz ? f[2] : 1 - f[2]);
+                out += w * get(i0[0] + dx, i0[1] + dy, i0[2] + dz);
+            }
+    return out * std::exp(blur_ * s.length_sq() / 4.0);
+}
+
+// --- the observations ---
+
+SearchReflections MakeSearchReflections(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
+                                        const gemmi::SpaceGroup &sg, double d_min, double d_max) {
+    const gemmi::GroupOps gops = sg.operations();
+    const gemmi::ReciprocalAsu asu(&sg);
+    struct Row { gemmi::Miller h; double i_eps; double s2; double mult; };
+    std::vector<Row> rows;
+    for (const MergedReflection &r : merged) {
+        if (!std::isfinite(r.I) || r.d < d_min || r.d > d_max)
+            continue;
+        gemmi::Miller h{{r.h, r.k, r.l}};
+        if (!asu.is_in(h))
+            h = asu.to_asu(h, gops).first;
+        if (gops.is_systematically_absent(h))
+            continue;
+        // The reflections the whole sphere holds for this one: its images under the point group and
+        // their Friedel mates, each counted once.
+        std::vector<gemmi::Miller> images;
+        for (const gemmi::Op &op : gops.sym_ops)
+            for (int sign : {1, -1}) {
+                gemmi::Miller k = rotated_index(op, h);
+                for (int &x : k)
+                    x *= sign;
+                if (std::find(images.begin(), images.end(), k) == images.end())
+                    images.push_back(k);
+            }
+        rows.push_back({h, r.I / gops.epsilon_factor_without_centering(h), cell.calculate_1_d2(h),
+                        static_cast<double>(images.size())});
+    }
+    std::stable_sort(rows.begin(), rows.end(), [](const Row &a, const Row &b) { return a.s2 < b.s2; });
+    SearchReflections out;
+    const size_t n = rows.size();
+    out.bins = static_cast<int>(std::min<size_t>(SEARCH_BINS, std::max<size_t>(n / 100, 1)));
+    std::vector<double> sum(out.bins, 0.0);
+    std::vector<int> count(out.bins, 0);
+    for (size_t i = 0; i < n; i++) {
+        const int b = static_cast<int>(i * out.bins / n);
+        out.hkl.push_back(rows[i].h);
+        out.multiplicity.push_back(rows[i].mult);
+        out.bin.push_back(b);
+        sum[b] += rows[i].i_eps;
+        count[b]++;
+    }
+    for (size_t i = 0; i < n; i++) {
+        const double mean = sum[out.bin[i]] / count[out.bin[i]];
+        out.e2.push_back(mean > 0 ? rows[i].i_eps / mean : 0.0);
+    }
+    return out;
+}
+
+// --- the rotation function ---
+
+std::vector<double> RotationFunction(const MolecularTransform &ft, const SearchReflections &refl,
+                                     const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                     const std::vector<gemmi::Mat33> &rotations, size_t nthreads) {
+    const gemmi::GroupOps gops = sg.operations();
+    const size_t nops = gops.sym_ops.size(), n = refl.hkl.size();
+    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
+    // Copy j of the model contributes |F_model(s)|^2 at s = (h R_j) as a Cartesian reciprocal vector.
+    std::vector<gemmi::Vec3> s(n * nops);
+    for (size_t i = 0; i < n; i++)
+        for (size_t j = 0; j < nops; j++)
+            s[i * nops + j] = frac_t.multiply(miller_vec(rotated_index(gops.sym_ops[j], refl.hkl[i])));
+    double w_mean = 0;
+    for (double e2 : refl.e2)
+        w_mean += e2;
+    w_mean /= static_cast<double>(n);
+    double w_ss = 0;
+    for (double e2 : refl.e2)
+        w_ss += (e2 - w_mean) * (e2 - w_mean);
+
+    std::vector<double> out(rotations.size(), 0.0);
+    ParallelFor(static_cast<int>(rotations.size()), nthreads, [&](int r) {
+        // The model turned by R has F_R(s) = F_model(R^T s).
+        const gemmi::Mat33 rt = rotations[r].transpose();
+        std::vector<double> ic(n, 0.0), bin_sum(refl.bins, 0.0);
+        std::vector<int> bin_n(refl.bins, 0);
+        for (size_t i = 0; i < n; i++) {
+            for (size_t j = 0; j < nops; j++)
+                ic[i] += std::norm(ft.At(rt.multiply(s[i * nops + j])));
+            bin_sum[refl.bin[i]] += ic[i];
+            bin_n[refl.bin[i]]++;
+        }
+        double x_mean = 0;
+        for (size_t i = 0; i < n; i++) {
+            ic[i] /= bin_sum[refl.bin[i]] / bin_n[refl.bin[i]];
+            x_mean += ic[i];
+        }
+        x_mean /= static_cast<double>(n);
+        double sxy = 0, sxx = 0;
+        for (size_t i = 0; i < n; i++) {
+            sxy += (ic[i] - x_mean) * (refl.e2[i] - w_mean);
+            sxx += (ic[i] - x_mean) * (ic[i] - x_mean);
+        }
+        out[r] = sxx > 0 && w_ss > 0 ? sxy / std::sqrt(sxx * w_ss) : 0.0;
+    });
+    return out;
+}
+
+// --- the translation function ---
+
+TranslationPeak TranslationFunction(const MolecularTransform &ft, const SearchReflections &refl,
+                                    const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                    const gemmi::Mat33 &rotation) {
+    // The full-symmetry translation function as one Fourier synthesis, following Crowther & Blow
+    // (1967) Acta Cryst. 23, 544-548: with the model at t, copy j adds G_j(h) exp(2 pi i (h R_j).t)
+    // exp(2 pi i h.t_j) to F(h), so |F(h)|^2 summed against E^2 - 1 over the sphere is, for each pair
+    // of copies j != k, a Fourier series in t with coefficient (E^2 - 1) G_j conj(G_k) exp(2 pi i
+    // h.(t_j - t_k)) at the index h R_j - h R_k. Every pair's terms go on one grid, and one FFT gives
+    // the function at every t. The pairs j = k are the rotation function's, constant in t.
+    const gemmi::GroupOps gops = sg.operations();
+    const size_t nops = gops.sym_ops.size();
+    TranslationPeak out;
+    if (nops < 2)
+        return out;   // P1: every position is the same
+    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
+    const gemmi::Mat33 rt = rotation.transpose();
+    // Per copy j: k = h R_j, and g = G_j(h) exp(2 pi i h.t_j), so a pair's coefficient is g_j conj(g_k).
+    struct Term { gemmi::Miller h; std::vector<gemmi::Miller> k; std::vector<std::complex<double>> g; double w; };
+    std::vector<Term> terms;
+    std::vector<double> bin_sum(refl.bins, 0.0);
+    std::vector<int> bin_n(refl.bins, 0);
+    std::vector<int> term_bin;
+    for (size_t i = 0; i < refl.hkl.size(); i++) {
+        if (gops.is_systematically_absent(refl.hkl[i]))
+            continue;
+        Term t{refl.hkl[i], {}, {}, refl.multiplicity[i] * (refl.e2[i] - 1.0)};
+        double sum = 0;
+        for (const gemmi::Op &op : gops.sym_ops) {
+            t.k.push_back(rotated_index(op, refl.hkl[i]));
+            const double ht = (t.h[0] * op.tran[0] + t.h[1] * op.tran[1] + t.h[2] * op.tran[2])
+                              / static_cast<double>(gemmi::Op::DEN);
+            t.g.push_back(ft.At(rt.multiply(frac_t.multiply(miller_vec(t.k.back())))) * std::polar(1.0, 2.0 * PI * ht));
+            sum += std::norm(t.g.back());
+        }
+        bin_sum[refl.bin[i]] += sum;
+        bin_n[refl.bin[i]]++;
+        term_bin.push_back(refl.bin[i]);
+        terms.push_back(std::move(t));
+    }
+    if (terms.empty())
+        return out;
+    int max_m[3] = {0, 0, 0};
+    for (const Term &t : terms)
+        for (size_t j = 0; j < nops; j++)
+            for (size_t k = 0; k < nops; k++)
+                for (int a = 0; a < 3; a++)
+                    max_m[a] = std::max(max_m[a], std::abs(t.k[j][a] - t.k[k][a]));
+    const int nu = good_even_size(2 * max_m[0] + 2), nv = good_even_size(2 * max_m[1] + 2),
+              nw = good_even_size(2 * max_m[2] + 2);
+    const int nwh = nw / 2 + 1;
+    std::vector<std::complex<double>> c(static_cast<size_t>(nu) * nv * nwh, 0.0);
+    for (size_t i = 0; i < terms.size(); i++) {
+        const Term &t = terms[i];
+        const double norm = bin_sum[term_bin[i]] / bin_n[term_bin[i]];
+        for (size_t j = 0; j < nops; j++)
+            for (size_t k = 0; k < nops; k++) {
+                if (j == k)
+                    continue;
+                const int m[3] = {t.k[j][0] - t.k[k][0], t.k[j][1] - t.k[k][1], t.k[j][2] - t.k[k][2]};
+                // Only l >= 0 is stored: the pair (k, j) carries the conjugate term at -m.
+                if (m[2] < 0)
+                    continue;
+                const std::complex<double> v = t.w / norm * t.g[j] * std::conj(t.g[k]);
+                const int u = (m[0] % nu + nu) % nu, vv = (m[1] % nv + nv) % nv;
+                c[u + static_cast<size_t>(nu) * (vv + static_cast<size_t>(nv) * m[2])] += v;
+            }
+    }
+    // MapFromFPhi sums F exp(-2 pi i m.x): with F = conj(C) that is the conjugate of the series, which
+    // is real.
+    gemmi::FPhiGrid<float> hkl;
+    hkl.unit_cell = cell;
+    hkl.spacegroup = &gemmi::get_spacegroup_p1();
+    hkl.axis_order = gemmi::AxisOrder::XYZ;
+    hkl.half_l = true;
+    hkl.set_size_without_checking(nu, nv, nwh);
+    for (size_t i = 0; i < c.size(); i++)
+        hkl.data[i] = std::complex<float>(std::conj(c[i]));
+    const gemmi::Grid<float> map = MapFromFPhi(hkl);
+    double sum = 0, sum2 = 0;
+    size_t best = 0;
+    for (size_t i = 0; i < map.data.size(); i++) {
+        sum += map.data[i];
+        sum2 += static_cast<double>(map.data[i]) * map.data[i];
+        if (map.data[i] > map.data[best])
+            best = i;
+    }
+    const double count = static_cast<double>(map.data.size());
+    const double mean = sum / count, sd = std::sqrt(std::max(0.0, sum2 / count - mean * mean));
+    out.z = sd > 0 ? (map.data[best] - mean) / sd : 0.0;
+    const size_t u = best % map.nu, v = (best / map.nu) % map.nv, w = best / (static_cast<size_t>(map.nu) * map.nv);
+    out.t = gemmi::Fractional(static_cast<double>(u) / map.nu, static_cast<double>(v) / map.nv,
+                              static_cast<double>(w) / map.nw);
+    return out;
+}
+
+// --- placements ---
+
+gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
+                                  const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg) {
+    gemmi::Structure st = model.st;
+    const gemmi::Position shift = cell.orthogonalize(t);
+    std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
+    for (gemmi::Position &p : pos)
+        p = gemmi::Position(rotation.multiply(p) + shift);
+    SetModelPositions(st.models[0], pos);
+    st.cell = cell;
+    st.spacegroup_hm = sg.xhm();
+    st.setup_cell_images();
+    return st;
+}
+
+double ClashFraction(const gemmi::Structure &placed) {
+    const gemmi::UnitCell &cell = placed.cell;
+    const gemmi::SpaceGroup *sg = placed.find_spacegroup();
+    std::vector<gemmi::Fractional> ca;
+    for (const gemmi::Chain &ch : placed.models[0].chains)
+        for (const gemmi::Residue &r : ch.residues)
+            for (const gemmi::Atom &a : r.atoms)
+                if (a.name == "CA")
+                    ca.push_back(cell.fractionalize(a.pos));
+    if (ca.empty() || sg == nullptr)
+        return 0.0;
+    gemmi::Fractional c0;
+    for (const gemmi::Fractional &f : ca)
+        c0 = c0 + f;
+    c0 = gemmi::Fractional(c0 * (1.0 / static_cast<double>(ca.size())));
+    double radius = 0;
+    for (const gemmi::Fractional &f : ca)
+        radius = std::max(radius, cell.orthogonalize_difference(f - c0).length());
+    constexpr double CLASH_A = 3.0;
+    std::vector<char> clash(ca.size(), 0);
+    for (const gemmi::Op &op : sg->operations().all_ops_sorted()) {
+        std::vector<gemmi::Fractional> mate;
+        for (const gemmi::Fractional &f : ca) {
+            const std::array<double, 3> x = op.apply_to_xyz({{f.x, f.y, f.z}});
+            mate.push_back(gemmi::Fractional(x[0], x[1], x[2]));
+        }
+        gemmi::Fractional cm;
+        for (const gemmi::Fractional &f : mate)
+            cm = cm + f;
+        cm = gemmi::Fractional(cm * (1.0 / static_cast<double>(mate.size())));
+        const gemmi::Fractional nearest(std::round(c0.x - cm.x), std::round(c0.y - cm.y), std::round(c0.z - cm.z));
+        for (int dx = -1; dx <= 1; dx++)
+            for (int dy = -1; dy <= 1; dy++)
+                for (int dz = -1; dz <= 1; dz++) {
+                    const gemmi::Fractional shift(nearest.x + dx, nearest.y + dy, nearest.z + dz);
+                    if (op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0)
+                        continue;
+                    if (cell.orthogonalize_difference(cm + shift - c0).length() > 2 * radius + CLASH_A)
+                        continue;
+                    for (size_t a = 0; a < ca.size(); a++) {
+                        if (clash[a])
+                            continue;
+                        for (const gemmi::Fractional &m : mate)
+                            if (cell.orthogonalize_difference(m + shift - ca[a]).length_sq() < CLASH_A * CLASH_A) {
+                                clash[a] = 1;
+                                break;
+                            }
+                    }
+                }
+    }
+    return static_cast<double>(std::count(clash.begin(), clash.end(), 1)) / static_cast<double>(ca.size());
+}
+
+std::vector<gemmi::SpaceGroup> SameRotationGroups(const gemmi::SpaceGroup &sg) {
+    const gemmi::GroupOps ops = sg.operations();
+    std::vector<gemmi::SpaceGroup> out{sg};
+    for (const gemmi::SpaceGroup &g : gemmi::spacegroup_tables::main) {
+        if (g.number == sg.number && g.xhm() == sg.xhm())
+            continue;
+        const gemmi::GroupOps gops = g.operations();
+        if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g)
+            continue;
+        out.push_back(g);
+    }
+    return out;
+}
+
+// --- the search and the decision ---
+
+MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
+                                                   Logger &logger, size_t nthreads) {
+    MolecularReplacementResult result;
+    result.model_path = model_path;
+    auto t0 = std::chrono::steady_clock::now();
+    const SearchModel model = PrepareSearchModel(model_path);
+    result.residues_read = model.residues_read;
+    result.residues_kept = model.residues_kept;
+    result.mass_da = model.mass_da;
+    if (!model.failure_reason.empty()) {
+        result.failure_reason = model.failure_reason;
+        return result;
+    }
+    if (merges.empty() || merges[0].merged == nullptr) {
+        result.failure_reason = "no merge to search";
+        return result;
+    }
+    double data_d_min = 1e9;
+    for (const MergedReflection &r : *merges[0].merged)
+        if (r.d > 0)
+            data_d_min = std::min<double>(data_d_min, r.d);
+    result.d_min_search = std::max(SEARCH_D_MIN, data_d_min);
+    result.d_max_search = SEARCH_D_MAX;
+    result.d_min_score = std::max(SCORE_D_MIN, data_d_min);
+    result.t_model = seconds_since(t0);
+    logger.Info("Model prior: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)",
+                model.residues_kept, model.residues_read, model.mass_da / 1000.0, model.radius_A);
+
+    t0 = std::chrono::steady_clock::now();
+    const MolecularTransform ft(model, result.d_min_search, nthreads);
+    result.t_transform = seconds_since(t0);
+    const std::vector<gemmi::Mat33> sphere = UniformRotations(ROTATION_STEP_DEG);
+
+    // One merge: the rotation function over its point group, then for each of its groups the
+    // translation function and the best placements through packing, the rigid body and R-free.
+    const auto search_merge = [&](size_t m) {
+        const MRMerge &merge = merges[m];
+        const gemmi::UnitCell cell(merge.cell);
+        const double d_min = result.d_min_score;
+        // --- rotation ---
+        t0 = std::chrono::steady_clock::now();
+        const SearchReflections refl = MakeSearchReflections(*merge.merged, cell, merge.space_group,
+                                                             result.d_min_search, result.d_max_search);
+        const std::vector<gemmi::Mat33> pg = CartesianPointGroup(merge.space_group, cell);
+        const std::vector<gemmi::Mat33> rotations = RotationAsu(sphere, pg);
+        const std::vector<double> rf = RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
+        double rf_mean = 0, rf_sd = 0;
+        for (double x : rf)
+            rf_mean += x;
+        rf_mean /= static_cast<double>(rf.size());
+        for (double x : rf)
+            rf_sd += (x - rf_mean) * (x - rf_mean);
+        rf_sd = std::sqrt(rf_sd / static_cast<double>(rf.size()));
+        // The highest orientations, one per peak: a rotation within two steps of one already taken
+        // (in any of its symmetry classes) belongs to that one's peak.
+        std::vector<size_t> order(rf.size());
+        std::iota(order.begin(), order.end(), 0);
+        std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) { return rf[a] > rf[b]; });
+        std::vector<size_t> peaks;
+        for (size_t i : order) {
+            if (std::none_of(peaks.begin(), peaks.end(), [&](size_t p) {
+                    return RotationDistanceDeg(rotations[i], rotations[p], pg) < 2 * ROTATION_STEP_DEG;
+                }))
+                peaks.push_back(i);
+            if (static_cast<int>(peaks.size()) == TRANSLATION_ROTATIONS)
+                break;
+        }
+        result.t_rotation += seconds_since(t0);
+        logger.Info("Model prior: {} in {}: {} reflections {:.1f}-{:.1f} A, {} orientations, best RF Z {:.1f}",
+                    merge.name, merge.space_group.xhm(), refl.hkl.size(), result.d_max_search, result.d_min_search,
+                    rotations.size(), rf_sd > 0 ? (rf[peaks.front()] - rf_mean) / rf_sd : 0.0);
+
+        // --- translation, for every group and every peak orientation ---
+        t0 = std::chrono::steady_clock::now();
+        const std::vector<gemmi::SpaceGroup> groups = SameRotationGroups(merge.space_group);
+        std::vector<TranslationPeak> tf_all(groups.size() * peaks.size());
+        ParallelFor(static_cast<int>(tf_all.size()), nthreads, [&](int i) {
+            tf_all[i] = TranslationFunction(ft, refl, cell, groups[i / peaks.size()], rotations[peaks[i % peaks.size()]]);
+        });
+        result.t_translation += seconds_since(t0);
+
+        for (size_t g = 0; g < groups.size(); g++) {
+            const gemmi::SpaceGroup &sg = groups[g];
+            MRCandidate cand;
+            cand.merge = static_cast<int>(m);
+            cand.space_group = sg;
+            const gemmi::GroupOps gops = sg.operations();
+            cand.vm = cell.volume / static_cast<double>(gops.order()) / model.mass_da;
+            cand.copies_likely = cand.vm >= VM_TWO_COPIES ? 2 : 1;
+            const TranslationPeak *tf = &tf_all[g * peaks.size()];
+            std::vector<size_t> by_tf(peaks.size());
+            std::iota(by_tf.begin(), by_tf.end(), 0);
+            std::stable_sort(by_tf.begin(), by_tf.end(), [&](size_t a, size_t b) { return tf[a].z > tf[b].z; });
+
+            // --- the best few through packing, the rigid body and R-free ---
+            t0 = std::chrono::steady_clock::now();
+            RigidBodyEngines engines;
+            bool have = false;
+            gemmi::Mat33 orientation;
+            // Where one copy cannot be the content of the cell - too big for it, or small enough for two -
+            // the group is not decided on here, and only its best translation-function placement is
+            // checked for packing.
+            const bool one_copy = cand.vm >= VM_MIN && cand.copies_likely == 1;
+            for (int k = 0; k < std::min<int>(one_copy ? RESCORED_PER_GROUP : 1, static_cast<int>(by_tf.size())); k++) {
+                const size_t i = by_tf[k];
+                auto placed = std::make_shared<gemmi::Structure>(
+                        PlaceSearchModel(model, rotations[peaks[i]], tf[i].t, cell, sg));
+                PlacementScore score{NAN, NAN};
+                if (one_copy) {
+                    if (k == 0)
+                        engines = make_engines(*placed, sg, *merge.merged, d_min, NULL_REPLICATES + 1, logger);
+                    score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads);
+                }
+                const double clash = ClashFraction(*placed);
+                if (!have || score.r_free < cand.r_free) {   // NaN (not refined) only ever on k = 0
+                    have = true;
+                    cand.rf_z = rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0;
+                    cand.tf_z = tf[i].z;
+                    cand.clash_fraction = clash;
+                    cand.r_work = score.r_work;
+                    cand.r_free = score.r_free;
+                    cand.placed = placed;
+                    orientation = rotations[peaks[i]];
+                }
+                logger.Info("Model prior:   {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
+                            sg.xhm(), rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash,
+                            score.r_work, score.r_free);
+            }
+            result.t_rescore += seconds_since(t0);
+
+            if (cand.vm < VM_MIN)
+                cand.verdict = "MATTHEWS";
+            else if (cand.copies_likely > 1)
+                cand.verdict = "MULTI_COPY";
+            else if (cand.clash_fraction > CLASH_MAX)
+                cand.verdict = "CLASH";
+            else if (!(cand.tf_z >= TF_Z_PLACED))
+                cand.verdict = "NOT_PLACED";
+            else
+                cand.verdict = "PLACED";
+            if (cand.verdict == "PLACED") {
+                t0 = std::chrono::steady_clock::now();
+                cand.null_sigma = null_sigma(*cand.placed, orientation, sg, *merge.merged, d_min, engines.get(),
+                                             nthreads);
+                result.t_null += seconds_since(t0);
+                if (!(cand.null_sigma >= NULL_SIGMA_MIN))
+                    cand.verdict = "NULL";
+            }
+            result.candidates.push_back(std::move(cand));
+        }
+    };
+
+    // --- the decision ---
+    // Within one merge the candidates are scored on the same reflections, so their R-free values can
+    // be compared directly: the placed candidate with the lowest R-free is chosen where it leads every
+    // other candidate of that merge that packs by R_FREE_MARGIN.
+    const auto choose_in = [&](int merge) -> int {
+        int best = -1;
+        for (size_t i = 0; i < result.candidates.size(); i++) {
+            const MRCandidate &c = result.candidates[i];
+            if (c.merge == merge && c.verdict == "PLACED" && (best < 0 || c.r_free < result.candidates[best].r_free))
+                best = static_cast<int>(i);
+        }
+        if (best < 0)
+            return -1;
+        for (size_t i = 0; i < result.candidates.size(); i++) {
+            const MRCandidate &c = result.candidates[i];
+            if (static_cast<int>(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX)
+                continue;
+            if (!(c.r_free - result.candidates[best].r_free >= R_FREE_MARGIN)) {
+                result.decision = fmt::format("UNDECIDED: {} and {} within {:.2f} in R-free",
+                                              result.candidates[best].space_group.xhm(), c.space_group.xhm(),
+                                              R_FREE_MARGIN);
+                return -2;
+            }
+        }
+        return best;
+    };
+    // A merge whose every group the model cannot be in at all - one copy does not fit the cell, or
+    // every placement clashes with its own symmetry mates.
+    const auto cannot_hold = [&](int merge) {
+        bool any = false;
+        for (const MRCandidate &c : result.candidates)
+            if (c.merge == merge) {
+                any = true;
+                if (c.verdict != "MATTHEWS" && c.verdict != "CLASH")
+                    return false;
+            }
+        return any;
+    };
+    // The twin hypothesis' merge is searched only where it can decide something: where the answer's
+    // merge cannot hold the model in any of its groups.
+    search_merge(0);
+    const int answer = choose_in(0);
+    if (answer >= 0) {
+        result.chosen = answer;
+    } else if (answer == -1 && merges.size() > 1 && cannot_hold(0)) {
+        search_merge(1);
+        const int twin = choose_in(1);
+        if (twin >= 0)
+            result.chosen = twin;
+        else if (twin == -1)
+            result.decision = "NOT_PLACED: no group of either merge holds a placed model";
+    } else if (answer == -1) {
+        result.decision = "NOT_PLACED: no group of the answer's merge holds a placed model";
+    }
+    if (result.chosen >= 0) {
+        const MRCandidate &c = result.candidates[result.chosen];
+        const gemmi::SpaceGroup &data = merges[0].space_group;
+        if (c.merge == 0 && c.space_group.xhm() == data.xhm())
+            result.decision = "AGREES";
+        else if (c.merge == 0 && enantiomorphs(data, c.space_group))
+            result.decision = "ENANTIOMORPH";
+        else if (c.merge == 0)
+            result.decision = "SPACE_GROUP";
+        else
+            result.decision = "TWIN_SUBGROUP";
+    }
+    result.ok = true;
+    logger.Info("Model prior: {} (model {:.1f} s, transform {:.2f} s, rotation {:.2f} s, translation {:.2f} s, "
+                "rescoring {:.2f} s, null {:.2f} s)",
+                result.chosen >= 0 ? result.decision + " - " + result.candidates[result.chosen].space_group.xhm()
+                                   : result.decision,
+                result.t_model, result.t_transform, result.t_rotation, result.t_translation, result.t_rescore,
+                result.t_null);
+    return result;
+}
@@ -0,0 +1,165 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+// Molecular replacement of ONE copy of a predicted model (rugnux --af-model, experimental): a global
+// search for the model's orientation and position in the cell, so that a model with no placement of
+// its own - an AlphaFold prediction - can be scored against the data, and the space groups the data
+// cannot tell apart can be scored against it.
+//
+//  - The model: residues with pLDDT < 70 removed, B-factors from pLDDT, hydrogens dropped.
+//  - Rotation: a uniform sampling of the point group's rotation ASU, each orientation scored by the
+//    correlation of the observed E^2 with the model's own symmetry-summed |F|^2 - read off the
+//    model's molecular transform, computed once by FFT and interpolated per orientation.
+//  - Translation: the full-symmetry translation function on one FFT per orientation and space group.
+//  - Rescoring: the best placements through the rigid body, bulk solvent, R-free, a packing check
+//    against the symmetry mates and a null of the same model in random orientations.
+//
+// Everything up to the rescoring runs on the CPU, in both builds; the rigid body runs on the GPU
+// where there is one.
+
+#include <complex>
+#include <memory>
+#include <string>
+#include <vector>
+
+#include "gemmi/model.hpp"
+#include "gemmi/recgrid.hpp"
+#include "gemmi/symmetry.hpp"
+#include "gemmi/unitcell.hpp"
+
+#include "../../common/Reflection.h"   // MergedReflection
+#include "../../common/UnitCell.h"
+
+class Logger;
+
+// The predicted model as the search uses it: trimmed, B-factors from pLDDT, centroid at the origin.
+struct SearchModel {
+    gemmi::Structure st;
+    int residues_read = 0;
+    int residues_kept = 0;
+    double mass_da = 0;        // of the residues kept, 110 Da each
+    double radius_A = 0;       // the farthest atom from the centroid
+    std::string failure_reason;
+};
+SearchModel PrepareSearchModel(const std::string &path);
+
+// Rotations covering SO(3) with neighbours about `step_deg` apart: the identity, then shells of
+// rotation angle kappa = step, 2 step, ... 180 deg, each with its axes on a golden-spiral lattice of
+// the sphere as dense as 2 sin(kappa/2) - the distance a change of axis moves a rotation - needs.
+std::vector<gemmi::Mat33> UniformRotations(double step_deg);
+
+// The rotations of the space group's point group, in Cartesian coordinates of `cell`.
+std::vector<gemmi::Mat33> CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell);
+
+// The rotations of `all` in the rotation ASU of `point_group`: the crystal puts a copy of the model
+// in every orientation S R, so only one rotation of each class {S R} is worth scoring - here the one
+// closest to the identity (largest trace), the classes' Dirichlet region about it.
+std::vector<gemmi::Mat33> RotationAsu(const std::vector<gemmi::Mat33> &all,
+                                      const std::vector<gemmi::Mat33> &point_group);
+
+// The angle in degrees between a and b, the nearest over the classes {S b}.
+double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b,
+                           const std::vector<gemmi::Mat33> &point_group);
+
+// The structure factors of one copy of the model, at its centroid, at any reciprocal vector.
+class MolecularTransform {
+public:
+    MolecularTransform(const SearchModel &model, double d_min, size_t nthreads);
+    // F at the Cartesian reciprocal vector s (1/A), interpolated trilinearly.
+    std::complex<double> At(const gemmi::Vec3 &s) const;
+
+private:
+    gemmi::FPhiGrid<float> grid_;   // half-l, on a P1 box `box_` A wide
+    double box_ = 0;
+    double blur_ = 0;
+};
+
+// The observed reflections the search scores against: E^2 in resolution bins, with each unique
+// reflection's multiplicity on the whole sphere.
+struct SearchReflections {
+    std::vector<gemmi::Miller> hkl;
+    std::vector<double> e2;
+    std::vector<double> multiplicity;
+    std::vector<int> bin;
+    int bins = 0;
+};
+SearchReflections MakeSearchReflections(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
+                                        const gemmi::SpaceGroup &sg, double d_min, double d_max);
+
+// The rotation function at each of `rotations`: the correlation of E^2 - 1 with the bin-normalised
+// sum over the point group's copies of |F_model|^2.
+std::vector<double> RotationFunction(const MolecularTransform &ft, const SearchReflections &refl,
+                                     const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                     const std::vector<gemmi::Mat33> &rotations, size_t nthreads);
+
+// The translation function of the model in orientation `rotation` in `sg`: its highest point as a
+// Z-score over every point of the function, and where it is (the fractional position of the
+// model's centroid).
+struct TranslationPeak {
+    double z = 0;
+    gemmi::Fractional t;
+};
+TranslationPeak TranslationFunction(const MolecularTransform &ft, const SearchReflections &refl,
+                                    const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg,
+                                    const gemmi::Mat33 &rotation);
+
+// The model placed: rotated about its centroid by `rotation`, the centroid at `t`, in `cell` and `sg`.
+gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t,
+                                  const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg);
+
+// The fraction of the model's C-alpha atoms within 3 A of a C-alpha of a symmetry mate.
+double ClashFraction(const gemmi::Structure &placed);
+
+// The space groups with the rotations and the centring of `sg`, in its setting: those merged
+// intensities cannot tell apart, apart from systematic absences. `sg` itself first.
+std::vector<gemmi::SpaceGroup> SameRotationGroups(const gemmi::SpaceGroup &sg);
+
+// One merge the model is searched against: the answer's, or the twin hypothesis' in a subgroup.
+struct MRMerge {
+    std::string name;                                // "MERGE" or "TWIN_SUBGROUP_MERGE"
+    const std::vector<MergedReflection> *merged = nullptr;
+    UnitCell cell;
+    gemmi::SpaceGroup space_group;                   // the group it was merged in
+};
+
+// One space group's best single-copy placement.
+struct MRCandidate {
+    int merge = 0;                                   // index into the merges searched
+    gemmi::SpaceGroup space_group;
+    double vm = 0;                                   // Matthews coefficient of one copy
+    int copies_likely = 1;
+    double rf_z = 0;                                 // of the orientation, over the rotation function
+    double tf_z = 0;
+    double clash_fraction = 0;
+    double r_work = NAN, r_free = NAN;               // after the rigid body, to d_min_score
+    double null_sigma = NAN;                         // R-work above its random-orientation null, in sd
+    std::string verdict;                             // PLACED, or why not
+    std::shared_ptr<gemmi::Structure> placed;
+};
+
+struct MolecularReplacementResult {
+    bool ok = false;
+    std::string failure_reason;
+    std::string model_path;
+    int residues_read = 0, residues_kept = 0;
+    double mass_da = 0;
+    double d_min_search = 0, d_max_search = 0;
+    double d_min_score = 0;                          // the limit R-work / R-free are computed to
+    std::vector<MRCandidate> candidates;
+    // The decision: the candidate the model chooses, or -1 where it chooses none, and in words why.
+    int chosen = -1;
+    std::string decision;
+    // Seconds per step.
+    double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0;
+};
+
+// Search `merges` (the answer's first) for one copy of the predicted model at `model_path`, in every
+// space group with the rotations of each merge's group, and decide. A candidate is PLACED where its
+// translation function stands out (TF Z >= 8), it packs, one copy fits the cell, and its R-work beats
+// its own null. The model chooses among the answer merge's groups where one placed candidate leads
+// every other by an R-free margin; it chooses the twin hypothesis' merge only where none of the
+// answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can.
+MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector<MRMerge> &merges,
+                                                   Logger &logger, size_t nthreads);
@@ -2562,6 +2562,81 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
        doc.sections.push_back(std::move(s));
    }

+    // ------------------------------------------------- MODEL PRIOR (--af-model, experimental)
+    if (result.model_prior.has_value()) {
+        const auto &mp = *result.model_prior;
+        ReportSection s;
+        s.title = "MODEL PRIOR (EXPERIMENTAL)";
+        Add(s, KeyText("AF_MODEL_FILE", mp.model_path));
+        if (!mp.ok) {
+            Add(s, KeyEnum("MODEL_PRIOR", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
+            Add(s, KeyText("MODEL_PRIOR_REASON", mp.failure_reason));
+        } else {
+            Add(s, KeyEnum("MODEL_PRIOR", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"}));
+            Add(s, KeyInt("AF_RESIDUES_READ", mp.residues_read));
+            Add(s, KeyInt("AF_RESIDUES_KEPT", mp.residues_kept));
+            Add(s, KeyReal("SEARCH_D_MIN", mp.d_min_search, "{:.2f}"));
+            Add(s, KeyReal("SCORE_D_MIN", mp.d_min_score, "{:.2f}"));
+            Add(s, KeyText("MODEL_PRIOR_DECISION", mp.decision));
+            if (mp.chosen >= 0)
+                Add(s, KeyText("MODEL_PRIOR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm()));
+            if (!result.model_prior_file_prefix.empty())
+                Add(s, KeyText("MODEL_PRIOR_FILE_PREFIX", result.model_prior_file_prefix));
+            Add(s, KeyReal("MODEL_PRIOR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation
+                                                          + mp.t_rescore + mp.t_null, "{:.2f}"));
+            Add(s, Blank());
+            ReportEntry e;
+            e.kind = ReportEntry::Kind::Table;
+            e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "RF_Z", "TF_Z", "CLASH", "R_WORK", "R_FREE",
+                               "NULL_SIGMA", "VERDICT"};
+            e.table.text_header = "  MERGE                SPACE_GROUP     VM   RF_Z   TF_Z  CLASH  R_WORK  R_FREE  NULL_SIGMA  VERDICT\n"
+                                  "  -------------------  -----------  -----  -----  -----  -----  ------  ------  ----------  ----------";
+            for (const auto &c : mp.candidates) {
+                const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE";
+                const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); };
+                e.table.text_rows.push_back(fmt::format(
+                        "  {:<19s}  {:<11s}  {:5.2f}  {:5.1f}  {:5.1f}  {:4.1f}%  {:>6s}  {:>6s}  {:>10s}  {}",
+                        merge, c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, 100 * c.clash_fraction, r(c.r_work),
+                        r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"),
+                        c.verdict));
+                std::vector<ReportValue> row;
+                row.push_back(KeyText("", merge).value);
+                row.push_back(KeyText("", c.space_group.xhm()).value);
+                row.push_back(KeyReal("", c.vm, "{:.2f}").value);
+                row.push_back(KeyReal("", c.rf_z, "{:.1f}").value);
+                row.push_back(KeyReal("", c.tf_z, "{:.1f}").value);
+                row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value);
+                row.push_back(KeyReal("", c.r_work, "{:.3f}").value);
+                row.push_back(KeyReal("", c.r_free, "{:.3f}").value);
+                row.push_back(KeyReal("", c.null_sigma, "{:+.1f}").value);
+                row.push_back(KeyText("", c.verdict).value);
+                e.table.cells.push_back(std::move(row));
+            }
+            Add(s, std::move(e));
+            Add(s, Blank());
+            Add(s, Prose(fmt::format(
+                    " One copy of the predicted model, trimmed to pLDDT >= 70, was searched for in every space group\n"
+                    " each merge cannot tell apart (rotation function on a 6 deg sampling, translation function,\n"
+                    " {:.1f}-{:.1f} A), and the best placements taken through the rigid body; R_WORK and R_FREE are\n"
+                    " to SCORE_D_MIN. PLACED needs TF_Z >= 8,\n"
+                    " at most 5% of the C-alpha atoms clashing with symmetry mates, one copy fitting the cell\n"
+                    " (MATTHEWS: Vm below 1.5; MULTI_COPY: Vm of 4 or more, not supported yet) and an R-work at\n"
+                    " least 3 sd better than the same model refitted from random orientations. The model chooses\n"
+                    " within a merge where its placed group leads every other that packs by 0.03 in R-free, and\n"
+                    " the twin hypothesis' merge - searched only then - where no group of the answer's merge can\n"
+                    " hold it at all. A group one copy cannot fill is checked for packing only, not refined.\n"
+                    " The answer's files take the enantiomorph from it and nothing else; another group it\n"
+                    " chooses is written as {}.mtz/.cif with the placed model.",
+                    mp.d_max_search, mp.d_min_search,
+                    result.model_prior_file_prefix.empty() ? std::string("<prefix>_mr") : result.model_prior_file_prefix)));
+            Add(s, Prose(fmt::format(
+                    " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n"
+                    " R-free {:.2f}, null {:.2f}.", mp.t_model, mp.t_transform, mp.t_rotation, mp.t_translation,
+                    mp.t_rescore, mp.t_null), true));
+        }
+        doc.sections.push_back(std::move(s));
+    }
+
    if (result.cancelled)
        Warn(doc, PathologyCode::CANCELLED,
             fmt::format("Processing was cancelled after {} images - this report describes an "
@@ -32,6 +32,7 @@
 #include "../image_analysis/scale_merge/TranslationalNCS.h"   // TranslationalNCSResult
 #include "../image_analysis/indexing/HarmonicContamination.h"  // HarmonicContaminationResult
 #include "../image_analysis/structure_refinement/ModelValidation.h"                  // ModelValidationResult
+#include "../image_analysis/structure_refinement/MolecularReplacement.h"            // MolecularReplacementResult
 #include "../image_analysis/indexing/LeftoverLattices.h"                 // LeftoverLatticeCensus
 #include "ReportDocument.h"                   // ReportWarning, PathologyCode
 #include "../image_analysis/scale_merge/SearchSpaceGroup.h"   // SearchSpaceGroupResult
@@ -196,6 +197,12 @@ struct ProcessConfig {
    // Fo-Fc maps (--model).
    std::string model_path;

+    // Experimental (--af-model): a predicted model (AlphaFold mmCIF, pLDDT in the B column) placed by
+    // molecular replacement, one copy, in every space group the merge cannot tell apart and in the twin
+    // hypothesis' subgroup. It never changes the answer's files beyond the enantiomorph label; a
+    // different group it chooses is written beside them as <prefix>_mr.
+    std::string af_model_path;
+
    // --finalist-ledger: report-only symmetry evidence table (and the twin-immune zones), folded from
    // a merge of its own in P1 that only this flag makes.
    bool finalist_ledger = false;
@@ -562,6 +569,11 @@ struct ProcessResult {
    // that was asked to validate looking like one that was not. Unset when no model was asked for.
    std::optional<ModelValidationResult> model_validation;

+    // What --af-model found, when one was given; failure_reason says why where it could not run.
+    std::optional<MolecularReplacementResult> model_prior;
+    // The file prefix the model's choice was written under (<prefix>_mr), empty where nothing was.
+    std::string model_prior_file_prefix;
+
    // Space group used for the final re-scale/merge and written to the master file: determined by the
    // search when the user did not fix one, otherwise the fixed group. The whole group, so it keeps
    // its setting - a number only ever names the reference one.
@@ -82,6 +82,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config,
        }
        if (!config.model_path.empty())
            add("--model", config.model_path);
+        if (!config.af_model_path.empty())
+            add("--af-model", config.af_model_path);
    }
    if (config.start_image != 0)
        add("-s", std::to_string(config.start_image));
@@ -3972,6 +3972,63 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
    }
    prepare_twin_file();

+    // --af-model (experimental): the predicted model placed by molecular replacement, one copy, in the
+    // groups the answer's merge cannot tell apart and in the twin hypothesis' subgroup. It takes the
+    // enantiomorph label into the answer, as --model does, and nothing else: a different group it
+    // chooses goes to files of its own beside the answer's, <prefix>_mr, from the merge it was chosen on.
+    struct ModelPriorFile {
+        std::vector<MergedReflection> merged;
+        UnitCell cell;
+        gemmi::SpaceGroup space_group;
+        std::optional<gemmi::Op> law;
+        TwinLawFraction fraction;
+        std::shared_ptr<gemmi::Structure> placed;
+    };
+    std::optional<ModelPriorFile> model_prior_file;
+    if (result.consensus_cell && write_files && !config_.af_model_path.empty() && !geometry_prepass
+        && !superseded) {
+        phase("Placing the predicted model");
+        Step("Model prior");
+        TimingMark mark(logger, "model prior");
+        std::vector<MRMerge> merges{{"MERGE", &sm.merged, *result.consensus_cell, experiment_.GetSpaceGroupOrP1()}};
+        if (twin_file)
+            merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group});
+        MolecularReplacementResult mr = RunMolecularReplacement(config_.af_model_path, merges, logger,
+                                                                static_cast<size_t>(config_.nthreads));
+        if (!mr.ok)
+            result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED,
+                                      "The predicted model was not placed: " + mr.failure_reason});
+        if (mr.chosen >= 0) {
+            const MRCandidate &c = mr.candidates[mr.chosen];
+            if (mr.decision == "ENANTIOMORPH") {
+                // A change of label only (see ValidateAgainstModel): the two groups index identically.
+                experiment_.SetSpaceGroup(c.space_group);
+                result.space_group = c.space_group;
+                if (twin_file && twin_file->space_group.is_enantiomorphic()) {
+                    gemmi::GroupOps e = twin_file->space_group.operations();
+                    e.change_basis_forward(twin_file->space_group.change_of_hand_op());
+                    if (const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e))
+                        twin_file->space_group = *other;
+                }
+                logger.Info("Model prior: the answer is relabelled {}, the enantiomorph the placed model chose",
+                            c.space_group.xhm());
+            } else if (mr.decision == "SPACE_GROUP") {
+                ModelPriorFile f{sm.merged, *result.consensus_cell, c.space_group, std::nullopt, {}, c.placed};
+                const gemmi::GroupOps gops = c.space_group.operations();
+                std::erase_if(f.merged, [&](const MergedReflection &r) {
+                    return gops.is_systematically_absent({{r.h, r.k, r.l}});
+                });
+                model_prior_file = std::move(f);
+            } else if (mr.decision == "TWIN_SUBGROUP" && twin_file) {
+                model_prior_file = ModelPriorFile{twin_file->merged, twin_file->cell, c.space_group, twin_file->law,
+                                                  twin_file->fraction, c.placed};
+            }
+        }
+        if (model_prior_file)
+            result.model_prior_file_prefix = config_.output_prefix + "_mr";
+        result.model_prior = std::move(mr);
+    }
+
    // The unmerged MTZ (written with the other output files, see there) is built from here on,
    // beside the reflection files: the indexing and the group it is written in are final once the
    // model has been read, and nothing after this changes the integrated observations it reads. It
@@ -4036,6 +4093,33 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
            WriteMmcifReflections(twin_file->merged, twin_file->cell, twin_x, twin_merge->statistics, twin_em,
                                  law, config_.provenance, prefix + ".cif", static_cast<size_t>(std::max(1, config_.nthreads)));
        }
+        // The predicted model's choice (see model_prior_file): its reflections, from the merge it was
+        // chosen on and labelled with the group it chose, and the model as it placed it.
+        if (model_prior_file) {
+            const ModelPriorFile &f = *model_prior_file;
+            DiffractionExperiment mr_x = experiment_;
+            mr_x.SetSpaceGroup(f.space_group);
+            TwinningAnalysisResult law;
+            const MergeStatistics *statistics = &sm.statistics;
+            ErrorModelReport mr_em = em_report;
+            if (f.law) {
+                law.twin_law_operator = f.law->as_hkl().triplet();
+                law.twin_law_fraction = std::isfinite(f.fraction.fraction) ? f.fraction.fraction : 0.5;
+                statistics = &twin_merge->statistics;
+                mr_em = ErrorModelReport{
+                        twin_merge->isa > 0 && twin_merge->isa_resolved ? fmt::format("{:.2f}", twin_merge->isa) : "?",
+                        std::string(), twin_merge->error_model_a > 0 ? fmt::format("{:.3f}", twin_merge->error_model_a)
+                                                                     : std::string(),
+                        twin_merge->error_model_b > 0 ? fmt::format("{:.4e}", twin_merge->error_model_b)
+                                                      : std::string()};
+            }
+            const std::string prefix = result.model_prior_file_prefix;
+            WriteMtzReflections(f.merged, f.cell, mr_x, config_.provenance, prefix + ".mtz");
+            WriteMmcifReflections(f.merged, f.cell, mr_x, *statistics, mr_em, law, config_.provenance, prefix + ".cif",
+                                  static_cast<size_t>(std::max(1, config_.nthreads)));
+            if (f.placed)
+                WritePlacedModel(*f.placed, f.cell, f.space_group, prefix, logger);
+        }
        writing_mark.Stop();

        // FINALIST LEDGER (--finalist-ledger), report-only and deliberately behind a flag: these
@@ -202,6 +202,7 @@ void print_usage() {
    std::cout << "   -z, --reference <file>          Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl;
    std::cout << "   --reference-column <label>      Reference column to use, by MTZ label; an mmCIF is read under the labels gemmi cif2mtz gives it (IMEAN, FP, I(+), FC, F-model, ...) (default: auto - F-model, else IMEAN/I, else FP/FOBS/F)" << std::endl;
    std::cout << "   --model <file>                  After merging, validate vs this atomic model: R-free, sigma_A-weighted 2mFo-DFc/mFo-DFc maps, the model as placed (written beside them), and - where the merge kept the Bijvoet split - an anomalous difference map with the strongest anomalous scatterers named by the atom they sit on. PDB or mmCIF, gzipped or not, recognised by content rather than by extension" << std::endl;
+    std::cout << "   --af-model <file>               Experimental: place this predicted model (AlphaFold mmCIF, pLDDT in the B column; residues below pLDDT 70 removed) by molecular replacement, one copy in the asymmetric unit, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. Reported in the MODEL PRIOR section; it may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as <prefix>_mr.mtz/.cif with the placed model <prefix>_mr_model.pdb. --mode mx only" << std::endl;
    std::cout << "   --export-unmerged               Write <prefix>_unmerged.mtz, an unmerged MTZ (POINTLESS column layout) of the integrated observations, for aimless / pointless / careless. On by default, whenever there is an output prefix. On a rotation run the partials of each reflection are summed into one full, written at the batch its rocking curve is centred on, with the summed rocking-curve fraction in FRACTIONCALC; an event that caught less of its rocking curve than --min-captured-fraction (or --min-partiality) is not written, exactly as in the merge, since a full is what the reading program is told these rows are. Intensities carry the deterministic per-reflection corrections and nothing else - Lorentz and polarization in the LP column, the sensor's angle-dependent efficiency in QE and the flight path's attenuation in FLIGHT, so raw counts are I / LP * QE * FLIGHT - while the partiality is not divided out and the per-image scale is not applied at all, since those programs scale the data themselves. Written in --mode mx and --mode scale, and with --no-merge" << std::endl;
    std::cout << "   --no-export-unmerged            Do not write <prefix>_unmerged.mtz. It is the largest file a run produces, so a run whose observations are not going to another scaling program - a regression battery, say - can skip writing it" << std::endl;
    std::cout << "   --export-unmerged-partials      Write <prefix>_unmerged_partials.mtz, the same observations with each partial as its own row (one batch per image) flagged for the reading program to sum, instead of summed here. Off by default and independent of --export-unmerged; both can be written in one run" << std::endl;
@@ -286,6 +287,7 @@ enum {
    OPT_REJECT_OUTLIERS,
    OPT_REFERENCE_COLUMN,
    OPT_MODEL,
+    OPT_AF_MODEL,
    OPT_EXPORT_UNMERGED,
    OPT_NO_EXPORT_UNMERGED,
    OPT_EXPORT_UNMERGED_PARTIALS,
@@ -348,6 +350,7 @@ static option long_options[] = {
    {"reference-mtz", required_argument, nullptr, 'z'},
    {"reference-column", required_argument, nullptr, OPT_REFERENCE_COLUMN},
    {"model", required_argument, nullptr, OPT_MODEL},
+    {"af-model", required_argument, nullptr, OPT_AF_MODEL},
    {"export-unmerged", no_argument, nullptr, OPT_EXPORT_UNMERGED},
    {"no-export-unmerged", no_argument, nullptr, OPT_NO_EXPORT_UNMERGED},
    {"export-unmerged-partials", no_argument, nullptr, OPT_EXPORT_UNMERGED_PARTIALS},
@@ -936,6 +939,7 @@ static int RunRugnux(int argc, char **argv) {
    std::string ref_mtz;
    std::string ref_column;
    std::string model_pdb;           // --model: PDB to validate merged intensities against (R-free + maps)
+    std::string af_model;            // --af-model: predicted model placed by molecular replacement (experimental)
    bool export_unmerged = true;            // on by default; --no-export-unmerged turns it off
    bool export_unmerged_partials = false;  // --export-unmerged-partials: the same, one row per image
    std::string dump_observations;   // diagnostic: dump unmerged -P rot3d fulls to this path
@@ -1160,6 +1164,9 @@ static int RunRugnux(int argc, char **argv) {
        case OPT_MODEL:
            model_pdb = optarg;
            break;
+        case OPT_AF_MODEL:
+            af_model = optarg;
+            break;
        case OPT_EXPORT_UNMERGED:
            export_unmerged = true;
            break;
@@ -1544,6 +1551,11 @@ static int RunRugnux(int argc, char **argv) {
        exit(EXIT_FAILURE);
    }

+    if (!af_model.empty() && mode != RugnuxMode::MX) {
+        logger.Error("--af-model is for --mode mx only");
+        exit(EXIT_FAILURE);
+    }
+
    input_file = argv[optind];
    const std::string broker_url = (optind == argc - 2) ? argv[optind + 1] : "";

@@ -3082,6 +3094,7 @@ static int RunRugnux(int argc, char **argv) {
    config.export_unmerged = export_unmerged;
    config.export_unmerged_partials = export_unmerged_partials;
    config.model_path = model_pdb;
+    config.af_model_path = af_model;
    // When merging, the merged reflections (.mtz/.cif) are the wanted output; skip the large
    // _process.h5 unless explicitly requested. Without merging, the _process.h5 is the only output.
    config.write_process_h5 = run_scaling ? write_process_h5_flag : true;
@@ -143,6 +143,7 @@ ADD_EXECUTABLE(jfjoch_test
        TCPImagePusherTest.cpp
        PreflightTest.cpp
        ModelValidationTest.cpp
+        MolecularReplacementTest.cpp
        ModelMaskGPUTest.cpp
        SearchSpaceGroupTest.cpp
        SearchSpaceGroupTwinTest.cpp
@@ -0,0 +1,166 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+
+#include <cmath>
+#include <random>
+
+#include <gemmi/dencalc.hpp>
+#include <gemmi/it92.hpp>
+
+#include "../image_analysis/structure_refinement/ModelFFT.h"
+#include "../image_analysis/structure_refinement/MolecularReplacement.h"
+
+namespace {
+    // A synthetic body: carbon atoms scattered in a 12 A sphere, centroid at the origin.
+    SearchModel synthetic_model() {
+        SearchModel m;
+        std::mt19937 rng(7);
+        std::uniform_real_distribution<double> u(-12.0, 12.0);
+        gemmi::Model model(1);
+        gemmi::Chain chain("A");
+        gemmi::Residue res;
+        res.name = "UNK";
+        res.seqid = gemmi::SeqId(1, ' ');
+        std::vector<gemmi::Vec3> pos;
+        while (pos.size() < 300) {
+            const gemmi::Vec3 p(u(rng), u(rng), u(rng));
+            if (p.length() < 12.0)
+                pos.push_back(p);
+        }
+        gemmi::Vec3 c;
+        for (const auto &p : pos)
+            c += p;
+        c *= 1.0 / pos.size();
+        for (const auto &p : pos) {
+            gemmi::Atom a;
+            a.name = "C";
+            a.element = gemmi::El::C;
+            a.pos = gemmi::Position(p - c);
+            a.b_iso = 15.0f;
+            a.occ = 1.0f;
+            res.atoms.push_back(a);
+            m.radius_A = std::max(m.radius_A, (p - c).length());
+        }
+        chain.residues.push_back(res);
+        model.chains.push_back(chain);
+        m.st.models.push_back(model);
+        m.residues_read = m.residues_kept = 1;
+        m.mass_da = 3600;
+        return m;
+    }
+
+    gemmi::Mat33 axis_angle(gemmi::Vec3 axis, double deg) {
+        axis = axis.normalized();
+        const double k = deg * M_PI / 180.0, s = std::sin(k), c1 = 1.0 - std::cos(k);
+        const double x = axis.x, y = axis.y, z = axis.z;
+        return {1 + c1 * (x * x - 1), -s * z + c1 * x * y, s * y + c1 * x * z,
+                s * z + c1 * x * y, 1 + c1 * (y * y - 1), -s * x + c1 * y * z,
+                -s * y + c1 * x * z, s * x + c1 * y * z, 1 + c1 * (z * z - 1)};
+    }
+
+    // |F|^2 of the placed crystal to 4 A, as merged reflections.
+    std::vector<MergedReflection> observed(const gemmi::Structure &crystal) {
+        gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
+        dc.d_min = 3.5;
+        dc.rate = 1.5;
+        dc.set_grid_cell_and_spacegroup(crystal);
+        dc.put_model_density_on_grid(crystal.models[0]);
+        const auto f = MapToFPhi(dc.grid).prepare_asu_data(3.5, 0, false, false, false);
+        std::vector<MergedReflection> out;
+        for (const auto &hv : f.v) {
+            MergedReflection r;
+            r.h = hv.hkl[0];
+            r.k = hv.hkl[1];
+            r.l = hv.hkl[2];
+            r.I = std::norm(hv.value);
+            r.sigma = 1;
+            r.F = std::abs(hv.value);
+            r.d = 1.0 / std::sqrt(crystal.cell.calculate_1_d2(hv.hkl));
+            out.push_back(r);
+        }
+        return out;
+    }
+}
+
+TEST_CASE("MolecularReplacement_UniformRotationsCoverSO3", "[MolecularReplacement]") {
+    const auto rotations = UniformRotations(6.0);
+    const double step = 6.0 * M_PI / 180.0;
+    // 8 pi^2 is the volume of SO(3) in the measure a 6 deg cube of Euler-like angles samples.
+    const double expected = 8 * M_PI * M_PI / (step * step * step);
+    CHECK(rotations.size() > 0.8 * expected);
+    CHECK(rotations.size() < 1.2 * expected);
+    // Every rotation has a sampled one nearby.
+    std::mt19937 rng(11);
+    std::normal_distribution<double> g;
+    const std::vector<gemmi::Mat33> identity{gemmi::Mat33()};
+    double worst = 0;
+    for (int i = 0; i < 100; i++) {
+        const gemmi::Mat33 r = axis_angle({g(rng), g(rng), g(rng)}, std::uniform_real_distribution<double>(0, 180)(rng));
+        double best = 180;
+        for (const auto &s : rotations)
+            best = std::min(best, RotationDistanceDeg(r, s, identity));
+        worst = std::max(worst, best);
+    }
+    CHECK(worst < 6.0);
+}
+
+TEST_CASE("MolecularReplacement_RotationAsu", "[MolecularReplacement]") {
+    const gemmi::UnitCell cell(40, 40, 60, 90, 90, 90);
+    const auto pg = CartesianPointGroup(gemmi::get_spacegroup_by_name("P 41 2 2"), cell);
+    REQUIRE(pg.size() == 8);
+    const auto all = UniformRotations(10.0);
+    const auto asu = RotationAsu(all, pg);
+    // One class in eight, give or take the boundaries.
+    CHECK(asu.size() > all.size() / 8 * 0.9);
+    CHECK(asu.size() < all.size() / 8 * 1.3);
+    // Every sampled rotation has a class member in the ASU close by.
+    for (size_t i = 0; i < all.size(); i += 37) {
+        double best = 180;
+        for (const auto &s : asu)
+            best = std::min(best, RotationDistanceDeg(s, all[i], pg));
+        CHECK(best < 10.0);
+    }
+}
+
+TEST_CASE("MolecularReplacement_SameRotationGroups", "[MolecularReplacement]") {
+    const auto groups = SameRotationGroups(gemmi::get_spacegroup_by_name("P 41"));
+    REQUIRE(groups.size() == 4);
+    CHECK(groups[0].xhm() == "P 41");
+    const auto i = SameRotationGroups(gemmi::get_spacegroup_by_name("I 41 2 2"));
+    CHECK(i.size() == 2);   // I 4 2 2, I 41 2 2
+}
+
+TEST_CASE("MolecularReplacement_SyntheticPlacement", "[MolecularReplacement]") {
+    const SearchModel model = synthetic_model();
+    const gemmi::UnitCell cell(42, 50, 58, 90, 90, 90);
+    const gemmi::SpaceGroup &sg = gemmi::get_spacegroup_by_name("P 21 21 21");
+    const gemmi::Mat33 truth = axis_angle({1, 2, 3}, 50);
+    const gemmi::Fractional t(0.17, 0.31, 0.09);
+    const gemmi::Structure crystal = PlaceSearchModel(model, truth, t, cell, sg);
+    const auto merged = observed(crystal);
+    const SearchReflections refl = MakeSearchReflections(merged, cell, sg, 4.0, 12.0);
+    REQUIRE(refl.hkl.size() > 300);
+    const MolecularTransform ft(model, 4.0, 4);
+
+    // The rotation function's best orientation is the true one, to the sampling.
+    const auto pg = CartesianPointGroup(sg, cell);
+    const auto rotations = RotationAsu(UniformRotations(6.0), pg);
+    const auto rf = RotationFunction(ft, refl, cell, sg, rotations, 4);
+    size_t best = 0;
+    for (size_t i = 1; i < rf.size(); i++)
+        if (rf[i] > rf[best])
+            best = i;
+    CHECK(RotationDistanceDeg(rotations[best], truth, pg) < 6.0);
+
+    // At the true orientation the translation function peaks at the true position, up to the
+    // origins P 21 21 21 allows (a half in each coordinate), and stands far out of its own noise.
+    const TranslationPeak peak = TranslationFunction(ft, refl, cell, sg, truth);
+    CHECK(peak.z > 10.0);
+    for (double d : {peak.t.x - t.x, peak.t.y - t.y, peak.t.z - t.z}) {
+        const double half = 2 * d - std::round(2 * d);
+        CHECK(std::fabs(half) < 0.06);
+    }
+}
+