// 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 <functional>
#include <numeric>
#include <optional>
#include <map>
#include <mutex>
#include <tuple>
#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 "ModelValidation.h"
#include "RigidBodyRefine.h"
#ifdef JFJOCH_USE_CUDA
#include "MolecularReplacementGPU.h"
#include "RigidBodyGPU.h"
#include "RigidBodyGPUEngine.h"
#include "ModelStructureFactorsGPU.h"
#include "../../common/CUDAWrapper.h"
#include "../../common/JFJochException.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;
// Hill-climbing steps per step size of a rotation-function peak's refinement.
constexpr int PEAK_REFINE_ITERATIONS = 4;
// 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.
constexpr double VM_MIN = 1.5;
// Copies are searched for one after another up to the largest copy number at least this probable
// (MatthewsProbabilities); a copy is kept where its translation function stands out beside the copies
// already placed, it packs, and the R-free of all of them drops by at least COPY_R_FREE_GAIN.
constexpr double MATTHEWS_P_MIN = 0.01;
constexpr double COPY_R_FREE_GAIN = 0.01;
// The translation-function peaks of a further copy its assembly correlation is computed for.
constexpr size_t COPY_TRIALS = 10;
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;
// 3 sd of the null, one-sided, allowing for its mean and sd being estimated from NULL_REPLICATES draws:
// a new draw lies t(n-1) sqrt(1 + 1/n) sample sds above their mean at that probability, 4.28 x 1.054 for
// n = 9 - the gate model validation's fit is read at (ModelValidation.cpp, MODEL_FIT_SIGMA).
constexpr double NULL_SIGMA_MIN = 4.51;

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)};
}

// A rotation by `deg` about Cartesian axis 0, 1 or 2.
gemmi::Mat33 axis_rotation(int axis, double deg) {
    const double c = std::cos(deg * PI / 180.0), s = std::sin(deg * PI / 180.0);
    if (axis == 0)
        return {1, 0, 0, 0, c, -s, 0, s, c};
    if (axis == 1)
        return {c, 0, s, 0, 1, 0, -s, 0, c};
    return {c, -s, 0, s, c, 0, 0, 0, 1};
}

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;
    double k_sol = NAN, b_sol = NAN;
};

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 on the CPU, as model validation makes them.
void model_factors(const gemmi::Structure &st, double d_min, gemmi::AsuData<std::complex<float>> &fmodel,
                   gemmi::AsuData<std::complex<float>> &fmask, 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);
    fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0);
    PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
    fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
}

// R-work and R-free of the scaled Fcalc over the reflections to d_min.
PlacementScore r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, const gemmi::Scaling<float> &scaling,
                         const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged, double d_min) {
    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.k_sol = scaling.k_sol;
    out.b_sol = scaling.b_sol;
    out.r_work = den_w > 0 ? num_w / den_w : 1;
    out.r_free = den_f > 0 ? num_f / den_f : 1;
    return out;
}

// Model validation's structure-factor engines on the GPU, one per group, copy count and resolution, made
// where the first placement of that kind is scored - or none, where there is no GPU, the cell is too small
// for its gridding or one would take over a quarter of the card. Decided on the card's total memory, so
// the same input always takes the same path.
#ifdef JFJOCH_USE_CUDA
class ScoreEngines {
public:
    ModelStructureFactorsGPU *get(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, double d_min) {
        if (get_gpu_count() == 0)
            return nullptr;
        size_t atoms = 0;
        for (const gemmi::Chain &ch : st.models[0].chains)
            for (const gemmi::Residue &r : ch.residues)
                atoms += r.atoms.size();
        std::lock_guard lock(m_);
        auto &e = engines_[{sg.xhm(), atoms, d_min}];
        if (!e.tried) {
            e.tried = true;
            try {
                auto sf = std::make_unique<ModelStructureFactorsGPU>(RigidBodyGPUEngine::CurrentDevice(), st.models[0],
                                                                     st.cell, sg, d_min);
                if (sf->Supported() && sf->DeviceBytes() <= sf->DeviceTotalMemory() / 4) {
                    sf->Reserve();
                    e.sf = std::move(sf);
                }
            } catch (const JFJochException &) {
                cuda_clear_error();
            }
        }
        return e.sf.get();
    }

private:
    struct Entry {
        bool tried = false;
        std::unique_ptr<ModelStructureFactorsGPU> sf;
    };
    std::map<std::tuple<std::string, size_t, double>, Entry> engines_;
    std::mutex m_;
};
#else
class ScoreEngines {};
#endif

// The model's Fcalc and bulk-solvent mask to d_min - on the GPU where `engines` has one - scaled to the
// working amplitudes as model validation scales them (ModelValidation.cpp, fit_model), and the R-factors
// over the reflections to d_min. `solvent`: k_sol and b_sol held at its, not fitted.
PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg,
                               const std::vector<MergedReflection> &merged, double d_min, size_t nthreads,
                               ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr) {
    gemmi::AsuData<std::complex<float>> fmask, fmodel;
#ifdef JFJOCH_USE_CUDA
    if (ModelStructureFactorsGPU *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr)
        gpu->Compute(st.models[0], fmodel, fmask);
    else
#endif
        model_factors(st, d_min, fmodel, fmask, nthreads);

    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);
    if (solvent != nullptr) {
        // The bulk solvent of an earlier fit of the same crystal, held, as the rigid body holds it.
        scaling.k_sol = solvent->k_sol;
        scaling.b_sol = solvent->b_sol;
        scaling.fix_k_sol = true;
        scaling.fix_b_sol = true;
        scaling.fit_isotropic_b_approximately();
        scaling.fit_parameters();
    } else {
        FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads);
    }
    scaling.scale_data(fmodel, &fmask);
    return r_factors(fmodel, scaling, sg, merged, d_min);
}

// 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.
// `moved`: where given, the rotation the rigid body turned the model by about its centroid. `fixed`: the
// other copies, held where they are.
void refine_rigid_body(gemmi::Structure &st, const gemmi::SpaceGroup &sg, const std::vector<MergedReflection> &merged,
                       double d_min, RigidBodyGPUPool *pool, size_t nthreads, gemmi::Mat33 *moved = nullptr,
                       const gemmi::Model *fixed = nullptr) {
    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, fixed);
    if (!rb.converged)
        SetModelPositions(st.models[0], start);
    if (moved != nullptr)
        *moved = rb.converged ? rb.rotation : gemmi::Mat33();
}

PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg,
                                const std::vector<MergedReflection> &merged, double d_min, RigidBodyGPUPool *pool,
                                size_t nthreads, ScoreEngines *engines, gemmi::Mat33 *moved = nullptr) {
    refine_rigid_body(st, sg, merged, d_min, pool, nthreads, moved);
    return score_placement(st, sg, merged, d_min, nthreads, engines);
}

