From 60a8bf3981be283dd75516f3f2f75f79def9206c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 11 Oct 2026 13:35:05 +0200 Subject: [PATCH] rugnux --af-model (experimental): single-copy molecular replacement of a predicted model A predicted model (AlphaFold mmCIF, pLDDT in the B column) is trimmed at pLDDT 70, given B-factors from the pLDDT (Oeffner et al. 2022) and placed by molecular replacement, one copy in the ASU: - rotation function on a uniform 6 deg sampling of the point group's rotation ASU, each orientation scored by the correlation of E^2 with the symmetry-summed |F|^2 read off the model's molecular transform, computed once by FFT and interpolated (Navaza 1994); - full-symmetry Crowther & Blow translation function (one FFT per orientation and group, through MapFromFPhi/FFTW) for the 50 best orientations in every group with the merge's rotations; - the best three placements per group through RefineRigidBody (GPU where there is one), bulk solvent and R-free to 3 A, a C-alpha packing check against symmetry mates, Matthews, and the random-orientation null for a placed candidate. Decision: PLACED = TF Z >= 8, clash <= 5%, Vm in [1.5, 4), null >= 3 sd. Within the answer's merge the placed group with the lowest R-free is chosen where it leads every other packing group by 0.03; the twin hypothesis' merge is searched, and chosen, only where no group of the answer's merge can hold the model at all (Matthews, packing). An enantiomorph relabels the answer as --model does; another group goes to _mr.mtz/.cif plus _mr_model.pdb; the MODEL PRIOR report section carries the table. The search runs on the CPU in both builds; without --af-model nothing changes. Measured (RTX 5080, MR step): 3r6o twin I41 chosen over I4122 (Matthews + clash) in 1.8 s, solution 0.09 A CA rmsd from Phaser's; 11if P41 -> P43 1.2 s; 7ris P3121 -> P3221 1.7 s; 9zmu P6122 -> P6522 4.7 s; 11if's model against 3r6o's data NOT_PLACED. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/ACKNOWLEDGEMENT.md | 29 + docs/RUGNUX_REPORT.md | 34 +- .../structure_refinement/CMakeLists.txt | 1 + .../MolecularReplacement.cpp | 994 ++++++++++++++++++ .../MolecularReplacement.h | 165 +++ rugnux/ResultReport.cpp | 75 ++ rugnux/Rugnux.h | 12 + rugnux/RugnuxCommandLine.cpp | 2 + rugnux/RugnuxScaleMerge.cpp | 84 ++ rugnux/rugnux_cli.cpp | 13 + tests/CMakeLists.txt | 1 + tests/MolecularReplacementTest.cpp | 166 +++ 12 files changed, 1575 insertions(+), 1 deletion(-) create mode 100644 image_analysis/structure_refinement/MolecularReplacement.cpp create mode 100644 image_analysis/structure_refinement/MolecularReplacement.h create mode 100644 tests/MolecularReplacementTest.cpp diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index f75778070..bc92f289d 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -438,6 +438,35 @@ Problems Whose Variables Separate" (1973), SIAM J. Numer. Anal. 10, 413-432 for solving separable nonlinear least squares problems" (1975), BIT 15, 49-57 [doi:10.1007/BF01932995](https://doi.org/10.1007/BF01932995). +**Predicted models** — `--af-model` (experimental) places an AlphaFold prediction by molecular +replacement. J. Jumper, R. Evans, A. Pritzel et al., "Highly accurate protein structure prediction +with AlphaFold" (2021), Nature 596, 583-589 +[doi:10.1038/s41586-021-03819-2](https://doi.org/10.1038/s41586-021-03819-2); the models are those of +the AlphaFold Protein Structure Database: M. Varadi, S. Anyango, M. Deshpande et al., "AlphaFold +Protein Structure Database: massively expanding the structural coverage of protein-sequence space +with high-accuracy models" (2022), Nucleic Acids Res. 50, D439-D444 +[doi:10.1093/nar/gkab1061](https://doi.org/10.1093/nar/gkab1061). The model is prepared as +phenix.process_predicted_model prepares one - residues below pLDDT 70 removed, B-factors from the +coordinate error the pLDDT predicts: R. D. Oeffner, T. Croll, C. Millán et al., "Putting AlphaFold +models to work with phenix.process_predicted_model and ISOLDE" (2022), Acta Cryst. D78, 1303-1314 +[doi:10.1107/S2059798322010026](https://doi.org/10.1107/S2059798322010026). + +**Molecular replacement** — the search scores each orientation of the model against the data by +interpolating its molecular transform, computed once by FFT, as Navaza's AMoRe does: J. Navaza, +"AMoRe: an automated package for molecular replacement" (1994), Acta Cryst. A50, 157-163 +[doi:10.1107/S0108767393007597](https://doi.org/10.1107/S0108767393007597). The translation +function is the full-symmetry one, all pairs of symmetry copies on one Fourier synthesis: +R. A. Crowther and D. M. Blow, "A method of positioning a known molecule in an unknown crystal +structure" (1967), Acta Cryst. 23, 544-548 +[doi:10.1107/S0365110X67003172](https://doi.org/10.1107/S0365110X67003172). A placement is called +a solution, as Phaser calls one, where its translation-function Z-score reaches 8: A. J. McCoy, +R. W. Grosse-Kunstleve, P. D. Adams, M. D. Winn, L. C. Storoni and R. J. Read, "Phaser +crystallographic software" (2007), J. Appl. Cryst. 40, 658-674 +[doi:10.1107/S0021889807021206](https://doi.org/10.1107/S0021889807021206). Whether one copy fits +the cell is read off the Matthews coefficient: B. W. Matthews, "Solvent content of protein +crystals" (1968), J. Mol. Biol. 33, 491-497 +[doi:10.1016/0022-2836(68)90205-2](https://doi.org/10.1016/0022-2836(68)90205-2). + ## Software and computing methods Decoding bitshuffle+LZ4 images on the GPU, rather than decompressing them on the host and uploading diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index e7a6d5fbf..516d1c699 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -109,7 +109,8 @@ The header block above section 1 records **how the result was produced**: `RUGNU Rates, per-image costs and progress remain on stdout only. Sections, in order: the `SUMMARY`, then `1. DATA SET AND GEOMETRY`, `2. CRYSTAL`, -`3. MERGED DATA`, `4. DIAGNOSTICS`, and `5. MODEL VALIDATION` only with `--model` +`3. MERGED DATA`, `4. DIAGNOSTICS`, `5. MODEL VALIDATION` only with `--model`, and the unnumbered +`MODEL PRIOR (EXPERIMENTAL)` only with `--af-model` ([below](#model-prior---af-model-experimental)) (see [The summary](#the-summary-and-what-the-run-decided) below). The numbering is contiguous, and the fifth section appearing renumbers nothing; a stage that did not run states that inside its section — `MERGE= NOT_PERFORMED` — rather than the section disappearing. @@ -958,3 +959,34 @@ the scaling has mistreated. files are **byte for byte** what a run with no model would have written — same space group, same indexing, same `.mtz`, `.cif`, `.hkl` and `_unmerged.mtz`. A rejected model is therefore safe to try: it costs the null's compute and changes nothing else. + +## Model prior (`--af-model`, experimental) + +`--af-model ` takes a predicted model - an AlphaFold mmCIF, with the pLDDT in the B-factor +column - and places it by molecular replacement: residues below pLDDT 70 removed, B-factors from the +pLDDT, then a rotation function over the point group (6 deg sampling, 4-12 A), a translation function +for the 50 best orientations in every space group the merge cannot tell apart, and the best three +placements per group through the rigid body, bulk solvent and R-free (to `SCORE_D_MIN`, 3 A or the +data's limit). One copy in the asymmetric unit only, for now. + +The table has one row per group and merge. A row is `PLACED` where its translation-function Z-score +(`TF_Z`) reaches 8, at most 5% of its C-alpha atoms clash with symmetry mates (`CLASH`), one copy fits +the cell (`MATTHEWS`: Vm below 1.5 A^3/Da; `MULTI_COPY`: Vm of 4 or more, which is reported but not +supported yet) and its R-work stands at least 3 sd above the same model refitted from random +orientations (`NULL_SIGMA`, or `NULL` where it does not). `MODEL_PRIOR_DECISION=` is then: + +- `AGREES` - the placed group with the lowest R-free is the answer's, ahead of every other group of + the merge that packs by 0.03 in R-free; +- `ENANTIOMORPH` - it is the answer's enantiomorph, and the answer's files take that label (a label + only, as with `--model`; no reflection moves); +- `SPACE_GROUP` - another group with the answer's rotations; written as `_mr.mtz`/`.cif` from + the answer's merge, the reflections that group extinguishes left out; +- `TWIN_SUBGROUP` - no group of the answer's merge can hold the model at all (`MATTHEWS`/`CLASH` in + every row) and a group of the twin hypothesis' subgroup can; written as `_mr.mtz`/`.cif` + from the twin hypothesis' merge, with its law and fraction. The twin hypothesis' merge is searched + only in this case; +- `UNDECIDED` or `NOT_PLACED` - nothing is written. + +With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `_mr_model.pdb`/`.cif`. +Apart from the enantiomorph label, the answer's own files are never changed by the model. + diff --git a/image_analysis/structure_refinement/CMakeLists.txt b/image_analysis/structure_refinement/CMakeLists.txt index 4f6c6b3d1..9dcef0e26 100644 --- a/image_analysis/structure_refinement/CMakeLists.txt +++ b/image_analysis/structure_refinement/CMakeLists.txt @@ -2,6 +2,7 @@ ADD_LIBRARY(JFJochStructureRefinement STATIC ModelFFT.cpp ModelFFT.h ModelGrid.cpp ModelGrid.h ModelScaling.cpp ModelScaling.h + MolecularReplacement.cpp MolecularReplacement.h ModelValidation.cpp ModelValidation.h RigidBodyRefine.cpp RigidBodyRefine.h $<$:RigidBodyGPU.cpp> diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp new file mode 100644 index 000000000..c57fe5eb9 --- /dev/null +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -0,0 +1,994 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "MolecularReplacement.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // read_structure_gz +#include // remove_hydrogens +#include // DensityCalculator +#include // IT92 x-ray form factors +#include // Scaling (bulk solvent + anisotropic B) + +#include "../../common/JFJochMath.h" // PI +#include "../../common/Logger.h" +#include "../../common/ParallelFor.h" +#include "ModelFFT.h" +#include "ModelGrid.h" +#include "ModelScaling.h" +#include "RigidBodyRefine.h" +#ifdef JFJOCH_USE_CUDA +#include "RigidBodyGPU.h" +#endif + +namespace { + +using Table = gemmi::IT92; + +// The search's resolution range: the low-resolution end is where the bulk solvent, which the search +// does not model, stops dominating; the high-resolution end is where a predicted model still agrees +// with the crystal's structure in its details. +constexpr double SEARCH_D_MIN = 4.0; +constexpr double SEARCH_D_MAX = 12.0; +constexpr int SEARCH_BINS = 10; +constexpr double ROTATION_STEP_DEG = 6.0; +// Orientations the translation function is run for, per merge, and placements per space group taken +// through the rigid body. +constexpr int TRANSLATION_ROTATIONS = 50; +constexpr int RESCORED_PER_GROUP = 3; +// The resolution the placements are refined and scored to: whether a model sits in the right place +// is a low-resolution question, the rigid body's ladder ends at 3.5 A anyway, and every R-factor of +// the comparison is made to the same limit. +constexpr double SCORE_D_MIN = 3.0; + +// The decision. A translation-function peak this many standard deviations above the function's mean +// is a solution, where the model packs and one copy fits the cell; the McCoy et al. (2007) J. Appl. +// Cryst. 40, 658-674 convention for a single-copy search, taken over to this Z-score's own null on +// the Crowther & Blow function (measured: the highest peak of a wrong orientation reads 4-6). +constexpr double TF_Z_PLACED = 8.0; +constexpr double CLASH_MAX = 0.05; +// Matthews (1968) J. Mol. Biol. 33, 491-497: protein crystals have Vm of about 1.7 to 3.5 A^3/Da. +// Below VM_MIN one copy does not fit the cell; at or above VM_TWO_COPIES two copies would. +constexpr double VM_MIN = 1.5; +constexpr double VM_TWO_COPIES = 4.0; +constexpr double DALTON_PER_RESIDUE = 110.0; +// How far below every other candidate of its merge the chosen one's R-free must be. +constexpr double R_FREE_MARGIN = 0.03; +// The null (as in model validation): the same model refitted from random orientations about its +// centroid, at a coarse limit, and how far above that the real placement's R-work has to sit. +constexpr int NULL_REPLICATES = 9; +constexpr unsigned NULL_SEED = 20261011; +constexpr double NULL_D_MIN = 3.5; +constexpr double NULL_SIGMA_MIN = 3.0; + +double seconds_since(std::chrono::steady_clock::time_point t0) { + return std::chrono::duration(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 u(0.0, 1.0); + const double u1 = u(rng), t2 = 2 * PI * u(rng), t3 = 2 * PI * u(rng); + const double r1 = std::sqrt(1 - u1), r2 = std::sqrt(u1); + const double x = r1 * std::sin(t2), y = r1 * std::cos(t2), z = r2 * std::sin(t3), w = r2 * std::cos(t3); + return {1 - 2 * (y * y + z * z), 2 * (x * y - z * w), 2 * (x * z + y * w), + 2 * (x * y + z * w), 1 - 2 * (x * x + z * z), 2 * (y * z - x * w), + 2 * (x * z - y * w), 2 * (y * z + x * w), 1 - 2 * (x * x + y * y)}; +} + +double trace(const gemmi::Mat33 &m) { + return m.a[0][0] + m.a[1][1] + m.a[2][2]; +} + +gemmi::Vec3 miller_vec(const std::array &h) { + return {static_cast(h[0]), static_cast(h[1]), static_cast(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(h[0]) + 512) * 1048576L + (static_cast(h[1]) + 512) * 1024L + (h[2] + 512); +} + +int good_even_size(int n) { + while (n % 2 != 0 || !gemmi::has_small_factorization(n)) + ++n; + return n; +} + +// --- placements scored against the data: R-work / R-free after the bulk solvent and the scale --- + +struct PlacementScore { + double r_work = 1, r_free = 1; +}; + +gemmi::AsuData> working_amplitudes(const std::vector &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> fobs; + fobs.unit_cell_ = cell; + fobs.spacegroup_ = &sg; + for (const MergedReflection &r : merged) { + if (std::isnan(r.F) || r.rfree_flag || r.d < d_min) + continue; + gemmi::Miller h{{r.h, r.k, r.l}}; + if (!asu.is_in(h)) + h = asu.to_asu(h, gops).first; + fobs.v.push_back({h, {r.F, 1.0f}}); + } + fobs.ensure_asu(); + fobs.ensure_sorted(); + return fobs; +} + +// The model's Fcalc and bulk-solvent mask to d_min, scaled to the working amplitudes as model +// validation scales them (ModelValidation.cpp, fit_model), and the R-factors over the reflections +// to d_min. +PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, size_t nthreads) { + gemmi::DensityCalculator 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 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 leaders = OrbitLeaders(mask, nthreads); + PutMaskOnGrid(mask, st.models[0], leaders, nthreads); + gemmi::AsuData> fmask = MapToFPhi(mask).prepare_asu_data(d_min, 0); + PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads); + gemmi::AsuData> fmodel = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false); + + gemmi::AsuData> fobs = working_amplitudes(merged, st.cell, sg, d_min); + gemmi::Scaling scaling(st.cell, &sg); + scaling.use_solvent = true; + scaling.prepare_points(fmodel, fobs, &fmask); + FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads); + scaling.scale_data(fmodel, &fmask); + + const gemmi::GroupOps gops = sg.operations(); + const gemmi::ReciprocalAsu asu(&sg); + std::unordered_map> obs; + for (const MergedReflection &r : merged) { + if (std::isnan(r.F) || r.d < d_min) + continue; + gemmi::Miller h{{r.h, r.k, r.l}}; + if (!asu.is_in(h)) + h = asu.to_asu(h, gops).first; + obs[hkl_key(h)] = {r.F, r.rfree_flag}; + } + double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + for (const auto &hv : fmodel.v) { + const auto it = obs.find(hkl_key(hv.hkl)); + if (it == obs.end()) + continue; + const double fo = it->second.first, fc = std::abs(hv.value); + if (it->second.second) { + num_f += std::fabs(fo - fc); + den_f += fo; + } else { + num_w += std::fabs(fo - fc); + den_w += fo; + } + } + PlacementScore out; + out.r_work = den_w > 0 ? num_w / den_w : 1; + out.r_free = den_f > 0 ? num_f / den_f : 1; + return out; +} + +// The rigid-body engines for one model in one group, where there is a GPU. +struct RigidBodyEngines { +#ifdef JFJOCH_USE_CUDA + std::unique_ptr 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 &merged, double d_min, size_t engines, + Logger &logger) { + RigidBodyEngines out; +#ifdef JFJOCH_USE_CUDA + const double finest = RigidBodyLadder(d_min).back(); + size_t observations = 0; + for (const MergedReflection &r : merged) + if (!std::isnan(r.F) && r.d >= finest) + ++observations; + out.pool = RigidBodyGPUPool::Create(st.models[0], st.cell, sg, d_min, observations, engines, logger); +#else + (void) st; (void) sg; (void) merged; (void) d_min; (void) engines; (void) logger; +#endif + return out; +} + +// The placement through the rigid body, and its R-factors where the rigid body left it. A search +// placement is a few degrees and a grid step from its answer, well inside the rigid body's reach, so +// the step is taken without the R-free check model validation makes before committing one. +PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, RigidBodyGPUPool *pool, + size_t nthreads) { + const std::vector start = ModelPositions(st.models[0]); + Logger quiet = Logger::Buffered(); + const RigidBodyRefineResult rb = RefineRigidBody(st.models[0], st.cell, sg, + working_amplitudes(merged, st.cell, sg, d_min), d_min, quiet, + nthreads, pool); + if (!rb.converged) + SetModelPositions(st.models[0], start); + return score_placement(st, sg, merged, d_min, nthreads); +} + +// The placed model's R-work at NULL_D_MIN against the same model refitted from random orientations +// about its centroid, in standard deviations of the latter: the null model validation builds for a +// model that claims something (ModelValidation.cpp). An orientation the crystal cannot tell from the +// placement's - within the rigid body's reach of one of its symmetry classes - is drawn again. +double null_sigma(const gemmi::Structure &placed, const gemmi::Mat33 &orientation, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, RigidBodyGPUPool *pool, + size_t nthreads) { + const double null_d_min = std::max(d_min, NULL_D_MIN); + std::vector coarse; + for (const MergedReflection &r : merged) + if (r.d >= null_d_min) + coarse.push_back(r); + const double real = score_placement(placed, sg, coarse, null_d_min, nthreads).r_work; + + const std::vector pg = CartesianPointGroup(sg, placed.cell); + const double reach = RigidBodyReachDeg(placed.models[0], null_d_min); + std::vector draws; + std::mt19937 rng(NULL_SEED); + while (static_cast(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 r_work(NULL_REPLICATES); + std::atomic next{0}; + std::vector> workers; + for (size_t t = 0; t < std::min(std::max(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 pos = ModelPositions(rep.models[0]); + gemmi::Vec3 c; + for (const gemmi::Position &p : pos) + c += p; + c *= 1.0 / static_cast(pos.size()); + for (gemmi::Position &p : pos) + p = gemmi::Position(draws[i].multiply(gemmi::Vec3(p) - c) + c); + SetModelPositions(rep.models[0], pos); + r_work[i] = refine_placement(rep, sg, coarse, null_d_min, pool, nthreads).r_work; + } + })); + for (auto &w : workers) + w.get(); + const double mean = std::accumulate(r_work.begin(), r_work.end(), 0.0) / NULL_REPLICATES; + double s2 = 0; + for (double x : r_work) + s2 += (x - mean) * (x - mean); + const double sd = std::sqrt(s2 / (NULL_REPLICATES - 1)); + return sd > 0 ? (mean - real) / sd : 0.0; +} + +bool same_rotations_and_centring(const gemmi::GroupOps &a, const gemmi::GroupOps &b) { + if (a.sym_ops.size() != b.sym_ops.size() || a.cen_ops.size() != b.cen_ops.size()) + return false; + for (const gemmi::Op &op : a.sym_ops) + if (b.find_by_rotation(op.rot) == nullptr) + return false; + for (const gemmi::Op::Tran &c : a.cen_ops) + if (std::find(b.cen_ops.begin(), b.cen_ops.end(), c) == b.cen_ops.end()) + return false; + return true; +} + +bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) { + if (!a.is_enantiomorphic()) + return false; + gemmi::GroupOps e = a.operations(); + e.change_basis_forward(a.change_of_hand_op()); + const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e); + return other != nullptr && other->number == b.number; +} + +} // namespace + +// --- the model --- + +SearchModel PrepareSearchModel(const std::string &path) { + SearchModel out; + try { + out.st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect); + } catch (const std::exception &e) { + out.failure_reason = std::string("cannot read the model: ") + e.what(); + return out; + } + if (out.st.models.empty()) { + out.failure_reason = "the model has no atoms"; + return out; + } + out.st.models.resize(1); + gemmi::remove_hydrogens(out.st); + gemmi::Model &model = out.st.models[0]; + for (gemmi::Chain &ch : model.chains) { + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + if (a.b_iso > 100.0f) { + out.failure_reason = "the B-factor column holds values above 100, so it is not a pLDDT: " + "not a predicted model"; + return out; + } + out.residues_read += static_cast(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(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 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(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(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 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(std::lround(PI / step))); + std::vector 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(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 CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell) { + std::vector 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(op.rot[i][j]) / gemmi::Op::DEN; + out.push_back(cell.orth.mat.multiply(f).multiply(cell.frac.mat)); + } + return out; +} + +std::vector RotationAsu(const std::vector &all, + const std::vector &point_group) { + std::vector 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 &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 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 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(std::floor(u[0])), static_cast(std::floor(u[1])), + static_cast(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 v = grid_.data[i + static_cast(grid_.nu) * (j + static_cast(grid_.nv) * l)]; + return conj ? std::conj(std::complex(v)) : std::complex(v); + }; + std::complex 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 &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 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 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(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(std::min(SEARCH_BINS, std::max(n / 100, 1))); + std::vector sum(out.bins, 0.0); + std::vector count(out.bins, 0); + for (size_t i = 0; i < n; i++) { + const int b = static_cast(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 RotationFunction(const MolecularTransform &ft, const SearchReflections &refl, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, + const std::vector &rotations, size_t nthreads) { + const gemmi::GroupOps gops = sg.operations(); + const size_t nops = gops.sym_ops.size(), n = refl.hkl.size(); + const gemmi::Mat33 frac_t = cell.frac.mat.transpose(); + // Copy j of the model contributes |F_model(s)|^2 at s = (h R_j) as a Cartesian reciprocal vector. + std::vector s(n * nops); + for (size_t i = 0; i < n; i++) + for (size_t j = 0; j < nops; j++) + s[i * nops + j] = frac_t.multiply(miller_vec(rotated_index(gops.sym_ops[j], refl.hkl[i]))); + double w_mean = 0; + for (double e2 : refl.e2) + w_mean += e2; + w_mean /= static_cast(n); + double w_ss = 0; + for (double e2 : refl.e2) + w_ss += (e2 - w_mean) * (e2 - w_mean); + + std::vector out(rotations.size(), 0.0); + ParallelFor(static_cast(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 ic(n, 0.0), bin_sum(refl.bins, 0.0); + std::vector bin_n(refl.bins, 0); + for (size_t i = 0; i < n; i++) { + for (size_t j = 0; j < nops; j++) + ic[i] += std::norm(ft.At(rt.multiply(s[i * nops + j]))); + bin_sum[refl.bin[i]] += ic[i]; + bin_n[refl.bin[i]]++; + } + double x_mean = 0; + for (size_t i = 0; i < n; i++) { + ic[i] /= bin_sum[refl.bin[i]] / bin_n[refl.bin[i]]; + x_mean += ic[i]; + } + x_mean /= static_cast(n); + double sxy = 0, sxx = 0; + for (size_t i = 0; i < n; i++) { + sxy += (ic[i] - x_mean) * (refl.e2[i] - w_mean); + sxx += (ic[i] - x_mean) * (ic[i] - x_mean); + } + out[r] = sxx > 0 && w_ss > 0 ? sxy / std::sqrt(sxx * w_ss) : 0.0; + }); + return out; +} + +// --- the translation function --- + +TranslationPeak TranslationFunction(const MolecularTransform &ft, const SearchReflections &refl, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, + const gemmi::Mat33 &rotation) { + // The full-symmetry translation function as one Fourier synthesis, following Crowther & Blow + // (1967) Acta Cryst. 23, 544-548: with the model at t, copy j adds G_j(h) exp(2 pi i (h R_j).t) + // exp(2 pi i h.t_j) to F(h), so |F(h)|^2 summed against E^2 - 1 over the sphere is, for each pair + // of copies j != k, a Fourier series in t with coefficient (E^2 - 1) G_j conj(G_k) exp(2 pi i + // h.(t_j - t_k)) at the index h R_j - h R_k. Every pair's terms go on one grid, and one FFT gives + // the function at every t. The pairs j = k are the rotation function's, constant in t. + const gemmi::GroupOps gops = sg.operations(); + const size_t nops = gops.sym_ops.size(); + TranslationPeak out; + if (nops < 2) + return out; // P1: every position is the same + const gemmi::Mat33 frac_t = cell.frac.mat.transpose(); + const gemmi::Mat33 rt = rotation.transpose(); + // Per copy j: k = h R_j, and g = G_j(h) exp(2 pi i h.t_j), so a pair's coefficient is g_j conj(g_k). + struct Term { gemmi::Miller h; std::vector k; std::vector> g; double w; }; + std::vector terms; + std::vector bin_sum(refl.bins, 0.0); + std::vector bin_n(refl.bins, 0); + std::vector term_bin; + for (size_t i = 0; i < refl.hkl.size(); i++) { + if (gops.is_systematically_absent(refl.hkl[i])) + continue; + Term t{refl.hkl[i], {}, {}, refl.multiplicity[i] * (refl.e2[i] - 1.0)}; + double sum = 0; + for (const gemmi::Op &op : gops.sym_ops) { + t.k.push_back(rotated_index(op, refl.hkl[i])); + const double ht = (t.h[0] * op.tran[0] + t.h[1] * op.tran[1] + t.h[2] * op.tran[2]) + / static_cast(gemmi::Op::DEN); + t.g.push_back(ft.At(rt.multiply(frac_t.multiply(miller_vec(t.k.back())))) * std::polar(1.0, 2.0 * PI * ht)); + sum += std::norm(t.g.back()); + } + bin_sum[refl.bin[i]] += sum; + bin_n[refl.bin[i]]++; + term_bin.push_back(refl.bin[i]); + terms.push_back(std::move(t)); + } + if (terms.empty()) + return out; + int max_m[3] = {0, 0, 0}; + for (const Term &t : terms) + for (size_t j = 0; j < nops; j++) + for (size_t k = 0; k < nops; k++) + for (int a = 0; a < 3; a++) + max_m[a] = std::max(max_m[a], std::abs(t.k[j][a] - t.k[k][a])); + const int nu = good_even_size(2 * max_m[0] + 2), nv = good_even_size(2 * max_m[1] + 2), + nw = good_even_size(2 * max_m[2] + 2); + const int nwh = nw / 2 + 1; + std::vector> c(static_cast(nu) * nv * nwh, 0.0); + for (size_t i = 0; i < terms.size(); i++) { + const Term &t = terms[i]; + const double norm = bin_sum[term_bin[i]] / bin_n[term_bin[i]]; + for (size_t j = 0; j < nops; j++) + for (size_t k = 0; k < nops; k++) { + if (j == k) + continue; + const int m[3] = {t.k[j][0] - t.k[k][0], t.k[j][1] - t.k[k][1], t.k[j][2] - t.k[k][2]}; + // Only l >= 0 is stored: the pair (k, j) carries the conjugate term at -m. + if (m[2] < 0) + continue; + const std::complex v = t.w / norm * t.g[j] * std::conj(t.g[k]); + const int u = (m[0] % nu + nu) % nu, vv = (m[1] % nv + nv) % nv; + c[u + static_cast(nu) * (vv + static_cast(nv) * m[2])] += v; + } + } + // MapFromFPhi sums F exp(-2 pi i m.x): with F = conj(C) that is the conjugate of the series, which + // is real. + gemmi::FPhiGrid hkl; + hkl.unit_cell = cell; + hkl.spacegroup = &gemmi::get_spacegroup_p1(); + hkl.axis_order = gemmi::AxisOrder::XYZ; + hkl.half_l = true; + hkl.set_size_without_checking(nu, nv, nwh); + for (size_t i = 0; i < c.size(); i++) + hkl.data[i] = std::complex(std::conj(c[i])); + const gemmi::Grid 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(map.data[i]) * map.data[i]; + if (map.data[i] > map.data[best]) + best = i; + } + const double count = static_cast(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(map.nu) * map.nv); + out.t = gemmi::Fractional(static_cast(u) / map.nu, static_cast(v) / map.nv, + static_cast(w) / map.nw); + return out; +} + +// --- placements --- + +gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg) { + gemmi::Structure st = model.st; + const gemmi::Position shift = cell.orthogonalize(t); + std::vector pos = ModelPositions(st.models[0]); + for (gemmi::Position &p : pos) + p = gemmi::Position(rotation.multiply(p) + shift); + SetModelPositions(st.models[0], pos); + st.cell = cell; + st.spacegroup_hm = sg.xhm(); + st.setup_cell_images(); + return st; +} + +double ClashFraction(const gemmi::Structure &placed) { + const gemmi::UnitCell &cell = placed.cell; + const gemmi::SpaceGroup *sg = placed.find_spacegroup(); + std::vector ca; + for (const gemmi::Chain &ch : placed.models[0].chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + if (a.name == "CA") + ca.push_back(cell.fractionalize(a.pos)); + if (ca.empty() || sg == nullptr) + return 0.0; + gemmi::Fractional c0; + for (const gemmi::Fractional &f : ca) + c0 = c0 + f; + c0 = gemmi::Fractional(c0 * (1.0 / static_cast(ca.size()))); + double radius = 0; + for (const gemmi::Fractional &f : ca) + radius = std::max(radius, cell.orthogonalize_difference(f - c0).length()); + constexpr double CLASH_A = 3.0; + std::vector clash(ca.size(), 0); + for (const gemmi::Op &op : sg->operations().all_ops_sorted()) { + std::vector mate; + for (const gemmi::Fractional &f : ca) { + const std::array 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(mate.size()))); + const gemmi::Fractional nearest(std::round(c0.x - cm.x), std::round(c0.y - cm.y), std::round(c0.z - cm.z)); + for (int dx = -1; dx <= 1; dx++) + for (int dy = -1; dy <= 1; dy++) + for (int dz = -1; dz <= 1; dz++) { + const gemmi::Fractional shift(nearest.x + dx, nearest.y + dy, nearest.z + dz); + if (op == gemmi::Op::identity() && shift.x == 0 && shift.y == 0 && shift.z == 0) + continue; + if (cell.orthogonalize_difference(cm + shift - c0).length() > 2 * radius + CLASH_A) + continue; + for (size_t a = 0; a < ca.size(); a++) { + if (clash[a]) + continue; + for (const gemmi::Fractional &m : mate) + if (cell.orthogonalize_difference(m + shift - ca[a]).length_sq() < CLASH_A * CLASH_A) { + clash[a] = 1; + break; + } + } + } + } + return static_cast(std::count(clash.begin(), clash.end(), 1)) / static_cast(ca.size()); +} + +std::vector SameRotationGroups(const gemmi::SpaceGroup &sg) { + const gemmi::GroupOps ops = sg.operations(); + std::vector out{sg}; + for (const gemmi::SpaceGroup &g : gemmi::spacegroup_tables::main) { + if (g.number == sg.number && g.xhm() == sg.xhm()) + continue; + const gemmi::GroupOps gops = g.operations(); + if (!same_rotations_and_centring(ops, gops) || gemmi::find_spacegroup_by_ops(gops) != &g) + continue; + out.push_back(g); + } + return out; +} + +// --- the search and the decision --- + +MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &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(data_d_min, r.d); + result.d_min_search = std::max(SEARCH_D_MIN, data_d_min); + result.d_max_search = SEARCH_D_MAX; + result.d_min_score = std::max(SCORE_D_MIN, data_d_min); + result.t_model = seconds_since(t0); + logger.Info("Model prior: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)", + model.residues_kept, model.residues_read, model.mass_da / 1000.0, model.radius_A); + + t0 = std::chrono::steady_clock::now(); + const MolecularTransform ft(model, result.d_min_search, nthreads); + result.t_transform = seconds_since(t0); + const std::vector sphere = UniformRotations(ROTATION_STEP_DEG); + + // One merge: the rotation function over its point group, then for each of its groups the + // translation function and the best placements through packing, the rigid body and R-free. + const auto search_merge = [&](size_t m) { + const MRMerge &merge = merges[m]; + const gemmi::UnitCell cell(merge.cell); + const double d_min = result.d_min_score; + // --- rotation --- + t0 = std::chrono::steady_clock::now(); + const SearchReflections refl = MakeSearchReflections(*merge.merged, cell, merge.space_group, + result.d_min_search, result.d_max_search); + const std::vector pg = CartesianPointGroup(merge.space_group, cell); + const std::vector rotations = RotationAsu(sphere, pg); + const std::vector rf = RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads); + double rf_mean = 0, rf_sd = 0; + for (double x : rf) + rf_mean += x; + rf_mean /= static_cast(rf.size()); + for (double x : rf) + rf_sd += (x - rf_mean) * (x - rf_mean); + rf_sd = std::sqrt(rf_sd / static_cast(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 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 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(peaks.size()) == TRANSLATION_ROTATIONS) + break; + } + result.t_rotation += seconds_since(t0); + logger.Info("Model prior: {} in {}: {} reflections {:.1f}-{:.1f} A, {} orientations, best RF Z {:.1f}", + merge.name, merge.space_group.xhm(), refl.hkl.size(), result.d_max_search, result.d_min_search, + rotations.size(), rf_sd > 0 ? (rf[peaks.front()] - rf_mean) / rf_sd : 0.0); + + // --- translation, for every group and every peak orientation --- + t0 = std::chrono::steady_clock::now(); + const std::vector groups = SameRotationGroups(merge.space_group); + std::vector tf_all(groups.size() * peaks.size()); + ParallelFor(static_cast(tf_all.size()), nthreads, [&](int i) { + tf_all[i] = TranslationFunction(ft, refl, cell, groups[i / peaks.size()], rotations[peaks[i % peaks.size()]]); + }); + result.t_translation += seconds_since(t0); + + for (size_t g = 0; g < groups.size(); g++) { + const gemmi::SpaceGroup &sg = groups[g]; + MRCandidate cand; + cand.merge = static_cast(m); + cand.space_group = sg; + const gemmi::GroupOps gops = sg.operations(); + cand.vm = cell.volume / static_cast(gops.order()) / model.mass_da; + cand.copies_likely = cand.vm >= VM_TWO_COPIES ? 2 : 1; + const TranslationPeak *tf = &tf_all[g * peaks.size()]; + std::vector by_tf(peaks.size()); + std::iota(by_tf.begin(), by_tf.end(), 0); + std::stable_sort(by_tf.begin(), by_tf.end(), [&](size_t a, size_t b) { return tf[a].z > tf[b].z; }); + + // --- the best few through packing, the rigid body and R-free --- + t0 = std::chrono::steady_clock::now(); + RigidBodyEngines engines; + bool have = false; + gemmi::Mat33 orientation; + // Where one copy cannot be the content of the cell - too big for it, or small enough for two - + // the group is not decided on here, and only its best translation-function placement is + // checked for packing. + const bool one_copy = cand.vm >= VM_MIN && cand.copies_likely == 1; + for (int k = 0; k < std::min(one_copy ? RESCORED_PER_GROUP : 1, static_cast(by_tf.size())); k++) { + const size_t i = by_tf[k]; + auto placed = std::make_shared( + PlaceSearchModel(model, rotations[peaks[i]], tf[i].t, cell, sg)); + PlacementScore score{NAN, NAN}; + if (one_copy) { + if (k == 0) + engines = make_engines(*placed, sg, *merge.merged, d_min, NULL_REPLICATES + 1, logger); + score = refine_placement(*placed, sg, *merge.merged, d_min, engines.get(), nthreads); + } + const double clash = ClashFraction(*placed); + if (!have || score.r_free < cand.r_free) { // NaN (not refined) only ever on k = 0 + have = true; + cand.rf_z = rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0; + cand.tf_z = tf[i].z; + cand.clash_fraction = clash; + cand.r_work = score.r_work; + cand.r_free = score.r_free; + cand.placed = placed; + orientation = rotations[peaks[i]]; + } + logger.Info("Model prior: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + sg.xhm(), rf_sd > 0 ? (rf[peaks[i]] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash, + score.r_work, score.r_free); + } + result.t_rescore += seconds_since(t0); + + if (cand.vm < VM_MIN) + cand.verdict = "MATTHEWS"; + else if (cand.copies_likely > 1) + cand.verdict = "MULTI_COPY"; + else if (cand.clash_fraction > CLASH_MAX) + cand.verdict = "CLASH"; + else if (!(cand.tf_z >= TF_Z_PLACED)) + cand.verdict = "NOT_PLACED"; + else + cand.verdict = "PLACED"; + if (cand.verdict == "PLACED") { + t0 = std::chrono::steady_clock::now(); + cand.null_sigma = null_sigma(*cand.placed, orientation, sg, *merge.merged, d_min, engines.get(), + nthreads); + result.t_null += seconds_since(t0); + if (!(cand.null_sigma >= NULL_SIGMA_MIN)) + cand.verdict = "NULL"; + } + result.candidates.push_back(std::move(cand)); + } + }; + + // --- the decision --- + // Within one merge the candidates are scored on the same reflections, so their R-free values can + // be compared directly: the placed candidate with the lowest R-free is chosen where it leads every + // other candidate of that merge that packs by R_FREE_MARGIN. + const auto choose_in = [&](int merge) -> int { + int best = -1; + for (size_t i = 0; i < result.candidates.size(); i++) { + const MRCandidate &c = result.candidates[i]; + if (c.merge == merge && c.verdict == "PLACED" && (best < 0 || c.r_free < result.candidates[best].r_free)) + best = static_cast(i); + } + if (best < 0) + return -1; + for (size_t i = 0; i < result.candidates.size(); i++) { + const MRCandidate &c = result.candidates[i]; + if (static_cast(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX) + continue; + if (!(c.r_free - result.candidates[best].r_free >= R_FREE_MARGIN)) { + result.decision = fmt::format("UNDECIDED: {} and {} within {:.2f} in R-free", + result.candidates[best].space_group.xhm(), c.space_group.xhm(), + R_FREE_MARGIN); + return -2; + } + } + return best; + }; + // A merge whose every group the model cannot be in at all - one copy does not fit the cell, or + // every placement clashes with its own symmetry mates. + const auto cannot_hold = [&](int merge) { + bool any = false; + for (const MRCandidate &c : result.candidates) + if (c.merge == merge) { + any = true; + if (c.verdict != "MATTHEWS" && c.verdict != "CLASH") + return false; + } + return any; + }; + // The twin hypothesis' merge is searched only where it can decide something: where the answer's + // merge cannot hold the model in any of its groups. + search_merge(0); + const int answer = choose_in(0); + if (answer >= 0) { + result.chosen = answer; + } else if (answer == -1 && merges.size() > 1 && cannot_hold(0)) { + search_merge(1); + const int twin = choose_in(1); + if (twin >= 0) + result.chosen = twin; + else if (twin == -1) + result.decision = "NOT_PLACED: no group of either merge holds a placed model"; + } else if (answer == -1) { + result.decision = "NOT_PLACED: no group of the answer's merge holds a placed model"; + } + if (result.chosen >= 0) { + const MRCandidate &c = result.candidates[result.chosen]; + const gemmi::SpaceGroup &data = merges[0].space_group; + if (c.merge == 0 && c.space_group.xhm() == data.xhm()) + result.decision = "AGREES"; + else if (c.merge == 0 && enantiomorphs(data, c.space_group)) + result.decision = "ENANTIOMORPH"; + else if (c.merge == 0) + result.decision = "SPACE_GROUP"; + else + result.decision = "TWIN_SUBGROUP"; + } + result.ok = true; + logger.Info("Model prior: {} (model {:.1f} s, transform {:.2f} s, rotation {:.2f} s, translation {:.2f} s, " + "rescoring {:.2f} s, null {:.2f} s)", + result.chosen >= 0 ? result.decision + " - " + result.candidates[result.chosen].space_group.xhm() + : result.decision, + result.t_model, result.t_transform, result.t_rotation, result.t_translation, result.t_rescore, + result.t_null); + return result; +} diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h new file mode 100644 index 000000000..89323630a --- /dev/null +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -0,0 +1,165 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +// Molecular replacement of ONE copy of a predicted model (rugnux --af-model, experimental): a global +// search for the model's orientation and position in the cell, so that a model with no placement of +// its own - an AlphaFold prediction - can be scored against the data, and the space groups the data +// cannot tell apart can be scored against it. +// +// - The model: residues with pLDDT < 70 removed, B-factors from pLDDT, hydrogens dropped. +// - Rotation: a uniform sampling of the point group's rotation ASU, each orientation scored by the +// correlation of the observed E^2 with the model's own symmetry-summed |F|^2 - read off the +// model's molecular transform, computed once by FFT and interpolated per orientation. +// - Translation: the full-symmetry translation function on one FFT per orientation and space group. +// - Rescoring: the best placements through the rigid body, bulk solvent, R-free, a packing check +// against the symmetry mates and a null of the same model in random orientations. +// +// Everything up to the rescoring runs on the CPU, in both builds; the rigid body runs on the GPU +// where there is one. + +#include +#include +#include +#include + +#include "gemmi/model.hpp" +#include "gemmi/recgrid.hpp" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +#include "../../common/Reflection.h" // MergedReflection +#include "../../common/UnitCell.h" + +class Logger; + +// The predicted model as the search uses it: trimmed, B-factors from pLDDT, centroid at the origin. +struct SearchModel { + gemmi::Structure st; + int residues_read = 0; + int residues_kept = 0; + double mass_da = 0; // of the residues kept, 110 Da each + double radius_A = 0; // the farthest atom from the centroid + std::string failure_reason; +}; +SearchModel PrepareSearchModel(const std::string &path); + +// Rotations covering SO(3) with neighbours about `step_deg` apart: the identity, then shells of +// rotation angle kappa = step, 2 step, ... 180 deg, each with its axes on a golden-spiral lattice of +// the sphere as dense as 2 sin(kappa/2) - the distance a change of axis moves a rotation - needs. +std::vector UniformRotations(double step_deg); + +// The rotations of the space group's point group, in Cartesian coordinates of `cell`. +std::vector CartesianPointGroup(const gemmi::SpaceGroup &sg, const gemmi::UnitCell &cell); + +// The rotations of `all` in the rotation ASU of `point_group`: the crystal puts a copy of the model +// in every orientation S R, so only one rotation of each class {S R} is worth scoring - here the one +// closest to the identity (largest trace), the classes' Dirichlet region about it. +std::vector RotationAsu(const std::vector &all, + const std::vector &point_group); + +// The angle in degrees between a and b, the nearest over the classes {S b}. +double RotationDistanceDeg(const gemmi::Mat33 &a, const gemmi::Mat33 &b, + const std::vector &point_group); + +// The structure factors of one copy of the model, at its centroid, at any reciprocal vector. +class MolecularTransform { +public: + MolecularTransform(const SearchModel &model, double d_min, size_t nthreads); + // F at the Cartesian reciprocal vector s (1/A), interpolated trilinearly. + std::complex At(const gemmi::Vec3 &s) const; + +private: + gemmi::FPhiGrid grid_; // half-l, on a P1 box `box_` A wide + double box_ = 0; + double blur_ = 0; +}; + +// The observed reflections the search scores against: E^2 in resolution bins, with each unique +// reflection's multiplicity on the whole sphere. +struct SearchReflections { + std::vector hkl; + std::vector e2; + std::vector multiplicity; + std::vector bin; + int bins = 0; +}; +SearchReflections MakeSearchReflections(const std::vector &merged, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg, double d_min, double d_max); + +// The rotation function at each of `rotations`: the correlation of E^2 - 1 with the bin-normalised +// sum over the point group's copies of |F_model|^2. +std::vector RotationFunction(const MolecularTransform &ft, const SearchReflections &refl, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, + const std::vector &rotations, size_t nthreads); + +// The translation function of the model in orientation `rotation` in `sg`: its highest point as a +// Z-score over every point of the function, and where it is (the fractional position of the +// model's centroid). +struct TranslationPeak { + double z = 0; + gemmi::Fractional t; +}; +TranslationPeak TranslationFunction(const MolecularTransform &ft, const SearchReflections &refl, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, + const gemmi::Mat33 &rotation); + +// The model placed: rotated about its centroid by `rotation`, the centroid at `t`, in `cell` and `sg`. +gemmi::Structure PlaceSearchModel(const SearchModel &model, const gemmi::Mat33 &rotation, const gemmi::Fractional &t, + const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg); + +// The fraction of the model's C-alpha atoms within 3 A of a C-alpha of a symmetry mate. +double ClashFraction(const gemmi::Structure &placed); + +// The space groups with the rotations and the centring of `sg`, in its setting: those merged +// intensities cannot tell apart, apart from systematic absences. `sg` itself first. +std::vector SameRotationGroups(const gemmi::SpaceGroup &sg); + +// One merge the model is searched against: the answer's, or the twin hypothesis' in a subgroup. +struct MRMerge { + std::string name; // "MERGE" or "TWIN_SUBGROUP_MERGE" + const std::vector *merged = nullptr; + UnitCell cell; + gemmi::SpaceGroup space_group; // the group it was merged in +}; + +// One space group's best single-copy placement. +struct MRCandidate { + int merge = 0; // index into the merges searched + gemmi::SpaceGroup space_group; + double vm = 0; // Matthews coefficient of one copy + int copies_likely = 1; + double rf_z = 0; // of the orientation, over the rotation function + double tf_z = 0; + double clash_fraction = 0; + double r_work = NAN, r_free = NAN; // after the rigid body, to d_min_score + double null_sigma = NAN; // R-work above its random-orientation null, in sd + std::string verdict; // PLACED, or why not + std::shared_ptr placed; +}; + +struct MolecularReplacementResult { + bool ok = false; + std::string failure_reason; + std::string model_path; + int residues_read = 0, residues_kept = 0; + double mass_da = 0; + double d_min_search = 0, d_max_search = 0; + double d_min_score = 0; // the limit R-work / R-free are computed to + std::vector candidates; + // The decision: the candidate the model chooses, or -1 where it chooses none, and in words why. + int chosen = -1; + std::string decision; + // Seconds per step. + double t_model = 0, t_transform = 0, t_rotation = 0, t_translation = 0, t_rescore = 0, t_null = 0; +}; + +// Search `merges` (the answer's first) for one copy of the predicted model at `model_path`, in every +// space group with the rotations of each merge's group, and decide. A candidate is PLACED where its +// translation function stands out (TF Z >= 8), it packs, one copy fits the cell, and its R-work beats +// its own null. The model chooses among the answer merge's groups where one placed candidate leads +// every other by an R-free margin; it chooses the twin hypothesis' merge only where none of the +// answer's groups can hold the model at all (Matthews, packing) and a group of the subgroup can. +MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, + Logger &logger, size_t nthreads); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 6c46f28eb..f1f8f32ad 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2562,6 +2562,81 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, doc.sections.push_back(std::move(s)); } + // ------------------------------------------------- MODEL PRIOR (--af-model, experimental) + if (result.model_prior.has_value()) { + const auto &mp = *result.model_prior; + ReportSection s; + s.title = "MODEL PRIOR (EXPERIMENTAL)"; + Add(s, KeyText("AF_MODEL_FILE", mp.model_path)); + if (!mp.ok) { + Add(s, KeyEnum("MODEL_PRIOR", "NOT_PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); + Add(s, KeyText("MODEL_PRIOR_REASON", mp.failure_reason)); + } else { + Add(s, KeyEnum("MODEL_PRIOR", "PERFORMED", {"PERFORMED", "NOT_PERFORMED"})); + Add(s, KeyInt("AF_RESIDUES_READ", mp.residues_read)); + Add(s, KeyInt("AF_RESIDUES_KEPT", mp.residues_kept)); + Add(s, KeyReal("SEARCH_D_MIN", mp.d_min_search, "{:.2f}")); + Add(s, KeyReal("SCORE_D_MIN", mp.d_min_score, "{:.2f}")); + Add(s, KeyText("MODEL_PRIOR_DECISION", mp.decision)); + if (mp.chosen >= 0) + Add(s, KeyText("MODEL_PRIOR_SPACE_GROUP", mp.candidates[mp.chosen].space_group.xhm())); + if (!result.model_prior_file_prefix.empty()) + Add(s, KeyText("MODEL_PRIOR_FILE_PREFIX", result.model_prior_file_prefix)); + Add(s, KeyReal("MODEL_PRIOR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation + + mp.t_rescore + mp.t_null, "{:.2f}")); + Add(s, Blank()); + ReportEntry e; + e.kind = ReportEntry::Kind::Table; + e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "RF_Z", "TF_Z", "CLASH", "R_WORK", "R_FREE", + "NULL_SIGMA", "VERDICT"}; + e.table.text_header = " MERGE SPACE_GROUP VM RF_Z TF_Z CLASH R_WORK R_FREE NULL_SIGMA VERDICT\n" + " ------------------- ----------- ----- ----- ----- ----- ------ ------ ---------- ----------"; + for (const auto &c : mp.candidates) { + const std::string merge = c.merge == 0 ? "MERGE" : "TWIN_SUBGROUP_MERGE"; + const auto r = [](double v) { return std::isfinite(v) ? fmt::format("{:.3f}", v) : std::string("-"); }; + e.table.text_rows.push_back(fmt::format( + " {:<19s} {:<11s} {:5.2f} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}", + merge, c.space_group.xhm(), c.vm, c.rf_z, c.tf_z, 100 * c.clash_fraction, r(c.r_work), + r(c.r_free), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), + c.verdict)); + std::vector row; + row.push_back(KeyText("", merge).value); + row.push_back(KeyText("", c.space_group.xhm()).value); + row.push_back(KeyReal("", c.vm, "{:.2f}").value); + row.push_back(KeyReal("", c.rf_z, "{:.1f}").value); + row.push_back(KeyReal("", c.tf_z, "{:.1f}").value); + row.push_back(KeyReal("", c.clash_fraction, "{:.3f}").value); + row.push_back(KeyReal("", c.r_work, "{:.3f}").value); + row.push_back(KeyReal("", c.r_free, "{:.3f}").value); + row.push_back(KeyReal("", c.null_sigma, "{:+.1f}").value); + row.push_back(KeyText("", c.verdict).value); + e.table.cells.push_back(std::move(row)); + } + Add(s, std::move(e)); + Add(s, Blank()); + Add(s, Prose(fmt::format( + " One copy of the predicted model, trimmed to pLDDT >= 70, was searched for in every space group\n" + " each merge cannot tell apart (rotation function on a 6 deg sampling, translation function,\n" + " {:.1f}-{:.1f} A), and the best placements taken through the rigid body; R_WORK and R_FREE are\n" + " to SCORE_D_MIN. PLACED needs TF_Z >= 8,\n" + " at most 5% of the C-alpha atoms clashing with symmetry mates, one copy fitting the cell\n" + " (MATTHEWS: Vm below 1.5; MULTI_COPY: Vm of 4 or more, not supported yet) and an R-work at\n" + " least 3 sd better than the same model refitted from random orientations. The model chooses\n" + " within a merge where its placed group leads every other that packs by 0.03 in R-free, and\n" + " the twin hypothesis' merge - searched only then - where no group of the answer's merge can\n" + " hold it at all. A group one copy cannot fill is checked for packing only, not refined.\n" + " The answer's files take the enantiomorph from it and nothing else; another group it\n" + " chooses is written as {}.mtz/.cif with the placed model.", + mp.d_max_search, mp.d_min_search, + result.model_prior_file_prefix.empty() ? std::string("_mr") : result.model_prior_file_prefix))); + Add(s, Prose(fmt::format( + " Seconds: model {:.2f}, transform {:.2f}, rotation {:.2f}, translation {:.2f}, rigid body and\n" + " R-free {:.2f}, null {:.2f}.", mp.t_model, mp.t_transform, mp.t_rotation, mp.t_translation, + mp.t_rescore, mp.t_null), true)); + } + doc.sections.push_back(std::move(s)); + } + if (result.cancelled) Warn(doc, PathologyCode::CANCELLED, fmt::format("Processing was cancelled after {} images - this report describes an " diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index a476554ac..09179e607 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -32,6 +32,7 @@ #include "../image_analysis/scale_merge/TranslationalNCS.h" // TranslationalNCSResult #include "../image_analysis/indexing/HarmonicContamination.h" // HarmonicContaminationResult #include "../image_analysis/structure_refinement/ModelValidation.h" // ModelValidationResult +#include "../image_analysis/structure_refinement/MolecularReplacement.h" // MolecularReplacementResult #include "../image_analysis/indexing/LeftoverLattices.h" // LeftoverLatticeCensus #include "ReportDocument.h" // ReportWarning, PathologyCode #include "../image_analysis/scale_merge/SearchSpaceGroup.h" // SearchSpaceGroupResult @@ -196,6 +197,12 @@ struct ProcessConfig { // Fo-Fc maps (--model). std::string model_path; + // Experimental (--af-model): a predicted model (AlphaFold mmCIF, pLDDT in the B column) placed by + // molecular replacement, one copy, in every space group the merge cannot tell apart and in the twin + // hypothesis' subgroup. It never changes the answer's files beyond the enantiomorph label; a + // different group it chooses is written beside them as _mr. + std::string af_model_path; + // --finalist-ledger: report-only symmetry evidence table (and the twin-immune zones), folded from // a merge of its own in P1 that only this flag makes. bool finalist_ledger = false; @@ -562,6 +569,11 @@ struct ProcessResult { // that was asked to validate looking like one that was not. Unset when no model was asked for. std::optional model_validation; + // What --af-model found, when one was given; failure_reason says why where it could not run. + std::optional model_prior; + // The file prefix the model's choice was written under (_mr), empty where nothing was. + std::string model_prior_file_prefix; + // Space group used for the final re-scale/merge and written to the master file: determined by the // search when the user did not fix one, otherwise the fixed group. The whole group, so it keeps // its setting - a number only ever names the reference one. diff --git a/rugnux/RugnuxCommandLine.cpp b/rugnux/RugnuxCommandLine.cpp index c05d31737..3bb2a2fe9 100644 --- a/rugnux/RugnuxCommandLine.cpp +++ b/rugnux/RugnuxCommandLine.cpp @@ -82,6 +82,8 @@ std::string RugnuxCommandLine(const ProcessConfig &config, } if (!config.model_path.empty()) add("--model", config.model_path); + if (!config.af_model_path.empty()) + add("--af-model", config.af_model_path); } if (config.start_image != 0) add("-s", std::to_string(config.start_image)); diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 90b14acf2..4909b2247 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3968,6 +3968,63 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } prepare_twin_file(); + // --af-model (experimental): the predicted model placed by molecular replacement, one copy, in the + // groups the answer's merge cannot tell apart and in the twin hypothesis' subgroup. It takes the + // enantiomorph label into the answer, as --model does, and nothing else: a different group it + // chooses goes to files of its own beside the answer's, _mr, from the merge it was chosen on. + struct ModelPriorFile { + std::vector merged; + UnitCell cell; + gemmi::SpaceGroup space_group; + std::optional law; + TwinLawFraction fraction; + std::shared_ptr placed; + }; + std::optional model_prior_file; + if (result.consensus_cell && write_files && !config_.af_model_path.empty() && !geometry_prepass + && !superseded) { + phase("Placing the predicted model"); + Step("Model prior"); + TimingMark mark(logger, "model prior"); + std::vector merges{{"MERGE", &sm.merged, *result.consensus_cell, experiment_.GetSpaceGroupOrP1()}}; + if (twin_file) + merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group}); + MolecularReplacementResult mr = RunMolecularReplacement(config_.af_model_path, merges, logger, + static_cast(config_.nthreads)); + if (!mr.ok) + result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED, + "The predicted model was not placed: " + mr.failure_reason}); + if (mr.chosen >= 0) { + const MRCandidate &c = mr.candidates[mr.chosen]; + if (mr.decision == "ENANTIOMORPH") { + // A change of label only (see ValidateAgainstModel): the two groups index identically. + experiment_.SetSpaceGroup(c.space_group); + result.space_group = c.space_group; + if (twin_file && twin_file->space_group.is_enantiomorphic()) { + gemmi::GroupOps e = twin_file->space_group.operations(); + e.change_basis_forward(twin_file->space_group.change_of_hand_op()); + if (const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e)) + twin_file->space_group = *other; + } + logger.Info("Model prior: the answer is relabelled {}, the enantiomorph the placed model chose", + c.space_group.xhm()); + } else if (mr.decision == "SPACE_GROUP") { + ModelPriorFile f{sm.merged, *result.consensus_cell, c.space_group, std::nullopt, {}, c.placed}; + const gemmi::GroupOps gops = c.space_group.operations(); + std::erase_if(f.merged, [&](const MergedReflection &r) { + return gops.is_systematically_absent({{r.h, r.k, r.l}}); + }); + model_prior_file = std::move(f); + } else if (mr.decision == "TWIN_SUBGROUP" && twin_file) { + model_prior_file = ModelPriorFile{twin_file->merged, twin_file->cell, c.space_group, twin_file->law, + twin_file->fraction, c.placed}; + } + } + if (model_prior_file) + result.model_prior_file_prefix = config_.output_prefix + "_mr"; + result.model_prior = std::move(mr); + } + // The unmerged MTZ (written with the other output files, see there) is built from here on, // beside the reflection files: the indexing and the group it is written in are final once the // model has been read, and nothing after this changes the integrated observations it reads. It @@ -4032,6 +4089,33 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { WriteMmcifReflections(twin_file->merged, twin_file->cell, twin_x, twin_merge->statistics, twin_em, law, config_.provenance, prefix + ".cif", static_cast(std::max(1, config_.nthreads))); } + // The predicted model's choice (see model_prior_file): its reflections, from the merge it was + // chosen on and labelled with the group it chose, and the model as it placed it. + if (model_prior_file) { + const ModelPriorFile &f = *model_prior_file; + DiffractionExperiment mr_x = experiment_; + mr_x.SetSpaceGroup(f.space_group); + TwinningAnalysisResult law; + const MergeStatistics *statistics = &sm.statistics; + ErrorModelReport mr_em = em_report; + if (f.law) { + law.twin_law_operator = f.law->as_hkl().triplet(); + law.twin_law_fraction = std::isfinite(f.fraction.fraction) ? f.fraction.fraction : 0.5; + statistics = &twin_merge->statistics; + mr_em = ErrorModelReport{ + twin_merge->isa > 0 && twin_merge->isa_resolved ? fmt::format("{:.2f}", twin_merge->isa) : "?", + std::string(), twin_merge->error_model_a > 0 ? fmt::format("{:.3f}", twin_merge->error_model_a) + : std::string(), + twin_merge->error_model_b > 0 ? fmt::format("{:.4e}", twin_merge->error_model_b) + : std::string()}; + } + const std::string prefix = result.model_prior_file_prefix; + WriteMtzReflections(f.merged, f.cell, mr_x, config_.provenance, prefix + ".mtz"); + WriteMmcifReflections(f.merged, f.cell, mr_x, *statistics, mr_em, law, config_.provenance, prefix + ".cif", + static_cast(std::max(1, config_.nthreads))); + if (f.placed) + WritePlacedModel(*f.placed, f.cell, f.space_group, prefix, logger); + } writing_mark.Stop(); // FINALIST LEDGER (--finalist-ledger), report-only and deliberately behind a flag: these diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index 10f78955c..32057f5bd 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -202,6 +202,7 @@ void print_usage() { std::cout << " -z, --reference Reference data of the same crystal form: an MTZ, or a PDB structure-factor mmCIF (e.g. 1abc-sf.cif, first merged block with a usable column; R-free set from _refln.status f/o, else pdbx_r_free_flag). Gzipped or not, recognised by content. --reference-mtz is the same option" << std::endl; std::cout << " --reference-column