Molecular replacement: a model's chains searched by themselves; tNCS half-vector moves told apart by a paired test

A model of several chains is searched for as a whole and then, side by side, as one chain of each
molecule in it (chains with the same residue at every number they share are one molecule; the most
complete one is searched), with the PAE domains' machinery. The whole model's choice stands unless a
chain's search chose a placement better by the R-free margin (0.03), or the whole model chose none; the
chains that chose must agree on the space group. MR_SEARCH_MODEL names the winner ("chain A, one of 2
alike"). An experimental dimer from another crystal form, placed as one rigid unit where the target
crystal holds the two copies at another separation (tNCS), was placed 3.4 A off per chain (R-free
0.47); the chain search places both copies (R-free 0.28 before refinement, refined 0.27), and where the
dimer is right the whole model is kept (R-free 0.26 against the chain's 0.29). MR takes about twice
as long on such a model (0.55 -> 1.1 s).

With tNCS by half a lattice vector, the placement and its moves by (v + o)/2 are now told apart by a
paired per-reflection test of |F_obs - F_calc| on the weak reflections (3 standard errors, as the
equal-copies group comparison) instead of a 0.01 margin in R over them; UNDECIDED only where the best
is not ahead of every other by that.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
(cherry picked from commit 690ee4784ddd8e046fec19edb558f28808d0a825)

docs/CPU_DATA_ANALYSIS_DECISIONS.md

@@ -314,8 +314,12 @@ Translational NCS changes what the search can tell apart (Read, Adams & McCoy, 2
 pseudo-translation, one copy beside a translated twin of itself fits the data as well in a group whose
 screw axes the translation turns into the chosen group's, so such a choice is left undecided unless the
 copies, carried into the other group and refined there, fit clearly worse; and where the vector is half
-a lattice vector, the placement is compared with every move of itself the vector allows, on the weak
-reflections it modulates.
+a lattice vector, the placement is compared with every move of itself the vector allows, each refined,
+on the weak reflections it modulates, and left undecided only where the best is not ahead of every other
+by 3 standard errors of the paired per-reflection difference there. A model of several chains is searched
+for as a whole and, side by side, as one chain of each molecule in it: the chains of an experimental model
+from another crystal form need not sit in this crystal as they sat in that one, and the whole model's
+choice gives way to a chain's where the chain's fits better by the R-free margin.
 
 The choice between groups and merges is a hypothesis comparison like §14.5: the candidate the choice
 falls on has to beat a null of the same model placed from random orientations, and lead the other

docs/RUGNUX_CLI.md

@@ -313,7 +313,7 @@ serial stills; `fft` (GPU) / `fftw` (CPU) index **de novo** and suit strong rota
 | `--free-flags <file>` | Take only the R-free test set of this MTZ or structure-factor mmCIF - e.g. the deposited reflections of the `--model` - for the validation, the refinement and the merged files; nothing else of the file is used, so unlike `-z` it fixes no space group, cell or frame. Taken on the model's axes where the merge matches the file's intensities (`R_FREE_FLAGS=`, see [Free set and model bias](RUGNUX_MODEL_AND_REFERENCE.md#free-set-and-model-bias)) |
 | `--reference-column <label>` | Reference column to use, by MTZ label — an mmCIF under the labels `gemmi cif2mtz` gives it (default: auto — F-model, else IMEAN/I/…, else FP/FOBS/F) |
 | `--model <file>` | Validate the merged intensities against this atomic model (PDB or mmCIF, gzipped or not; the format is taken from the file's content) — R-work / R-free and maps (see [Validating against a model](RUGNUX_ADVANCED.md#validating-against-a-model-rugnux---model)). An experimental model is placed as a rigid body, and by molecular replacement where the rigid body does not place it; a prediction (AlphaFold or ModelCIF) goes to molecular replacement straight away (`MODEL_PLACEMENT=`). It also settles the frame the reflections are written in: the enantiomorph, and the [indexing ambiguity](RUGNUX_ADVANCED.md#the-indexing-ambiguity) where no `-z` did; another space group or the twin hypothesis' subgroup it chooses is written beside the answer as `<prefix>_mr.mtz`/`.cif`. It does not change the processing: the merge is the same with and without it. A placed model is then [refined](RUGNUX_MODEL_AND_REFERENCE.md#refinement) against the merge it was placed on (`<prefix>_refined_model.cif`, `<prefix>_refined_maps.mtz`, `<prefix>_blobs.json`, `REFINE_` keys) - seconds on a GPU, typically 2-15 minutes on the CPU - which decides nothing about the processing |
-| `--mr <file>` | As `--model`, but placed by molecular replacement straight away, without trying the rigid body first: as many copies as the cell holds, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. A prediction is trimmed to the residues of pLDDT 70 and above, and where the AlphaFold DB's predicted-aligned-error JSON sits beside it under its own name it is cut into rigid domains; any other PDB or mmCIF is taken as an experimental structure, its polymer only (no waters, ligands, hydrogens or alternative conformers). It may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as `<prefix>_mr.mtz`/`.cif` with the placed model `<prefix>_mr_model.pdb`. `--mode mx` only; not with `--model` |
+| `--mr <file>` | As `--model`, but placed by molecular replacement straight away, without trying the rigid body first: as many copies as the cell holds, in every space group the merge cannot tell apart and in the twin hypothesis' subgroup. A prediction is trimmed to the residues of pLDDT 70 and above, and where the AlphaFold DB's predicted-aligned-error JSON sits beside it under its own name it is cut into rigid domains; any other PDB or mmCIF is taken as an experimental structure, its polymer only (no waters, ligands, hydrogens or alternative conformers). A model of several chains is searched for as a whole and, beside it, one chain of each molecule in it by itself, and the better placement kept (`MR_SEARCH_MODEL=`). It may relabel the enantiomorph of the answer, and a different space group it chooses is written beside the answer as `<prefix>_mr.mtz`/`.cif` with the placed model `<prefix>_mr_model.pdb`. `--mode mx` only; not with `--model` |
 | `--no-refine` | Do not refine the placed model (for debugging): by default a model that `--model` or `--mr` places is refined |
 | `--refine-solvent` | Also build the solvent of the refined model: waters added at difference-map peaks and removed where the maps no longer hold them, waters and ions relabelled, and each site audited (the `REFINE_SOLVENT_` keys and `<prefix>_solvent.json`). Off by default for now: the model is refined with its solvent as given, and the report says `REFINE_SOLVENT=NOT_PERFORMED`. Also with `--mode refine` |
 | `--no-refine-solvent` | The default: no solvent build or audit (turns off an earlier `--refine-solvent`) |

docs/RUGNUX_MODEL_AND_REFERENCE.md

@@ -246,7 +246,9 @@ The safeguards behind 2 and 3:
 - **Molecular replacement's own rules.** Its choice of group has its own null (`NULL_SIGMA`); a
   placement must pack (`CLASH`, `MATTHEWS`); a choice that translational NCS could explain, or one
   where the merge's own group was searched to fewer copies, is left `UNDECIDED`; a single domain of
-  a prediction cut by its PAE can only confirm the answer's group, never change it; in P1 the first
+  a prediction cut by its PAE can only confirm the answer's group, never change it; a model of several
+  chains is searched for as a whole and, beside it, one chain of each molecule by itself, and a chain
+  replaces the whole only where it fits better by 0.03 in R-free or the whole was placed nowhere; in P1 the first
   copy's origin is fixed and it counts only on its rotation function. The details are in
   [Molecular replacement](RUGNUX_REPORT.md#molecular-replacement).
 - **The wrong-model test.** Rugnux's validation battery runs a set of datasets once without a model

docs/RUGNUX_REPORT.md

@@ -1320,8 +1320,16 @@ Where the AlphaFold DB's predicted-aligned-error file sits beside the model unde
 into rigid domains by average linkage on the PAE (domains a mean 10 A apart, of 40 residues or more;
 `MR_PAE_DOMAINS=`; a prediction only). Where the whole model is placed nowhere, each domain is searched by itself, side by
 side; the domains placed have to agree, and a domain alone can only confirm the answer's group
-(`AGREES`), never change it (`UNDECIDED` otherwise). `MR_SEARCH_MODEL=` names what the decision
-was made with.
+(`AGREES`), never change it (`UNDECIDED` otherwise).
+
+A model of several chains is searched for as a whole and, beside it, one chain of each molecule in it
+by itself (the most complete of the chains with the same residue at every number they share), since the
+crystal may hold the molecules in another arrangement than the model has them in - another crystal form,
+or a tNCS the model's crystal did not have. The chain's search places its own copies (`COPIES`) by the
+same rules. The whole model's choice stands unless a chain's search chose a placement better by 0.03 in
+R-free, or the whole model chose none; the chains that chose have to agree on the space group.
+`MR_SEARCH_MODEL=` names what the decision was made with: `the whole model`, a chain (`chain A`, or
+`chain A, one of 2 alike`) or a PAE domain; the table is that search's.
 
 The table has one row per group and merge. A row is `PLACED` where one of its copies stands out, at most 5% of the C-alpha atoms clash with another
 copy or a symmetry mate (`CLASH`), and one copy fits the cell (`MATTHEWS`: Vm below 1.5 A^3/Da). The
@@ -1342,8 +1350,9 @@ in both. The group the merge was made in, searched to fewer copies than the chos
 the tNCS vector v is half a lattice vector, the reflections with h.v half-integral are weak, and a
 placement moved by (v + o) / 2 (o an origin shift of the group) fits the strong ones exactly as well:
 the candidate the choice falls on is compared with every such move of itself, refined and packing, by
-R-work and R-free together over the weak reflections alone; a move better by 0.01 there is taken
-instead, and where the two best are within 0.01 the choice is `UNDECIDED`. The rotation and translation functions, the rigid
+R-work and R-free together over the weak reflections alone: the best there is taken where its
+|F_obs - F_calc| on those reflections is lower than every other's by 3 standard errors of the paired
+per-reflection difference, and the choice is `UNDECIDED` only where it is not. The rotation and translation functions, the rigid
 body and the structure factors the R-factors are taken on run on the
 GPU where there is one (`MR_SEARCH_DEVICE=`). `MR_DECISION=` is then:
 

image_analysis/structure_refinement/MolecularReplacement.cpp

@@ -114,18 +114,19 @@ constexpr double NULL_SIGMA_MIN = 4.51;
 constexpr double ELLG_MIN = 10.0;
 // tNCS by half a lattice vector v: copies v apart make the reflections with h.v half-integral weak, and a
 // placement moved by (v + o) / 2, o an origin shift of the group, fits the strong ones as well as the right one
-// does; only the weak ones tell them apart. A placement is compared with every such move of itself on those
-// reflections (R_tncs_weak), and one of them that fits them better by TNCS_SHIFT_MARGIN is taken instead; where the
-// two best are closer than that, the placement is not decided on. A vector is taken as half-integral within
-// TNCS_HALF_TOLERANCE.
-constexpr double TNCS_SHIFT_MARGIN = 0.01;
+// does; only the weak ones tell them apart. A placement is compared with every such move of itself, each refined,
+// on those reflections (R_tncs_weak), and the best is taken where it fits them better than every other by
+// PAIRED_Z_MIN standard errors of the paired difference (paired_z); where it does not, the placement is not
+// decided on. A vector is taken as half-integral within TNCS_HALF_TOLERANCE.
+// The modulation of the intensities by tNCS following Read, Adams & McCoy (2013) Acta Cryst. D69, 176-183
+constexpr double TNCS_HALF_TOLERANCE = 0.05;
 // The chosen copies carried into a group with tNCS rule it out only where they fit the data themselves: a
 // placement at R-free CARRY_R_FREE_MAX or above, little better than a model at random, is not the crystal, and
 // how it reads in another group says nothing about that group.
 constexpr double CARRY_R_FREE_MAX = 0.50;
-constexpr double TNCS_HALF_TOLERANCE = 0.05;
 // Two placements with as many copies within R_FREE_MARGIN of each other: the one that fits every reflection
-// better by this many standard errors of the paired difference (paired_r_z) is taken.
+// better by this many standard errors of the paired difference (paired_r_z) is taken; the same for a placement
+// and its moves by half a tNCS vector, on the reflections the vector weakens.
 constexpr double PAIRED_Z_MIN = 3.0;
 
 double seconds_since(std::chrono::steady_clock::time_point t0) {
@@ -263,7 +264,8 @@ PlacementScore r_factors(const gemmi::AsuData<std::complex<float>> &fmodel, cons
         if (it == obs.end())
             continue;
         const double fo = it->second.f, fc = std::abs(hv.value);
-        if (residuals != nullptr)
+        // With a tNCS vector, of the reflections it weakens only.
+        if (residuals != nullptr && (tncs == nullptr || it->second.weak))
             (*residuals)[hkl_key(hv.hkl)] = {fo, std::fabs(fo - fc)};
         if (it->second.free) {
             num_f += std::fabs(fo - fc);
@@ -696,14 +698,10 @@ PlacementScore carry_copies(const gemmi::Structure &assembly, int copies, const
     return score_placement(AssembleCopies(body), sg, merged, d_min, nthreads, scorers);
 }
 
-// Two placements of one merge compared on every reflection to d_min, work and free together: how far b's
-// R lies below a's, in standard errors of that difference, the reflections paired.
-double paired_r_z(const gemmi::Structure &a, const gemmi::SpaceGroup &a_sg, const gemmi::Structure &b,
-                  const gemmi::SpaceGroup &b_sg, const std::vector<MergedReflection> &merged, double d_min,
-                  size_t nthreads, ScoreEngines *scorers) {
-    std::unordered_map<long, std::pair<double, double>> ra, rb;
-    score_placement(a, a_sg, merged, d_min, nthreads, scorers, nullptr, nullptr, nullptr, &ra);
-    score_placement(b, b_sg, merged, d_min, nthreads, scorers, nullptr, nullptr, nullptr, &rb);
+// How far the residuals |F_obs - F_calc| of b lie below those of a (score_placement's `residuals`), over the
+// reflections both have, in standard errors of their paired difference.
+double paired_z(const std::unordered_map<long, std::pair<double, double>> &ra,
+                const std::unordered_map<long, std::pair<double, double>> &rb) {
     double sum_fo = 0, sum = 0, sum2 = 0;
     size_t n = 0;
     for (const auto &[key, x] : ra) {
@@ -723,6 +721,17 @@ double paired_r_z(const gemmi::Structure &a, const gemmi::SpaceGroup &a_sg, cons
     return sum / std::sqrt(var * static_cast<double>(n));
 }
 
+// Two placements of one merge compared on every reflection to d_min, work and free together: how far b's
+// R lies below a's, in standard errors of that difference, the reflections paired.
+double paired_r_z(const gemmi::Structure &a, const gemmi::SpaceGroup &a_sg, const gemmi::Structure &b,
+                  const gemmi::SpaceGroup &b_sg, const std::vector<MergedReflection> &merged, double d_min,
+                  size_t nthreads, ScoreEngines *scorers) {
+    std::unordered_map<long, std::pair<double, double>> ra, rb;
+    score_placement(a, a_sg, merged, d_min, nthreads, scorers, nullptr, nullptr, nullptr, &ra);
+    score_placement(b, b_sg, merged, d_min, nthreads, scorers, nullptr, nullptr, nullptr, &rb);
+    return paired_z(ra, rb);
+}
+
 // 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]);
@@ -1085,6 +1094,78 @@ SearchModel DomainModel(const SearchModel &model, const std::vector<int> &residu
     return out;
 }
 
+// --- the chains of a model of several ---
+
+namespace {
+
+// Whether two chains are copies of one molecule: the same residue at every number both have, and those at least
+// half of the shorter chain's.
+bool same_molecule(const gemmi::Chain &a, const gemmi::Chain &b) {
+    std::map<std::string, std::string> names;
+    for (const gemmi::Residue &r : a.residues)
+        names[r.seqid.str()] = r.name;
+    size_t shared = 0;
+    for (const gemmi::Residue &r : b.residues) {
+        const auto it = names.find(r.seqid.str());
+        if (it == names.end())
+            continue;
+        if (it->second != r.name)
+            return false;
+        shared++;
+    }
+    return 2 * shared >= std::min(a.residues.size(), b.residues.size());
+}
+
+}  // namespace
+
+std::vector<SearchModel> ChainModels(const SearchModel &model) {
+    const std::vector<gemmi::Chain> &chains = model.st.models[0].chains;
+    if (chains.size() < 2)
+        return {};
+    // The chains grouped into molecules, each group in model order.
+    std::vector<std::vector<size_t>> molecules;
+    for (size_t c = 0; c < chains.size(); c++) {
+        auto it = std::find_if(molecules.begin(), molecules.end(), [&](const std::vector<size_t> &m) {
+            return same_molecule(chains[m[0]], chains[c]);
+        });
+        if (it == molecules.end())
+            molecules.push_back({c});
+        else
+            it->push_back(c);
+    }
+    std::vector<SearchModel> out;
+    for (const std::vector<size_t> &m : molecules) {
+        // The most complete chain of the molecule.
+        size_t keep = m[0];
+        for (size_t c : m)
+            if (chains[c].residues.size() > chains[keep].residues.size())
+                keep = c;
+        SearchModel one = model;
+        one.residue_index.clear();
+        one.residues_kept = 0;
+        std::vector<gemmi::Chain> kept;
+        size_t i = 0;
+        for (size_t c = 0; c < chains.size(); c++) {
+            const size_t n = chains[c].residues.size();
+            if (c == keep) {
+                kept.push_back(chains[c]);
+                one.residues_kept = static_cast<int>(n);
+                if (model.predicted)
+                    one.residue_index.assign(model.residue_index.begin() + i, model.residue_index.begin() + i + n);
+            }
+            i += n;
+        }
+        one.st.models[0].chains = std::move(kept);
+        finish_search_model(one);
+        if (!model.predicted)
+            one.error_A = model.error_A;   // not from the B-factors (PrepareSearchModel)
+        one.description = m.size() > 1 ? fmt::format("chain {}, one of {} alike", chains[keep].name, m.size())
+                                       : fmt::format("chain {}", chains[keep].name);
+        out.push_back(std::move(one));
+    }
+    return out;
+}
+
 // --- rotations ---
 
 std::vector<gemmi::Mat33> UniformRotations(double step_deg) {
@@ -3310,8 +3391,8 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
         // by side, each taken where the choice comes to it.
         std::vector<double> nulls(result.candidates.size(), NAN);
         bool first = true;
-        // With tNCS by half a lattice vector, candidate i against its moves by half that vector (TNCS_SHIFT_MARGIN):
-        // the better taken into it, and whether the two best are too close to tell apart.
+        // With tNCS by half a lattice vector, candidate i against its moves by half that vector (TNCS_HALF_TOLERANCE):
+        // the best taken into it, and whether it cannot be told from the others.
         // The tNCS vector is the answer's merge's, in its cell.
         const std::optional<gemmi::Fractional> half = merge == 0 ? half_tncs_vector(result.tncs) : std::nullopt;
         std::vector<char> shift_checked(result.candidates.size(), 0);
@@ -3325,10 +3406,13 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
                 std::shared_ptr<gemmi::Structure> placed;
                 PlacementScore score;
                 double clash = 0;
+                std::unordered_map<long, std::pair<double, double>> residuals;   // on the weak reflections
             };
-            std::vector<Option> options{{gemmi::Fractional(0, 0, 0), c.placed,
-                                         score_placement(*c.placed, c.space_group, merged, d_min, nthreads, &scorers,
-                                                         nullptr, nullptr, &v), c.clash_fraction}};
+            std::vector<Option> options(1);
+            options[0].placed = c.placed;
+            options[0].clash = c.clash_fraction;
+            options[0].score = score_placement(*c.placed, c.space_group, merged, d_min, nthreads, &scorers, nullptr,
+                                               nullptr, &v, &options[0].residuals);
             for (int o = 0; o < 8; o++) {
                 const gemmi::Fractional shift((v.x + 0.5 * (o & 1)) / 2, (v.y + 0.5 * ((o >> 1) & 1)) / 2,
                                               (v.z + 0.5 * ((o >> 2) & 1)) / 2);
@@ -3344,12 +3428,12 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
                 for (gemmi::Position &p : pos)
                     p = gemmi::Position(p + by);
                 SetModelPositions(moved->models[0], pos);
-                Option opt{shift, moved, {}, ClashFraction(*moved, c.copies)};
+                Option opt{shift, moved, {}, ClashFraction(*moved, c.copies), {}};
                 if (opt.clash <= 2 * CLASH_MAX) {
                     refine_rigid_body(*moved, c.space_group, merged, d_min, rigid_body_pool(merge), nthreads);
                     opt.clash = ClashFraction(*moved, c.copies);
                     opt.score = score_placement(*moved, c.space_group, merged, d_min, nthreads, &scorers, nullptr,
-                                                nullptr, &v);
+                                                nullptr, &v, &opt.residuals);
                 }
                 logger.Info("Molecular replacement: {} moved by {:.2f} {:.2f} {:.2f}: clash {:.1f}%, R-free {:.3f}, R on the "
                             "reflections tNCS weakens {:.3f} against {:.3f}", c.space_group.xhm(), shift.x, shift.y,
@@ -3364,12 +3448,26 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec
             std::stable_sort(packing.begin(), packing.end(), [](const Option *a, const Option *b) {
                 return a->score.r_tncs_weak < b->score.r_tncs_weak;
             });
-            if (packing.empty() || packing[0] == &options[0])
-                return packing.size() > 1 && packing[1]->score.r_tncs_weak - packing[0]->score.r_tncs_weak < TNCS_SHIFT_MARGIN;
+            if (packing.empty())
+                return false;
+            // The best on the weak reflections is taken where it fits them better than every other that packs by
+            // PAIRED_Z_MIN standard errors of the paired difference: the moves differ on those reflections alone, and
+            // with thousands of them a difference too small to read in R is still measured.
             const Option &b = *packing[0];
-            if (options[0].score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN
-                || (packing.size() > 1 && packing[1]->score.r_tncs_weak - b.score.r_tncs_weak < TNCS_SHIFT_MARGIN))
+            double z = INFINITY;
+            for (size_t k = 1; k < packing.size(); k++)
+                z = std::fmin(z, paired_z(packing[k]->residuals, b.residuals));
+            if (packing.size() > 1)
+                logger.Info("Molecular replacement: {} on the reflections tNCS weakens: {} fits them better than the {} "
+                            "other(s) that pack by {:.1f} standard errors or more", c.space_group.xhm(),
+                            &b == &options[0] ? std::string("the placement")
+                                              : fmt::format("the move by {:.2f} {:.2f} {:.2f}", b.shift.x, b.shift.y,
+                                                            b.shift.z),
+                            packing.size() - 1, z);
+            if (!(z >= PAIRED_Z_MIN))
                 return true;
+            if (&b == &options[0])
+                return false;
             logger.Info("Molecular replacement: {} taken moved by {:.2f} {:.2f} {:.2f}", c.space_group.xhm(), b.shift.x,
                         b.shift.y, b.shift.z);
             c.placed = b.placed;
@@ -3807,33 +3905,94 @@ MolecularReplacementResult run_molecular_replacement(const RugnuxModel &input, c
     result.model_name = input.name;
     result.pae_domains = static_cast<int>(domains.size());
     result.t_model += t_model;
+    // Search models besides the whole one, each on a thread of its own, its log kept and replayed in order, and with
+    // an equal share of the budget on every card.
+    const auto side_by_side = [&](const std::vector<SearchModel> &parts) {
+        std::vector<size_t> other_shares;
+        for (size_t b : other_budgets)
+            other_shares.push_back(b / parts.size());
+        std::vector<Logger> logs;
+        for (size_t i = 0; i < parts.size(); i++)
+            logs.push_back(Logger::Buffered());
+        std::vector<std::future<MolecularReplacementResult>> running;
+        for (size_t i = 0; i < parts.size(); i++)
+            running.push_back(std::async(std::launch::async, [&, i] {
+                return search_model(parts[i], merges, logs[i], nthreads, gpu_budget / parts.size(), other_shares);
+            }));
+        std::vector<MolecularReplacementResult> out;
+        for (size_t i = 0; i < parts.size(); i++) {
+            out.push_back(running[i].get());
+            logs[i].ReplayInto(logger);
+        }
+        return out;
+    };
+    // The search model the decision was made with taken over from another search, the whole model's times kept.
+    const auto take = [&](MolecularReplacementResult chosen, double seconds) {
+        chosen.model_name = input.name;
+        chosen.pae_domains = result.pae_domains;
+        chosen.t_model = result.t_model;
+        chosen.t_transform = result.t_transform;
+        chosen.t_rotation = result.t_rotation;
+        chosen.t_translation = result.t_translation;
+        chosen.t_rescore = result.t_rescore;
+        chosen.t_null = result.t_null;
+        chosen.t_domains = result.t_domains + seconds;
+        result = std::move(chosen);
+    };
+    // A model of several chains: one chain of each molecule in it searched by itself as well, beside the whole - the
+    // crystal may hold the molecules in another arrangement than the model has them in (another crystal form, a tNCS
+    // that is not the model's). The whole model's choice stands unless a chain's search chose a placement that fits
+    // better by R_FREE_MARGIN in R-free, or the whole model chose none; the chains that chose have to agree on the
+    // space group, and the one with the lowest R-free is taken.
+    const std::vector<SearchModel> chains = result.ok ? ChainModels(model) : std::vector<SearchModel>();
+    if (!chains.empty()) {
+        const auto t_chains = std::chrono::steady_clock::now();
+        logger.Info("Molecular replacement: {} chains in the model, {} molecule(s): each searched by itself as well",
+                    model.st.models[0].chains.size(), chains.size());
+        std::vector<MolecularReplacementResult> found = side_by_side(chains);
+        const double chain_seconds = seconds_since(t_chains);
+        const auto summary = [](const MolecularReplacementResult &r) {
+            if (r.chosen < 0)
+                return r.ok ? r.decision : r.failure_reason;
+            const MRCandidate &c = r.candidates[r.chosen];
+            return fmt::format("{} {}, {} {}, R-free {:.3f}", r.decision, c.space_group.xhm(), c.copies,
+                               c.copies == 1 ? "copy" : "copies", c.r_free);
+        };
+        logger.Info("Molecular replacement:   {}: {}", result.search_model, summary(result));
+        int best = -1;
+        bool agree = true;
+        for (size_t i = 0; i < found.size(); i++) {
+            logger.Info("Molecular replacement:   {}: {}", found[i].search_model, summary(found[i]));
+            if (!found[i].ok || found[i].chosen < 0)
+                continue;
+            const MRCandidate &c = found[i].candidates[found[i].chosen];
+            if (best >= 0 && c.space_group.xhm() != found[best].candidates[found[best].chosen].space_group.xhm())
+                agree = false;
+            if (best < 0 || c.r_free < found[best].candidates[found[best].chosen].r_free)
+                best = static_cast<int>(i);
+        }
+        const double whole_r_free = result.chosen >= 0 ? result.candidates[result.chosen].r_free : INFINITY;
+        if (best >= 0 && agree && found[best].candidates[found[best].chosen].r_free <= whole_r_free - R_FREE_MARGIN) {
+            take(std::move(found[best]), chain_seconds);
+        } else {
+            result.t_domains += chain_seconds;
+        }
+        logger.Info("Molecular replacement: decided with {}{}", result.search_model,
+                    best >= 0 && !agree ? " (the chains that chose disagree on the space group)" : "");
+    }
     // Where the whole model is placed nowhere and its prediction falls apart into domains whose relative
     // placement it is not confident of, each domain is searched by itself; the domains that are placed have
     // to agree on the space group, and the one with the lowest R-free is taken where that is the answer's
     // own.
     if (result.ok && result.chosen < 0 && result.decision.starts_with("NOT_PLACED") && domains.size() >= 2) {
         const auto t_domains = std::chrono::steady_clock::now();
-        // Each domain on a thread of its own, its log kept and replayed in order, and with an equal share of the
-        // budget on every card.
-        std::vector<size_t> other_shares;
-        for (size_t b : other_budgets)
-            other_shares.push_back(b / domains.size());
-        std::vector<Logger> logs;
-        for (size_t i = 0; i < domains.size(); i++)
-            logs.push_back(Logger::Buffered());
-        std::vector<std::future<MolecularReplacementResult>> running;
-        for (size_t i = 0; i < domains.size(); i++)
-            running.push_back(std::async(std::launch::async, [&, i] {
-                return search_model(DomainModel(model, domains[i]), merges, logs[i], nthreads,
-                                    gpu_budget / domains.size(), other_shares);
-            }));
+        std::vector<SearchModel> parts;
+        for (const std::vector<int> &d : domains)
+            parts.push_back(DomainModel(model, d));
         std::vector<MolecularReplacementResult> found;
-        for (size_t i = 0; i < domains.size(); i++) {
-            MolecularReplacementResult r = running[i].get();
-            logs[i].ReplayInto(logger);
+        for (MolecularReplacementResult &r : side_by_side(parts))
             if (r.ok && (r.chosen >= 0 || r.decision.starts_with("UNDECIDED")))
                 found.push_back(std::move(r));
-        }
         const double domain_seconds = seconds_since(t_domains);
         int best = -1;
         bool agree = true;
@@ -3853,19 +4012,9 @@ MolecularReplacementResult run_molecular_replacement(const RugnuxModel &input, c
             result.decision = fmt::format("UNDECIDED: a PAE domain alone chooses {}",
                                           found[best].candidates[found[best].chosen].space_group.xhm());
         } else if (best >= 0 && agree) {
-            MolecularReplacementResult chosen = std::move(found[best]);
-            chosen.model_name = input.name;
-            chosen.pae_domains = result.pae_domains;
-            chosen.t_model = result.t_model;
-            chosen.t_transform = result.t_transform;
-            chosen.t_rotation = result.t_rotation;
-            chosen.t_translation = result.t_translation;
-            chosen.t_rescore = result.t_rescore;
-            chosen.t_null = result.t_null;
-            result = std::move(chosen);
-            result.t_domains = domain_seconds;
+            take(std::move(found[best]), domain_seconds);
         } else {
-            result.t_domains = domain_seconds;
+            result.t_domains += domain_seconds;
             if (!found.empty())
                 result.decision = "UNDECIDED: the PAE domains placed disagree";
         }

image_analysis/structure_refinement/MolecularReplacement.h

@@ -67,6 +67,11 @@ std::vector<std::vector<int>> PaeDomains(const std::vector<std::vector<double>>
 // The search model of one domain: the residues of `model` at those positions, recentred.
 SearchModel DomainModel(const SearchModel &model, const std::vector<int> &residues);
 
+// The search models of a model of several chains: one chain of each molecule in it, the most complete, recentred
+// ("chain A"). Two chains are one molecule where they have the same residue at every number both have, and those
+// are at least half of the shorter chain's. Empty for a model of one chain.
+std::vector<SearchModel> ChainModels(const SearchModel &model);
+
 // Rotations covering SO(3) with neighbours about `step_deg` apart: the identity, then shells of
 // rotation angle kappa = step, 2 step, ... 180 deg, each with its axes on a golden-spiral lattice of
 // the sphere as dense as 2 sin(kappa/2) - the distance a change of axis moves a rotation - needs.
@@ -348,7 +353,7 @@ struct MolecularReplacementResult {
     std::string model_name;                          // RugnuxModel::name
     int residues_read = 0, residues_kept = 0;
     double mass_da = 0;
-    std::string search_model;                        // the whole model, or the PAE domain the decision was made with
+    std::string search_model;                        // the whole model, or the chain or PAE domain the decision was made with
     bool predicted_model = true;                     // SearchModel::predicted
     int pae_domains = 0;                             // the rigid domains of the prediction's PAE; 0 without a PAE file
     double d_min_search = 0, d_max_search = 0;
@@ -376,7 +381,7 @@ struct MolecularReplacementResult {
     } hand_check;
     // Seconds per step.
     double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0;
-    double t_domains = 0;   // the PAE domains' own searches, wall clock (they run side by side)
+    double t_domains = 0;   // the chains' and the PAE domains' own searches, wall clock (each kind side by side)
     // Searched on another merge of the crystal and carried over to this one (CarryMolecularReplacement): the
     // chosen placement's R-free on the merge it was found on, NaN where this merge was searched itself.
     double r_free_searched = NAN;

tests/MolecularReplacementTest.cpp

@@ -526,6 +526,57 @@ TEST_CASE("MolecularReplacement_PaeDomains", "[MolecularReplacement]") {
     CHECK(d.description == "PAE domain of 50 residues (1-50)");
 }
 
+// A model of three chains: A and B one molecule (B one residue shorter, at the start), C another. One chain of each
+// molecule, the most complete, recentred and keeping the experimental structure's coordinate error.
+TEST_CASE("MolecularReplacement_ChainModels", "[MolecularReplacement]") {
+    const auto chain = [](const std::string &name, int first, int last, const std::vector<std::string> &names,
+                          double x0) {
+        gemmi::Chain ch(name);
+        for (int i = first; i <= last; i++) {
+            gemmi::Residue r;
+            r.name = names[i % names.size()];
+            r.seqid = gemmi::SeqId(i, ' ');
+            gemmi::Atom a;
+            a.name = "CA";
+            a.element = gemmi::El::C;
+            a.pos = gemmi::Position(x0 + 3.8 * i, 0.0, 0.0);
+            a.b_iso = 30.0f;
+            a.occ = 1.0f;
+            r.atoms.push_back(a);
+            ch.residues.push_back(r);
+        }
+        return ch;
+    };
+    SearchModel m;
+    m.predicted = false;
+    m.error_A = 0.67;
+    gemmi::Model model(1);
+    model.chains.push_back(chain("A", 1, 30, {"ALA", "GLY", "LYS"}, 0.0));
+    model.chains.push_back(chain("B", 2, 30, {"ALA", "GLY", "LYS"}, 200.0));
+    model.chains.push_back(chain("C", 1, 20, {"SER"}, 400.0));
+    m.st.models.push_back(model);
+    m.residues_read = m.residues_kept = 79;
+
+    const std::vector<SearchModel> parts = ChainModels(m);
+    REQUIRE(parts.size() == 2);
+    CHECK(parts[0].description == "chain A, one of 2 alike");
+    CHECK(parts[0].residues_kept == 30);
+    CHECK(parts[0].residues_read == 79);
+    CHECK(parts[0].mass_da == Catch::Approx(30 * 110.0));
+    CHECK(parts[0].error_A == Catch::Approx(0.67));
+    CHECK(parts[0].st.models[0].chains.size() == 1);
+    CHECK(parts[0].st.models[0].chains[0].residues.front().atoms[0].pos.x == Catch::Approx(-3.8 * 14.5));
+    CHECK(parts[1].description == "chain C");
+    CHECK(parts[1].residues_kept == 20);
+
+    // The same residue numbers with another residue at one of them: another molecule.
+    m.st.models[0].chains[1].residues[5].name = "TRP";
+    CHECK(ChainModels(m).size() == 3);
+    // One chain: nothing besides the whole model.
+    m.st.models[0].chains.resize(1);
+    CHECK(ChainModels(m).empty());
+}
+
 // A placement of a group carried into a subgroup's merge: in P 21 21 21 read as P 1 21 1 on the same cell each
 // copy becomes two, the second its image under a two-fold the subgroup does not have; in C 2 2 21 read as
 // P 1 21 1 on the primitive monoclinic cell half the volume holds the same two.