// The fractional position of the model's centroid.
gemmi::Fractional centroid(const gemmi::Structure &st) {
    const std::vector<gemmi::Position> pos = ModelPositions(st.models[0]);
    gemmi::Vec3 c;
    for (const gemmi::Position &p : pos)
        c += p;
    return st.cell.fractionalize(gemmi::Position(c * (1.0 / static_cast<double>(pos.size()))));
}

// 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.
// `translate`: where given, each replicate is placed where the translation function puts the model in its
// orientation (the centroids, fractional, of the model in each of the orientations given), as the search
// placed the real one, instead of at the real placement's centroid; the search chose that position among
// all others for the real orientation, and a null that is not given the same choice reads high.
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, ScoreEngines *scorers,
                  const std::function<std::vector<gemmi::Fractional>(const std::vector<gemmi::Mat33> &)> &translate) {
    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, scorers).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<gemmi::Fractional> where;
    if (translate) {
        std::vector<gemmi::Mat33> turned;
        for (const gemmi::Mat33 &r : draws)
            turned.push_back(r.multiply(orientation));
        where = translate(turned);
    }
    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());
                const gemmi::Vec3 to = translate ? gemmi::Vec3(placed.cell.orthogonalize(where[i])) : c;
                for (gemmi::Position &p : pos)
                    p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + to);
                SetModelPositions(rep.models[0], pos);
                r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads, scorers).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;
}

// Whether b is a on a shifted origin, the shift on the twelfths every crystallographic one lies on.
bool same_up_to_origin(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
    if (a.number != b.number)
        return false;
    const gemmi::GroupOps bops = b.operations();
    for (int i = 0; i < 12; i++)
        for (int j = 0; j < 12; j++)
            for (int k = 0; k < 12; k++) {
                gemmi::Op shift = gemmi::Op::identity();
                shift.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12, k * gemmi::Op::DEN / 12};
                gemmi::GroupOps aops = a.operations();
                aops.change_basis_forward(shift);
                if (aops.is_same_as(bops))
                    return true;
            }
    return false;
}

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<gemmi::Vec3> CopyVectors(const SearchReflections &refl, const gemmi::UnitCell &cell,
                                     const gemmi::SpaceGroup &sg) {
    const gemmi::GroupOps gops = sg.operations();
    const size_t nops = gops.sym_ops.size();
    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
    std::vector<gemmi::Vec3> s(refl.hkl.size() * nops);
    for (size_t i = 0; i < refl.hkl.size(); 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])));
    return s;
}

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 size_t n = refl.hkl.size();
    const std::vector<gemmi::Vec3> s = CopyVectors(refl, cell, sg);
    const size_t nops = s.size() / n;
    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_sum = 0, x_ss = 0, xy = 0;
        for (size_t i = 0; i < n; i++) {
            const double x = ic[i] / (bin_sum[refl.bin[i]] / bin_n[refl.bin[i]]);
            x_sum += x;
            x_ss += x * x;
            xy += x * (refl.e2[i] - w_mean);
        }
        const double sxx = x_ss - x_sum * x_sum / static_cast<double>(n);
        out[r] = sxx > 0 && w_ss > 0 ? xy / std::sqrt(sxx * w_ss) : 0.0;
    });
    return out;
}

// --- the translation function ---
//
// 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), g_j = G_j(h) exp(2 pi i h.t_j), 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.

TranslationSetup MakeTranslationSetup(const SearchReflections &refl, const gemmi::UnitCell &cell,
                                      const gemmi::SpaceGroup &sg, bool known_part) {
    const gemmi::GroupOps gops = sg.operations();
    const int nops = static_cast<int>(gops.sym_ops.size());
    const gemmi::Mat33 frac_t = cell.frac.mat.transpose();
    TranslationSetup out;
    out.ops = nops;
    out.bins = refl.bins;
    for (int j = 0; j < nops; j++)
        for (int k = 0; k < nops; k++)
            if (j != k) {
                out.pair_j.push_back(j);
                out.pair_k.push_back(k);
            }
    std::vector<gemmi::Miller> rotated;   // terms x ops
    for (size_t i = 0; i < refl.hkl.size(); i++) {
        const gemmi::Miller &h = refl.hkl[i];
        if (gops.is_systematically_absent(h))
            continue;
        for (const gemmi::Op &op : gops.sym_ops) {
            rotated.push_back(rotated_index(op, h));
            out.s.push_back(frac_t.multiply(miller_vec(rotated.back())));
            const double ht = (h[0] * op.tran[0] + h[1] * op.tran[1] + h[2] * op.tran[2])
                              / static_cast<double>(gemmi::Op::DEN);
            out.phase.push_back(std::polar(1.0, 2.0 * PI * ht));
        }
        out.weight.push_back(refl.multiplicity[i] * (refl.e2[i] - 1.0));
        out.e2.push_back(refl.e2[i]);
        out.multiplicity.push_back(refl.multiplicity[i]);
        out.bin.push_back(refl.bin[i]);
    }
    int max_m[3] = {0, 0, 0};
    for (size_t t = 0; t < out.Terms(); t++)
        for (size_t p = 0; p < out.pair_j.size(); p++)
            for (int a = 0; a < 3; a++)
                max_m[a] = std::max(max_m[a], std::abs(rotated[t * nops + out.pair_j[p]][a]
                                                       - rotated[t * nops + out.pair_k[p]][a]));
    if (known_part)
        for (const gemmi::Miller &k : rotated)
            for (int a = 0; a < 3; a++)
                max_m[a] = std::max(max_m[a], std::abs(k[a]));
    out.nu = good_even_size(2 * max_m[0] + 2);
    out.nv = good_even_size(2 * max_m[1] + 2);
    out.nw = good_even_size(2 * max_m[2] + 2);
    for (size_t t = 0; t < out.Terms(); t++)
        for (size_t p = 0; p < out.pair_j.size(); p++) {
            const gemmi::Miller &kj = rotated[t * nops + out.pair_j[p]], &kk = rotated[t * nops + out.pair_k[p]];
            const int l = kj[2] - kk[2];
            out.target_u.push_back(l < 0 ? -1 : ((kj[0] - kk[0]) % out.nu + out.nu) % out.nu);
            out.target_v.push_back(((kj[1] - kk[1]) % out.nv + out.nv) % out.nv);
            out.target_l.push_back(l);
        }
    if (known_part)
        for (size_t t = 0; t < out.Terms(); t++)
            for (int j = 0; j < nops; j++) {
                const gemmi::Miller &k = rotated[t * nops + j];
                const auto add = [&](int sign) {
                    out.cross_term.push_back(static_cast<int>(t));
                    out.cross_op.push_back(j);
                    out.cross_u.push_back(((sign * k[0]) % out.nu + out.nu) % out.nu);
                    out.cross_v.push_back(((sign * k[1]) % out.nv + out.nv) % out.nv);
                    out.cross_l.push_back(sign * k[2]);
                    out.cross_conj.push_back(sign < 0);
                };
                if (k[2] >= 0)
                    add(1);
                if (k[2] <= 0)
                    add(-1);
            }
    return out;
}

std::vector<std::complex<double>> KnownStructureFactors(const MolecularTransform &ft, const TranslationSetup &setup,
                                                        const gemmi::UnitCell &cell,
                                                        const std::vector<gemmi::Mat33> &rotations,
                                                        const std::vector<gemmi::Fractional> &positions) {
    const size_t terms = setup.Terms(), nops = setup.ops;
    std::vector<std::complex<double>> out(terms, 0.0);
    for (size_t c = 0; c < rotations.size(); c++) {
        const gemmi::Mat33 rt = rotations[c].transpose();
        const gemmi::Position x = cell.orthogonalize(positions[c]);
        for (size_t t = 0; t < terms; t++)
            for (size_t j = 0; j < nops; j++) {
                const gemmi::Vec3 &sj = setup.s[t * nops + j];
                // (h R_j).t = s_j.x for the Cartesian s_j and x
                out[t] += ft.At(rt.multiply(sj)) * setup.phase[t * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x));
            }
    }
    return out;
}

double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell,
                           const std::vector<std::complex<double>> &known, const gemmi::Mat33 &rotation,
                           const gemmi::Fractional &t) {
    const size_t terms = setup.Terms(), nops = setup.ops;
    const gemmi::Mat33 rt = rotation.transpose();
    const gemmi::Position x = cell.orthogonalize(t);
    std::vector<double> ic(terms), bin_sum(setup.bins, 0.0);
    std::vector<int> bin_n(setup.bins, 0);
    for (size_t i = 0; i < terms; i++) {
        std::complex<double> f = known[i];
        for (size_t j = 0; j < nops; j++) {
            const gemmi::Vec3 &sj = setup.s[i * nops + j];
            f += ft.At(rt.multiply(sj)) * setup.phase[i * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x));
        }
        ic[i] = std::norm(f);
        bin_sum[setup.bin[i]] += ic[i];
        bin_n[setup.bin[i]]++;
    }
    double w = 0, sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0;
    for (size_t i = 0; i < terms; i++) {
        const double m = setup.multiplicity[i], a = ic[i] / (bin_sum[setup.bin[i]] / bin_n[setup.bin[i]]),
                     b = setup.e2[i];
        w += m;
        sx += m * a;
        sy += m * b;
        sxx += m * a * a;
        syy += m * b * b;
        sxy += m * a * b;
    }
    const double vx = sxx - sx * sx / w, vy = syy - sy * sy / w;
    return vx > 0 && vy > 0 ? (sxy - sx * sy / w) / std::sqrt(vx * vy) : 0.0;
}

TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup,
                                    const gemmi::Mat33 &rotation, const std::vector<std::complex<double>> *known) {
    TranslationPeak out;
    const size_t terms = setup.Terms(), nops = setup.ops, npairs = setup.pair_j.size();
    if (terms == 0 || (nops < 2 && known == nullptr))
        return out;   // one copy in P1: every position is the same
    const gemmi::Mat33 rt = rotation.transpose();
    std::vector<std::complex<double>> g(terms * nops);
    std::vector<double> bin_sum(setup.bins, 0.0);
    std::vector<int> bin_n(setup.bins, 0);
    for (size_t t = 0; t < terms; t++) {
        double sum = 0;
        for (size_t j = 0; j < nops; j++) {
            g[t * nops + j] = ft.At(rt.multiply(setup.s[t * nops + j])) * setup.phase[t * nops + j];
            sum += std::norm(g[t * nops + j]);
        }
        bin_sum[setup.bin[t]] += sum;
        bin_n[setup.bin[t]]++;
    }
    const int nu = setup.nu, nv = setup.nv, nwh = setup.nw / 2 + 1;
    std::vector<std::complex<double>> c(static_cast<size_t>(nu) * nv * nwh, 0.0);
    for (size_t t = 0; t < terms; t++) {
        const double w = setup.weight[t] / (bin_sum[setup.bin[t]] / bin_n[setup.bin[t]]);
        for (size_t p = 0; p < npairs; p++) {
            const size_t e = t * npairs + p;
            if (setup.target_u[e] < 0)
                continue;
            c[setup.target_u[e] + static_cast<size_t>(nu) * (setup.target_v[e] + static_cast<size_t>(nv) * setup.target_l[e])]
                    += w * g[t * nops + setup.pair_j[p]] * std::conj(g[t * nops + setup.pair_k[p]]);
        }
    }
    if (known != nullptr)
        for (size_t e = 0; e < setup.cross_term.size(); e++) {
            const size_t t = setup.cross_term[e];
            const double w = setup.weight[t] / (bin_sum[setup.bin[t]] / bin_n[setup.bin[t]]);
            const std::complex<double> v = std::conj((*known)[t]) * g[t * nops + setup.cross_op[e]];
            c[setup.cross_u[e] + static_cast<size_t>(nu) * (setup.cross_v[e] + static_cast<size_t>(nv) * setup.cross_l[e])]
                    += w * (setup.cross_conj[e] ? std::conj(v) : 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 = ft.Grid().unit_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;
}

// --- translational NCS ---

TranslationalNcs DetectTranslationalNcs(const std::vector<MergedReflection> &merged, const gemmi::UnitCell &cell,
                                        const gemmi::SpaceGroup &sg) {
    constexpr double D_MIN = 5.0, D_MAX = 10.0, MIN_DISTANCE_A = 15.0, MIN_HEIGHT = 0.2;
    TranslationalNcs out;
    const SearchReflections refl = MakeSearchReflections(merged, cell, sg, D_MIN, D_MAX);
    const gemmi::GroupOps gops = sg.operations();
    int max_m[3] = {0, 0, 0};
    std::vector<std::pair<gemmi::Miller, double>> terms;   // the whole sphere, l >= 0
    for (size_t i = 0; i < refl.hkl.size(); i++)
        for (const gemmi::Op &op : gops.sym_ops)
            for (int sign : {1, -1}) {
                gemmi::Miller k = rotated_index(op, refl.hkl[i]);
                for (int &x : k)
                    x *= sign;
                if (k[2] < 0)
                    continue;
                terms.push_back({k, refl.e2[i]});
                for (int a = 0; a < 3; a++)
                    max_m[a] = std::max(max_m[a], std::abs(k[a]));
            }
    if (terms.empty())
        return out;
    // Sampled at a third of D_MIN at least.
    gemmi::FPhiGrid<float> hkl;
    hkl.unit_cell = cell;
    hkl.spacegroup = &gemmi::get_spacegroup_p1();
    hkl.axis_order = gemmi::AxisOrder::XYZ;
    hkl.half_l = true;
    const int nu = good_even_size(std::max(2 * max_m[0] + 2, static_cast<int>(std::ceil(3 * cell.a / D_MIN))));
    const int nv = good_even_size(std::max(2 * max_m[1] + 2, static_cast<int>(std::ceil(3 * cell.b / D_MIN))));
    const int nw = good_even_size(std::max(2 * max_m[2] + 2, static_cast<int>(std::ceil(3 * cell.c / D_MIN))));
    hkl.set_size_without_checking(nu, nv, nw / 2 + 1);
    for (const auto &[k, e2] : terms) {
        const size_t idx = ((k[0] % nu + nu) % nu) + static_cast<size_t>(nu) * (((k[1] % nv + nv) % nv)
                                                                               + static_cast<size_t>(nv) * k[2]);
        hkl.data[idx] = std::complex<float>(static_cast<float>(e2), 0.0f);
    }
    const gemmi::Grid<float> map = MapFromFPhi(hkl);
    const double origin = map.data[0];
    for (int w = 0; w < map.nw; w++)
        for (int v = 0; v < map.nv; v++)
            for (int u = 0; u < map.nu; u++) {
                const double value = map.data[u + static_cast<size_t>(map.nu) * (v + static_cast<size_t>(map.nv) * w)];
                if (!(value > out.height * origin))
                    continue;
                const gemmi::Fractional f(static_cast<double>(u) / map.nu, static_cast<double>(v) / map.nv,
                                          static_cast<double>(w) / map.nw);
                bool near_lattice = false;
                for (const gemmi::Op::Tran &c : gops.cen_ops) {
                    gemmi::Fractional d(f.x - static_cast<double>(c[0]) / gemmi::Op::DEN,
                                        f.y - static_cast<double>(c[1]) / gemmi::Op::DEN,
                                        f.z - static_cast<double>(c[2]) / gemmi::Op::DEN);
                    d = gemmi::Fractional(d.x - std::round(d.x), d.y - std::round(d.y), d.z - std::round(d.z));
                    if (cell.orthogonalize_difference(d).length() < MIN_DISTANCE_A)
                        near_lattice = true;
                }
                if (near_lattice)
                    continue;
                out.height = value / origin;
                out.vector = f;
            }
    out.present = out.height >= MIN_HEIGHT;
    return out;
}

bool TncsRelatesGroups(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b, const gemmi::Fractional &v) {
    constexpr double TOLERANCE = 0.05;
    const auto mod1 = [](double x) { return x - std::round(x); };
    const gemmi::GroupOps bops = b.operations();
    for (const gemmi::Op &op : a.operations().sym_ops) {
        const gemmi::Op *other = bops.find_by_rotation(op.rot);
        if (other == nullptr)
            return false;
        // The screw parts, (1/n) sum over k of R^k t, and the same of v.
        double wa[3] = {0, 0, 0}, wb[3] = {0, 0, 0}, wv[3] = {0, 0, 0};
        gemmi::Op power = gemmi::Op::identity();
        int n = 0;
        do {
            for (int i = 0; i < 3; i++) {
                double ta = 0, tb = 0, tv = 0;
                for (int j = 0; j < 3; j++) {
                    const double r = static_cast<double>(power.rot[i][j]) / gemmi::Op::DEN;
                    ta += r * op.tran[j] / gemmi::Op::DEN;
                    tb += r * other->tran[j] / gemmi::Op::DEN;
                    tv += r * (j == 0 ? v.x : j == 1 ? v.y : v.z);
                }
                wa[i] += ta;
                wb[i] += tb;
                wv[i] += tv;
            }
            power = power.combine(gemmi::Op{op.rot, {0, 0, 0}, 'x'});
            n++;
        } while (!(power.rot == gemmi::Op::identity().rot) && n < 6);
        bool bridged = false;
        for (int m = 0; m < n && !bridged; m++) {
            bridged = true;
            for (int i = 0; i < 3; i++)
                if (std::fabs(mod1((wa[i] - wb[i] - m * wv[i]) / n)) > TOLERANCE)
                    bridged = false;
        }
        if (!bridged)
            return false;
    }
    return true;
}

// --- 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, size_t copies) {
    const gemmi::UnitCell &cell = placed.cell;
    const gemmi::SpaceGroup *sg = placed.find_spacegroup();
    const size_t chains = placed.models[0].chains.size();
    copies = std::max<size_t>(1, std::min(copies, chains));
    // Each copy's C-alpha atoms, fractional, its centroid and its radius.
    struct Copy {
        std::vector<gemmi::Fractional> ca;
        gemmi::Fractional c0;
        double radius = 0;
    };
    std::vector<Copy> body(copies);
    for (size_t i = 0; i < chains; i++)
        for (const gemmi::Residue &r : placed.models[0].chains[i].residues)
            for (const gemmi::Atom &a : r.atoms)
                if (a.name == "CA")
                    body[i * copies / chains].ca.push_back(cell.fractionalize(a.pos));
    size_t total = 0;
    for (Copy &b : body) {
        total += b.ca.size();
        for (const gemmi::Fractional &f : b.ca)
            b.c0 = b.c0 + f;
        if (!b.ca.empty())
            b.c0 = gemmi::Fractional(b.c0 * (1.0 / static_cast<double>(b.ca.size())));
        for (const gemmi::Fractional &f : b.ca)
            b.radius = std::max(b.radius, cell.orthogonalize_difference(f - b.c0).length());
    }
    if (total == 0 || sg == nullptr)
        return 0.0;
    constexpr double CLASH_A = 3.0;
    // Every C-alpha of every image that can reach a copy - each copy under each operator and the lattice
    // shifts about it - binned in CLASH_A cubes, sorted by cube; a C-alpha clashes where one of another
    // image (not its own copy unmoved) sits within CLASH_A of it in its own cube or a neighbouring one.
    struct Image {
        long key;
        gemmi::Position x;
        int copy;
        bool unmoved;
    };
    const auto cube = [](const gemmi::Position &x, int dx, int dy, int dz) {
        return (static_cast<long>(std::floor(x.x / CLASH_A)) + dx + 4096) * 67108864L
               + (static_cast<long>(std::floor(x.y / CLASH_A)) + dy + 4096) * 8192L
               + (static_cast<long>(std::floor(x.z / CLASH_A)) + dz + 4096);
    };
    std::vector<Image> images;
    for (size_t b = 0; b < copies; b++)
        for (const gemmi::Op &op : sg->operations().all_ops_sorted()) {
            std::vector<gemmi::Fractional> mate;
            for (const gemmi::Fractional &f : body[b].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())));
            // The lattice shifts that bring this image within reach of some copy.
            std::vector<gemmi::Fractional> shifts;
            for (size_t a = 0; a < copies; a++) {
                const gemmi::Fractional nearest(std::round(body[a].c0.x - cm.x), std::round(body[a].c0.y - cm.y),
                                                std::round(body[a].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 (cell.orthogonalize_difference(cm + shift - body[a].c0).length()
                                    <= body[a].radius + body[b].radius + CLASH_A
                                && std::none_of(shifts.begin(), shifts.end(), [&](const gemmi::Fractional &o) {
                                       return o.x == shift.x && o.y == shift.y && o.z == shift.z;
                                   }))
                                shifts.push_back(shift);
                        }
            }
            for (const gemmi::Fractional &shift : shifts) {
                const bool unmoved = op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0;
                for (const gemmi::Fractional &m : mate) {
                    const gemmi::Position x = cell.orthogonalize(m + shift);
                    images.push_back({cube(x, 0, 0, 0), x, static_cast<int>(b), unmoved});
                }
            }
        }
    std::sort(images.begin(), images.end(), [](const Image &p, const Image &q) { return p.key < q.key; });
    size_t clashing = 0;
    for (size_t a = 0; a < copies; a++)
        for (const gemmi::Fractional &f : body[a].ca) {
            const gemmi::Position x = cell.orthogonalize(f);
            bool clash = false;
            for (int ex = -1; ex <= 1 && !clash; ex++)
                for (int ey = -1; ey <= 1 && !clash; ey++)
                    for (int ez = -1; ez <= 1 && !clash; ez++) {
                        const long k = cube(x, ex, ey, ez);
                        auto it = std::lower_bound(images.begin(), images.end(), k,
                                                   [](const Image &p, long v) { return p.key < v; });
                        for (; it != images.end() && it->key == k; ++it)
                            if (!(it->unmoved && it->copy == static_cast<int>(a)) && x.dist_sq(it->x) < CLASH_A * CLASH_A) {
                                clash = true;
                                break;
                            }
                    }
            clashing += clash;
        }
    return static_cast<double>(clashing) / static_cast<double>(total);
}

gemmi::Structure AssembleCopies(const std::vector<gemmi::Structure> &copies) {
    static const std::string names = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789";
    gemmi::Structure out = copies.at(0);
    out.models[0].chains.clear();
    size_t n = 0;
    for (const gemmi::Structure &c : copies)
        for (gemmi::Chain ch : c.models[0].chains) {
            ch.name = names[n % names.size()] + (n < names.size() ? std::string() : std::to_string(n / names.size()));
            out.models[0].chains.push_back(std::move(ch));
            n++;
        }
    return out;
}

// The density of the Matthews coefficient over protein crystals: a Gumbel (extreme-value) distribution
// - the family Kantardjieff & Rupp fit to this highly skewed peak - with their mode (2.34 A^3/Da) and
// mean (2.69 A^3/Da) for 10471 protein crystal forms, which make its scale (mean - mode) / Euler's
// constant; without their resolution dependence.
// Following Kantardjieff & Rupp (2003) Protein Sci. 12, 1865-1871
std::vector<double> MatthewsProbabilities(double asu_volume, double mass_da) {
    constexpr double MODE = 2.34, MEAN = 2.69, EULER_GAMMA = 0.5772156649;
    constexpr double SCALE = (MEAN - MODE) / EULER_GAMMA;
    std::vector<double> p;
    for (int n = 1; asu_volume / (n * mass_da) >= VM_MIN; n++) {
        const double z = (asu_volume / (n * mass_da) - MODE) / SCALE;
        p.push_back(std::exp(-z - std::exp(-z)));
    }
    const double sum = std::accumulate(p.begin(), p.end(), 0.0);
    for (double &x : p)
        x /= sum;
    return p;
}

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) {
        const gemmi::GroupOps gops = g.operations();
        if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g)
            continue;
        // The table lists some groups a second time on another origin (P 21 21 2 and "P 21212(a)",
        // I 2 2 2 and "I 2 2 2a"): the same crystal, which would tie with itself in the choice.
        if (std::none_of(out.begin(), out.end(), [&](const gemmi::SpaceGroup &k) { return same_up_to_origin(k, g); }))
            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);
#ifdef JFJOCH_USE_CUDA
    // The rotation and translation functions on the GPU where there is one; the rest of the search is
    // the same either way.
    std::unique_ptr<MolecularReplacementGPU> gpu;
    if (get_gpu_count() > 0) {
        try {
            gpu = std::make_unique<MolecularReplacementGPU>(ft);
        } catch (const JFJochException &e) {
            result.failure_reason = std::string("GPU: ") + e.what();
            return result;
        }
    }
    result.gpu_search = gpu != nullptr;
#endif
    result.t_transform = seconds_since(t0);
    result.tncs = DetectTranslationalNcs(*merges[0].merged, gemmi::UnitCell(merges[0].cell), merges[0].space_group);
    logger.Info("Model prior: native Patterson's highest peak off the lattice {:.2f} of the origin at {:.3f} {:.3f} {:.3f}{}",
                result.tncs.height, result.tncs.vector.x, result.tncs.vector.y, result.tncs.vector.z,
                result.tncs.present ? " - tNCS" : "");
    const std::vector<gemmi::Mat33> sphere = UniformRotations(ROTATION_STEP_DEG);
    ScoreEngines scorers;

    std::vector<gemmi::Mat33> orientations;   // of each candidate's placement, before the rigid body
    // 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);
#ifdef JFJOCH_USE_CUDA
        const std::vector<double> rf = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, rotations)
                                           : RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
#else
        const std::vector<double> rf = RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads);
#endif
        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;
        }
        // Each peak refined to the function's local maximum, off the sampling: a large model turned by
        // half a sampling step is a few angstroms out at its surface, which costs the translation
        // function much of its peak. Hill climbing about the three axes, the step halved from half the
        // sampling to an eighth.
        std::vector<gemmi::Mat33> peak_rotations;
        std::vector<double> peak_rf;
        for (size_t p : peaks) {
            peak_rotations.push_back(rotations[p]);
            peak_rf.push_back(rf[p]);
        }
        for (double step = ROTATION_STEP_DEG / 2; step >= ROTATION_STEP_DEG / 8 - 1e-9; step /= 2) {
            std::vector<char> active(peaks.size(), 1);
            for (int iter = 0; iter < PEAK_REFINE_ITERATIONS
                               && std::find(active.begin(), active.end(), 1) != active.end(); iter++) {
                std::vector<gemmi::Mat33> trial;
                for (size_t k = 0; k < peaks.size(); k++)
                    if (active[k])
                        for (int axis = 0; axis < 3; axis++)
                            for (double sign : {1.0, -1.0})
                                trial.push_back(axis_rotation(axis, sign * step).multiply(peak_rotations[k]));
#ifdef JFJOCH_USE_CUDA
                const std::vector<double> v = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, trial)
                                                  : RotationFunction(ft, refl, cell, merge.space_group, trial, nthreads);
#else
                const std::vector<double> v = RotationFunction(ft, refl, cell, merge.space_group, trial, nthreads);
#endif
                size_t next = 0;
                for (size_t k = 0; k < peaks.size(); k++) {
                    if (!active[k])
                        continue;
                    size_t best = next;
                    for (size_t j = next; j < next + 6; j++)
                        if (v[j] > v[best])
                            best = j;
                    if (v[best] > peak_rf[k]) {
                        peak_rf[k] = v[best];
                        peak_rotations[k] = trial[best];
                    } else {
                        active[k] = 0;
                    }
                    next += 6;
                }
            }
        }
        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 ? (peak_rf.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;
        for (const gemmi::SpaceGroup &sg : groups) {
            const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
            std::vector<TranslationPeak> tf(peaks.size());
#ifdef JFJOCH_USE_CUDA
            if (gpu)
                tf = gpu->TranslationFunction(setup, peak_rotations);
            else
#endif
                ParallelFor(static_cast<int>(peaks.size()), nthreads, [&](int i) {
                    tf[i] = TranslationFunction(ft, setup, peak_rotations[i]);
                });
            tf_all.insert(tf_all.end(), tf.begin(), tf.end());
        }
        result.t_translation += seconds_since(t0);

        // Each group searched on a thread of its own - they share nothing but the read-only search and the
        // score engines - its log kept and replayed in order.
        struct GroupResult {
            MRCandidate cand;
            gemmi::Mat33 orientation;
            std::vector<gemmi::Mat33> copy_rotation;   // of every copy placed
            Logger log = Logger::Buffered();
        };
        const auto search_group = [&](size_t g, const std::vector<gemmi::Mat33> &seeds) {
            GroupResult out;
            Logger &glog = out.log;
            const gemmi::SpaceGroup &sg = groups[g];
            MRCandidate cand;
            cand.merge = static_cast<int>(m);
            cand.space_group = sg;
            const gemmi::GroupOps gops = sg.operations();
            const double asu_volume = cell.volume / static_cast<double>(gops.order());
            cand.vm = asu_volume / model.mass_da;
            const std::vector<double> pm = MatthewsProbabilities(asu_volume, model.mass_da);
            for (size_t n = 0; n < pm.size(); n++) {
                if (pm[n] > pm[cand.copies_likely - 1])
                    cand.copies_likely = static_cast<int>(n) + 1;
                if (pm[n] >= MATTHEWS_P_MIN)
                    cand.copies_max = static_cast<int>(n) + 1;
            }
            // The orientations searched: the rotation function's peaks and the seeds, which are tried first.
            std::vector<gemmi::Mat33> rots = peak_rotations;
            std::vector<double> rfs = peak_rf;
            std::vector<TranslationPeak> tfs(tf_all.begin() + g * peaks.size(), tf_all.begin() + (g + 1) * 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 tfs[a].z > tfs[b].z; });
            if (!seeds.empty()) {
                const TranslationSetup setup = MakeTranslationSetup(refl, cell, sg);
                std::vector<TranslationPeak> tf_seed(seeds.size());
#ifdef JFJOCH_USE_CUDA
                const std::vector<double> rf_seed = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, seeds)
                                                        : RotationFunction(ft, refl, cell, merge.space_group, seeds, nthreads);
                if (gpu)
                    tf_seed = gpu->TranslationFunction(setup, seeds);
                else
#else
                const std::vector<double> rf_seed = RotationFunction(ft, refl, cell, merge.space_group, seeds, nthreads);
#endif
                    for (size_t k = 0; k < seeds.size(); k++)
                        tf_seed[k] = TranslationFunction(ft, setup, seeds[k]);
                std::vector<size_t> first;
                for (size_t k = 0; k < seeds.size(); k++) {
                    first.push_back(rots.size());
                    rots.push_back(seeds[k]);
                    rfs.push_back(rf_seed[k]);
                    tfs.push_back(tf_seed[k]);
                }
                by_tf.insert(by_tf.begin(), first.begin(), first.end());
            }
            const TranslationPeak *tf = tfs.data();

            // --- the first copy: the best few through packing, the rigid body and R-free ---
            RigidBodyEngines engines;
            bool have = false;
            gemmi::Mat33 orientation;
            // Where not even one copy fits the cell, only the best translation-function placement is
            // checked for packing.
            const bool fits = cand.vm >= VM_MIN;
            // Each copy placed: its orientation and its centroid, where the rigid body left them.
            std::vector<gemmi::Structure> copies;
            std::vector<gemmi::Mat33> copy_rotation;
            std::vector<gemmi::Fractional> copy_position;
            // Down the list until RESCORED_PER_GROUP placements that pack have been refined: a placement
            // clashing with its own mates over twice the limit is not refined - the rigid body moves a
            // search placement by a few degrees and a fraction of an angstrom, and the clash fraction moves
            // by a percent or two with it - and does not count, so that a group is CLASH only where no
            // orientation the rotation function offers packs in it.
            int refined = 0;
            for (size_t k = 0; k < by_tf.size() && refined < (fits ? RESCORED_PER_GROUP : 1); k++) {
                const size_t i = by_tf[k];
                auto placed = std::make_shared<gemmi::Structure>(
                        PlaceSearchModel(model, rots[i], tf[i].t, cell, sg));
                double clash = ClashFraction(*placed);
                PlacementScore score{NAN, NAN, NAN, NAN};
                gemmi::Mat33 moved;
                if (fits && clash <= 2 * CLASH_MAX) {
                    if (engines.get() == nullptr)
                        engines = make_engines(*placed, sg, *merge.merged, d_min, 1, logger);
                    score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads, &scorers,
                                             &moved);
                    clash = ClashFraction(*placed);
                    refined++;
                } else if (fits && have) {
                    continue;
                } else if (!fits) {
                    refined++;
                }
                // The refined placement with the lowest R-free; the best translation-function one where
                // none was refined.
                if (!have || (std::isfinite(score.r_free) && !(cand.r_free <= score.r_free))) {
                    have = true;
                    cand.rf_z = rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0;
                    cand.tf_z = tf[i].z;
                    cand.copy_tf_z = {tf[i].z};
                    cand.clash_fraction = clash;
                    cand.r_work = score.r_work;
                    cand.r_free = score.r_free;
                    cand.k_sol = score.k_sol;
                    cand.b_sol = score.b_sol;
                    cand.placed = placed;
                    cand.copies = 1;
                    orientation = rots[i];
                    copies = {*placed};
                    copy_rotation = {moved.multiply(rots[i])};
                    copy_position = {centroid(*placed)};
                }
                glog.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 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash,
                            score.r_work, score.r_free);
            }

            // --- further copies, one at a time, with the copies placed as a known partial structure ---
            // Each copy is refined on its own and all of them scored together.
            std::optional<TranslationSetup> setup2;
            // The first `count` copies refined in turn, each with all the others fixed, and the assembly scored.
            const auto refine_beside_others = [&](int count) {
                for (int c = 0; c < count; c++) {
                    std::vector<gemmi::Structure> rest;
                    for (size_t o = 0; o < copies.size(); o++)
                        if (static_cast<int>(o) != c)
                            rest.push_back(copies[o]);
                    const gemmi::Structure fixed = AssembleCopies(rest);
                    gemmi::Mat33 moved;
                    refine_rigid_body(copies[c], sg, *merge.merged, d_min, engines.get(), nthreads, &moved,
                                      &fixed.models[0]);
                    copy_rotation[c] = moved.multiply(copy_rotation[c]);
                    copy_position[c] = centroid(copies[c]);
                }
                const gemmi::Structure assembly = AssembleCopies(copies);
                const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol};
                const PlacementScore score = score_placement(assembly, sg, *merge.merged, d_min, nthreads, &scorers,
                                                             &held);
                cand.clash_fraction = ClashFraction(assembly, cand.copies);
                cand.r_work = score.r_work;
                cand.r_free = score.r_free;
                cand.k_sol = score.k_sol;
                cand.b_sol = score.b_sol;
                cand.placed = std::make_shared<gemmi::Structure>(assembly);
                glog.Info("Model prior:   {:<12} {} copies, refined beside each other: R-work {:.3f} R-free {:.3f}",
                            sg.xhm(), cand.copies, score.r_work, score.r_free);
            };
            while (fits && cand.copies < cand.copies_max && cand.clash_fraction <= CLASH_MAX
                   && std::isfinite(cand.r_free)) {
                // The copies placed are searched for again in their own orientations: copies related by a
                // translation alone (tNCS) share one.
                std::vector<gemmi::Mat33> rots2 = rots;
                rots2.insert(rots2.end(), copy_rotation.begin(), copy_rotation.end());
                if (!setup2)
                    setup2 = MakeTranslationSetup(refl, cell, sg, true);
                const std::vector<std::complex<double>> known = KnownStructureFactors(ft, *setup2, cell, copy_rotation,
                                                                                      copy_position);
                std::vector<TranslationPeak> tf2(rots2.size());
#ifdef JFJOCH_USE_CUDA
                if (gpu)
                    tf2 = gpu->TranslationFunction(*setup2, rots2, &known);
                else
#endif
                    ParallelFor(static_cast<int>(rots2.size()), nthreads, [&](int i) {
                        tf2[i] = TranslationFunction(ft, *setup2, rots2[i], &known);
                    });
                std::vector<size_t> by_tf2(rots2.size());
                std::iota(by_tf2.begin(), by_tf2.end(), 0);
                std::stable_sort(by_tf2.begin(), by_tf2.end(), [&](size_t a, size_t b) { return tf2[a].z > tf2[b].z; });

                struct Trial {
                    size_t i = 0;
                    double clash = 1;
                    double cc = 0;
                    PlacementScore score{NAN, NAN, NAN, NAN};
                    gemmi::Structure copy, assembly;
                    gemmi::Mat33 rotation;
                };
                std::optional<Trial> best;
                const int n = cand.copies + 1;
                // The highest peaks that do not sit on a copy already placed (or a mate of one), ranked by
                // the correlation of the assembly with E^2, and the best few of those through the rigid body.
                std::vector<Trial> trials;
                for (size_t k = 0; k < by_tf2.size() && trials.size() < COPY_TRIALS; k++) {
                    Trial trial;
                    trial.i = by_tf2[k];
                    trial.copy = PlaceSearchModel(model, rots2[trial.i], tf2[trial.i].t, cell, sg);
                    std::vector<gemmi::Structure> all = copies;
                    all.push_back(trial.copy);
                    trial.clash = ClashFraction(AssembleCopies(all), n);
                    if (trial.clash > 2 * CLASH_MAX)
                        continue;
                    trials.push_back(std::move(trial));
                }
                ParallelFor(static_cast<int>(trials.size()), nthreads, [&](int k) {
                    trials[k].cc = AssemblyCorrelation(ft, *setup2, cell, known, rots2[trials[k].i],
                                                       tf2[trials[k].i].t);
                });
                std::stable_sort(trials.begin(), trials.end(), [](const Trial &a, const Trial &b) { return a.cc > b.cc; });
                if (trials.size() > RESCORED_PER_GROUP)
                    trials.resize(RESCORED_PER_GROUP);
                const gemmi::Structure others = AssembleCopies(copies);
                for (Trial &trial : trials) {
                    std::vector<gemmi::Structure> all = copies;
                    gemmi::Mat33 moved;
                    refine_rigid_body(trial.copy, sg, *merge.merged, d_min, engines.get(), nthreads, &moved,
                                      &others.models[0]);
                    all.push_back(trial.copy);
                    trial.assembly = AssembleCopies(all);
                    trial.clash = ClashFraction(trial.assembly, n);
                    const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol};
                    trial.score = score_placement(trial.assembly, sg, *merge.merged, d_min, nthreads, &scorers, &held);
                    trial.rotation = moved.multiply(rots2[trial.i]);
                    glog.Info("Model prior:   {:<12} copy {} TF Z {:5.1f} CC {:.3f} clash {:4.1f}%  R-work {:.3f} R-free {:.3f}",
                                sg.xhm(), n, tf2[trial.i].z, trial.cc, 100 * trial.clash, trial.score.r_work,
                                trial.score.r_free);
                    // The first, in that order, that stands out, packs and lowers R-free.
                    if (tf2[trial.i].z >= TF_Z_PLACED && trial.clash <= CLASH_MAX
                        && trial.score.r_free <= cand.r_free - COPY_R_FREE_GAIN) {
                        best = trial;
                        break;
                    }
                }
                if (!best)
                    break;
                copies.push_back(best->copy);
                copy_rotation.push_back(best->rotation);
                copy_position.push_back(centroid(best->copy));
                cand.copies = n;
                cand.copy_tf_z.push_back(tf2[best->i].z);
                cand.tf_z = std::max(cand.tf_z, tf2[best->i].z);
                cand.clash_fraction = best->clash;
                cand.r_work = best->score.r_work;
                cand.r_free = best->score.r_free;
                cand.k_sol = best->score.k_sol;
                cand.b_sol = best->score.b_sol;
                cand.placed = std::make_shared<gemmi::Structure>(best->assembly);
                // The first copy, refined alone, refined again beside the second: next to a copy related to it
                // by a translation alone, one copy by itself fits the data no better than at random.
                if (n == 2)
                    refine_beside_others(1);
            }
            // Every copy but the last refined again beside all the others, which each was refined beside
            // only some of.
            if (cand.copies >= 3)
                refine_beside_others(cand.copies - 1);
            // The copies' R-factors were taken with the bulk solvent of the first copy alone held; all of
            // them scored with the solvent fitted to them.
            if (cand.copies >= 2) {
                const PlacementScore score = score_placement(*cand.placed, sg, *merge.merged, d_min, nthreads, &scorers);
                cand.r_work = score.r_work;
                cand.r_free = score.r_free;
                cand.k_sol = score.k_sol;
                cand.b_sol = score.b_sol;
            }
                        
            if (cand.vm < VM_MIN)
                cand.verdict = "MATTHEWS";
            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";
            out.copy_rotation = copy_rotation;
            out.cand = std::move(cand);
            out.orientation = orientation;
            return out;
        };
        // Their time in wall-clock seconds, counted as rescoring (the further copies' translation functions
        // included).
        const auto t_groups = std::chrono::steady_clock::now();
        std::vector<std::future<GroupResult>> running;
        for (size_t g = 0; g < groups.size(); g++)
            running.push_back(std::async(std::launch::async, search_group, g, std::vector<gemmi::Mat33>()));
        std::vector<GroupResult> found;
        for (auto &r : running)
            found.push_back(r.get());
        // Every group the best one leads searched again with the best one's copies' orientations tried
        // first: the groups share the rotations, so a copy the search placed in one is in the same
        // orientation in the others, and one group's search failing - next to copies related by a
        // translation, the first copy alone fits no better than at random, and the copy the search keeps
        // can be the wrong one - is then not taken for evidence against it. The better of the two searches
        // is kept.
        int lead = -1;
        for (size_t g = 0; g < found.size(); g++)
            if (found[g].cand.verdict == "PLACED" && (lead < 0 || found[g].cand.r_free < found[lead].cand.r_free))
                lead = static_cast<int>(g);
        if (lead >= 0) {
            std::vector<std::future<GroupResult>> again(found.size());
            for (size_t g = 0; g < found.size(); g++)
                if (static_cast<int>(g) != lead && found[g].cand.vm >= VM_MIN)
                    again[g] = std::async(std::launch::async, search_group, g, found[lead].copy_rotation);
            for (size_t g = 0; g < found.size(); g++)
                if (again[g].valid()) {
                    GroupResult r = again[g].get();
                    found[g].log.Info("Model prior:   {:<12} searched again from {}'s copies: {} copies, R-free {:.3f}, was {:.3f}",
                                      groups[g].xhm(), groups[lead].xhm(), r.cand.copies, r.cand.r_free,
                                      found[g].cand.r_free);
                    if (std::isfinite(r.cand.r_free) && !(found[g].cand.r_free <= r.cand.r_free)) {
                        r.log.ReplayInto(found[g].log);
                        r.log = found[g].log;
                        found[g] = std::move(r);
                    }
                }
        }
        result.t_rescore += seconds_since(t_groups);
        for (GroupResult &gr : found) {
            gr.log.ReplayInto(logger);
            orientations.push_back(gr.orientation);
            result.candidates.push_back(std::move(gr.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.
    // The null is built for the candidate the choice would fall on, and only for it: a placed candidate
    // that does not beat its own null is set aside, and the choice falls on the next.
    const auto choose_in = [&](int merge) -> int {
        int best = -1;
        while (true) {
            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;
            MRCandidate &c = result.candidates[best];
            const auto t_null = std::chrono::steady_clock::now();
            RigidBodyEngines engines = make_engines(*c.placed, c.space_group, *merges[merge].merged,
                                                    result.d_min_score, NULL_REPLICATES, logger);
            // One copy: its replicates placed by the translation function, as it was.
            std::function<std::vector<gemmi::Fractional>(const std::vector<gemmi::Mat33> &)> translate;
            if (c.copies == 1)
                translate = [&](const std::vector<gemmi::Mat33> &turned) {
                    const SearchReflections refl = MakeSearchReflections(*merges[merge].merged,
                                                                         gemmi::UnitCell(merges[merge].cell),
                                                                         merges[merge].space_group,
                                                                         result.d_min_search, result.d_max_search);
                    const TranslationSetup setup = MakeTranslationSetup(refl, gemmi::UnitCell(merges[merge].cell),
                                                                        c.space_group);
                    std::vector<TranslationPeak> tf(turned.size());
#ifdef JFJOCH_USE_CUDA
                    if (gpu)
                        tf = gpu->TranslationFunction(setup, turned);
                    else
#endif
                        ParallelFor(static_cast<int>(turned.size()), nthreads, [&](int i) {
                            tf[i] = TranslationFunction(ft, setup, turned[i]);
                        });
                    std::vector<gemmi::Fractional> out;
                    for (const TranslationPeak &p : tf)
                        out.push_back(p.t);
                    return out;
                };
            c.null_sigma = null_sigma(*c.placed, orientations[best], c.space_group, *merges[merge].merged,
                                      result.d_min_score, engines.get(), nthreads, &scorers, translate);
            result.t_null += seconds_since(t_null);
            if (c.null_sigma >= NULL_SIGMA_MIN)
                break;
            c.verdict = "NULL";
        }
        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;
            // With tNCS one copy alone fits no better than at random, so a group whose search stopped at
            // fewer copies than the chosen one's has not been shown unable to hold them.
            if (result.tncs.present && c.copies < result.candidates[best].copies
                && TncsRelatesGroups(c.space_group, result.candidates[best].space_group, result.tncs.vector)) {
                result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS",
                                              c.space_group.xhm(), result.candidates[best].space_group.xhm());
                return -2;
            }
            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;
}

void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector<MergedReflection> &merged,
                           const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads) {
    const MRCandidate &c = result.candidates[result.chosen];
    const auto t0 = std::chrono::steady_clock::now();
    Logger quiet = Logger::Buffered();
    ModelValidationResult v;
    try {
        v = ReadModelAnomalousMap(merged, *c.placed, data_space_group, c.k_sol, c.b_sol, quiet, nthreads);
    } catch (const std::runtime_error &e) {   // a GPU failure: no reading, so no relabelling
        logger.Warning("Model prior: the anomalous map could not be read ({})", e.what());
    }
    auto &h = result.hand_check;
    h.performed = true;
    h.model_fits = c.null_sigma >= NULL_SIGMA_MIN;
    h.fit_sigma = c.null_sigma;
    h.adopted = v.ok && h.model_fits && !v.anomalous_hands_disagree;
    h.anomalous_inverted = v.anomalous_hands_disagree;
    h.anomalous_pairs = v.anomalous_pairs;
    h.anomalous_scatterer_mean_sigma = v.anomalous_scatterer_mean_sigma;
    h.deepest_site = v.anomalous_deepest_site;
    h.deepest_sigma = v.anomalous_deepest_sigma;
    h.seconds = seconds_since(t0);
    if (!h.adopted) {
        result.decision = "ENANTIOMORPH_VETOED";
        logger.Info("Model prior: the enantiomorph {} is not taken into the answer - {} ({:.3f} s)",
                    c.space_group.xhm(),
                    h.anomalous_inverted ? fmt::format("the anomalous density at the placed model's atoms is "
                                                       "inverted ({} at {:.1f} sigma)", h.deepest_site, h.deepest_sigma)
                                         : fmt::format("the placed model does not fit against its null ({:+.1f} sd)",
                                                       h.fit_sigma),
                    h.seconds);
    }
}