From 0d33d7a40a6af1c9c9894025978f058aea55d6dd Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 11:41:30 +0200 Subject: [PATCH] rugnux --af-model: likelihood judging, PAE domains, twin comparison, tNCS carry test, faster multi-copy search - Likelihood (Read 2001; Z-score over random positions as McCoy et al. 2005): log-likelihood gain of Rice/Woolfson E distributions, sigma_A from the pLDDT coordinate error with its scale refined, to 2.5 A (or the 20000 lowest-resolution reflections). Further-copy TF peaks ranked by the gain on the search's reflections (peaks that lower it dropped); a copy whose TF Z is below 8 also stands out where its refined likelihood Z-score reaches 8 - only where the copies placed are as many as Matthews makes most likely. The likelihood is computed only where a translation function does not stand out. Copy structure factors exact from one FFT of the cell's grid, or read off the molecular transform where it reaches the likelihood's resolution. - PAE domains: the AlphaFold DB's predicted_aligned_error JSON beside the model (same name pattern) is read automatically; average linkage on the symmetrised PAE (10 A, domains >= 40 residues). Where the whole model is placed nowhere, each domain is searched by itself side by side; placed domains must agree, and a domain alone may only confirm the answer's group. - Twin comparison: where the answer's merge chose and the twin hypothesis' merge exists, the chosen placement read as a crystal in the subgroup (its symmetry images as copies, the subgroup's setting found from the two cells) against the model searched for afresh in the subgroup, both against the data twinned by the law (fraction refined on the working set, Yeates 1997); the subgroup is taken by a 0.03 R-free lead plus its own null/margin -> TWIN_SUBGROUP, written to p_mr with law and fraction. With --af-model rugnux makes the subgroup merge wherever a metric re-ask promoted the indexed group. - tNCS rule: a group with fewer copies beside tNCS (screws bridged by the vector, or the enantiomorph) no longer leaves the choice undecided where the chosen copies, carried into it and refined, fit worse by 0.03 in R-free. - Speed: the molecular transform's FFT on the GPU (cuFFT, serialised), no rigid body for trials that lower the likelihood, re-search only with tNCS, PAE domains searched in parallel. Survey (193 sets, MR replay on stored merges): single GREEN 58 -> 68, multi 48 -> 58, partial 8 -> 9; 4bwl (and, in full runs, 2wnq and a withdrawn entry of the same protein) TWIN_SUBGROUP P21 instead of a wrong AGREES; a P21212 set with tNCS now SPACE_GROUP P21212 instead of UNDECIDED; no wrong call. CPU and GPU builds make the same decisions on 18 demo/multi-copy/twin sets. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/ACKNOWLEDGEMENT.md | 16 +- docs/RUGNUX_REPORT.md | 53 +- .../MolecularReplacement.cpp | 1019 +++++++++++++++-- .../MolecularReplacement.h | 117 +- .../MolecularReplacementGPU.cu | 65 ++ .../MolecularReplacementGPU.h | 7 + rugnux/ResultReport.cpp | 43 +- rugnux/RugnuxScaleMerge.cpp | 15 +- rugnux/rugnux_cli.cpp | 2 +- tests/MolecularReplacementTest.cpp | 142 ++- 10 files changed, 1341 insertions(+), 138 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index b42dccda8..d6264b1b6 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -473,7 +473,21 @@ complex crystals" (2003), Protein Sci. 12, 1865-1871 [doi:10.1110/ps.0350503](https://doi.org/10.1110/ps.0350503). A copy refined beside the others has the others' placements read off by least-squares superposition: W. Kabsch, "A solution for the best rotation to relate two sets of vectors" (1976), Acta Cryst. A32, 922-923 -[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873). +[doi:10.1107/S0567739476001873](https://doi.org/10.1107/S0567739476001873). The placements are judged by +their log-likelihood gain, from the Rice and Woolfson distributions of the observed normalised amplitudes +given the model's, with sigma_A from the model's expected coordinate error: R. J. Read, "Pushing the +boundaries of molecular replacement with maximum likelihood" (2001), Acta Cryst. D57, 1373-1382 +[doi:10.1107/S0907444901012471](https://doi.org/10.1107/S0907444901012471); and the gain is read as a +Z-score over the same orientation at other positions, as the likelihood-enhanced translation function's +is: A. J. McCoy, R. W. Grosse-Kunstleve, L. C. Storoni and R. J. Read, "Likelihood-enhanced fast +translation functions" (2005), Acta Cryst. D61, 458-464 +[doi:10.1107/S0907444905001617](https://doi.org/10.1107/S0907444905001617). The logarithm of the Bessel +function I0 in it is evaluated with the polynomial approximations of M. Abramowitz and I. A. Stegun, +"Handbook of Mathematical Functions" (1964), National Bureau of Standards, 9.8.1-9.8.2. A predicted +model is cut into rigid domains by its predicted aligned error, as Oeffner et al. (above) do before +molecular replacement; the error itself is AlphaFold's: 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). ## Software and computing methods diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index eff93826a..923b68e6b 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -980,14 +980,30 @@ merge cannot tell apart, and the best three placements per group that pack throu solvent and R-free (to `SCORE_D_MIN`, 3 A or the data's limit). Further copies follow one at a time, up to the largest copy number with a Matthews probability of 1% (`COPIES_LIKELY`: the most probable): the translation function again with the copies already placed as a known partial structure, its peaks -ranked by the correlation of all the copies with the data, each new copy refined beside the others held -fixed, and kept where its `TF_Z` reaches 8, it packs and R-free drops by 0.01 (`COPIES`, and -`MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched side by side; every group the best -one leads is searched once more with the best one's copies' orientations tried first, and the better -search kept. +ranked by the likelihood gain of the new copy beside them (peaks that lower the likelihood dropped), +each new copy refined beside the others held fixed, and kept where it stands out (below), it packs and +R-free drops by 0.01 (`COPIES`, and `MODEL_PRIOR_COPIES=` for the group chosen). The groups are searched +side by side; with tNCS, every group that placed fewer copies than the best one is searched once more +with the best one's copies' orientations tried first, and the better search kept. -The table has one row per group and merge. A row is `PLACED` where one of its copies' translation- -function Z-scores (`TF_Z`, the highest) reaches 8, at most 5% of the C-alpha atoms clash with another +The likelihood is the log-likelihood gain of the Rice/Woolfson distributions of E_obs given E_calc, with +sigma_A falling off with the coordinate error the pLDDTs predict and its scale refined against the data, +summed to 2.5 A (or the data's limit, or the 20000 lowest-resolution reflections in a large cell). A copy +stands out where its translation-function Z-score reaches 8, or - only where the copies placed are as +many as the cell most probably holds - where its likelihood after the rigid body reaches a Z-score of 8 +over the same copy in the same orientation at random positions (`LLG_Z`; taken only where `TF_Z` is +below 8). That second way in is what places a single copy in a small cell, whose few reflections keep +the translation function's own Z-score low. + +Where the AlphaFold DB's predicted-aligned-error file sits beside the model under its own name +(`AF--F1-model_v.cif` and `AF--F1-predicted_aligned_error_v.json`), the prediction is cut +into rigid domains by average linkage on the PAE (domains a mean 10 A apart, of 40 residues or more; +`AF_PAE_DOMAINS=`). 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). `MODEL_PRIOR_SEARCH_MODEL=` names what the decision +was made with. + +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 null - the same model refitted from random orientations, each placed by the translation function as the search placed the real one (all copies turned together about their centroid where there are several) - @@ -995,9 +1011,12 @@ is built for the candidate the choice falls on; its R-work has to beat the null' (`NULL_SIGMA`: 3 sd allowing for the null being nine draws, as `MODEL_FIT_SIGMA`), and one that does not is marked `NULL` and the choice falls on the next. A placement clashing over twice the limit is not refined. Translational NCS is looked for in the native Patterson (a peak off the lattice of a fifth of -the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's and -whose search placed fewer copies leaves the choice `UNDECIDED`, since one copy beside a translated twin -of itself fits the data no better than at random. The rotation and translation functions, the rigid +the origin's): with it, a group whose screw axes the tNCS translation turns into the chosen group's (or +the chosen group's enantiomorph) and whose search placed fewer copies leaves the choice `UNDECIDED`, +since one copy beside a translated twin of itself fits the data no better than at random - unless the +chosen copies, carried into that group as they are and refined there, fit worse by 0.03 in R-free: +where the translation does make up the difference between the groups, the copies build the same crystal +in both. 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 (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` is then: @@ -1012,10 +1031,16 @@ GPU where there is one (`MODEL_PRIOR_SEARCH_DEVICE=`). `MODEL_PRIOR_DECISION=` i `ENANTIOMORPH_VETOED` and the label stays; - `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; +- `TWIN_SUBGROUP` - the twin hypothesis' subgroup: either no group of the answer's merge can hold the + model at all (`MATTHEWS`/`CLASH` in every row) and a group of the subgroup can, or the twin comparison + chose it. The comparison runs where the answer's merge chose a group and the twin hypothesis' merge + is there: on that merge, against its data as twinned by its law (the fraction refined on the working + set), the chosen placement with all its symmetry images - the group's crystal read in the subgroup - + is scored against the model searched for afresh in the subgroup; the subgroup is taken where its own + placement leads by 0.03 in R-free and passes the subgroup's own choice (its null, its lead over the + subgroup's other groups). Written as `_mr.mtz`/`.cif` from the twin hypothesis' merge, with + its law and fraction. With `--af-model`, the twin hypothesis' merge is made wherever a metric re-ask + promoted the indexed lattice's group, so that the comparison can run; - `UNDECIDED` or `NOT_PLACED` - nothing is written. With `SPACE_GROUP` and `TWIN_SUBGROUP` the placed model is written as `_mr_model.pdb`/`.cif`. diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 6d1be24f6..23e1b4600 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -17,6 +17,9 @@ #include #include #include +#include + +#include #include // read_structure_gz #include // remove_hydrogens @@ -62,6 +65,10 @@ constexpr int PEAK_REFINE_ITERATIONS = 4; // 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 likelihood a placement is judged by is summed over the reflections to LLG_D_MIN - finer than the +// search, where the information a small cell's few reflections hold is - or over as many as LLG_TERMS. +constexpr double LLG_D_MIN = 2.5; +constexpr size_t LLG_TERMS = 20000; // 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. @@ -147,6 +154,9 @@ int good_even_size(int n) { struct PlacementScore { double r_work = 1, r_free = 1; double k_sol = NAN, b_sol = NAN; + // Against the data as twinned by a law (score_placement's `law`): at the fraction that gives the lowest + // R-work. + double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN; }; gemmi::AsuData> working_amplitudes(const std::vector &merged, @@ -272,9 +282,61 @@ class ScoreEngines {}; // The model's Fcalc and bulk-solvent mask to d_min - on the GPU where `engines` has one - scaled to the // working amplitudes as model validation scales them (ModelValidation.cpp, fit_model), and the R-factors // over the reflections to d_min. `solvent`: k_sol and b_sol held at its, not fitted. +// R-work and R-free against the data as twinned by `law`: |F_twin(h)|^2 = (1-a)|F(h)|^2 + a|F(law h)|^2 of the +// scaled Fcalc, one overall scale refitted at each fraction a, a scanned from 0 to 0.5 in steps of 0.01 and +// chosen on the working set. +// Following Yeates (1997) Methods Enzymol. 276, 344-358 +void twin_r_factors(const gemmi::AsuData> &fmodel, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, const gemmi::Op &law, + PlacementScore &out) { + const gemmi::GroupOps gops = sg.operations(); + const gemmi::ReciprocalAsu asu(&sg); + std::unordered_map model_i; + for (const auto &hv : fmodel.v) + model_i[hkl_key(hv.hkl)] = std::norm(hv.value); + struct Term { double fo, i1, i2; bool free; }; + std::vector terms; + 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; + const auto a = model_i.find(hkl_key(h)); + const auto b = model_i.find(hkl_key(asu.to_asu(law.apply_to_hkl(h), gops).first)); + if (a != model_i.end() && b != model_i.end()) + terms.push_back({r.F, a->second, b->second, r.rfree_flag}); + } + double best = INFINITY; + for (int step = 0; step <= 50; step++) { + const double a = 0.01 * step; + double fo_fc = 0, fc_fc = 0; + for (const Term &t : terms) + if (!t.free) { + const double fc = std::sqrt((1 - a) * t.i1 + a * t.i2); + fo_fc += t.fo * fc; + fc_fc += fc * fc; + } + const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0; + double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + for (const Term &t : terms) { + const double d = std::fabs(t.fo - k * std::sqrt((1 - a) * t.i1 + a * t.i2)); + (t.free ? num_f : num_w) += d; + (t.free ? den_f : den_w) += t.fo; + } + if (den_w > 0 && den_f > 0 && num_w / den_w < best) { + best = num_w / den_w; + out.r_work_twin = best; + out.r_free_twin = num_f / den_f; + out.twin_fraction = a; + } + } +} + PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGroup &sg, const std::vector &merged, double d_min, size_t nthreads, - ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr) { + ScoreEngines *engines = nullptr, const PlacementScore *solvent = nullptr, + const gemmi::Op *law = nullptr) { gemmi::AsuData> fmask, fmodel; #ifdef JFJOCH_USE_CUDA if (ModelStructureFactorsGPU *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr) @@ -299,7 +361,10 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads); } scaling.scale_data(fmodel, &fmask); - return r_factors(fmodel, scaling, sg, merged, d_min); + PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min); + if (law != nullptr) + twin_r_factors(fmodel, sg, merged, d_min, *law, out); + return out; } // The rigid-body engines for one model in one group, where there is a GPU. @@ -359,6 +424,28 @@ PlacementScore refine_placement(gemmi::Structure &st, const gemmi::SpaceGroup &s return score_placement(st, sg, merged, d_min, nthreads, engines); } +// The copies of `assembly` (its chains split evenly, in order) as a crystal of group `sg`, each refined beside +// the others, and scored. +PlacementScore carry_copies(const gemmi::Structure &assembly, int copies, const gemmi::SpaceGroup &sg, + const std::vector &merged, double d_min, size_t nthreads, + ScoreEngines *scorers, Logger &logger) { + std::vector body = SplitCopies(assembly, copies); + for (gemmi::Structure &b : body) { + b.spacegroup_hm = sg.xhm(); + b.setup_cell_images(); + } + RigidBodyEngines engines = make_engines(body[0], sg, merged, d_min, 1, logger); + for (int k = 0; k < copies; k++) { + std::vector rest; + for (int o = 0; o < copies; o++) + if (o != k) + rest.push_back(body[o]); + const gemmi::Structure fixed = AssembleCopies(rest); + refine_rigid_body(body[k], sg, merged, d_min, engines.get(), nthreads, nullptr, &fixed.models[0]); + } + return score_placement(AssembleCopies(body), sg, merged, d_min, nthreads, scorers); +} + // The fractional position of the model's centroid. gemmi::Fractional centroid(const gemmi::Structure &st) { const std::vector pos = ModelPositions(st.models[0]); @@ -475,6 +562,32 @@ bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) { // --- the model --- +namespace { + +// The atoms' centroid moved to the origin, and the radius, mass and coordinate error of what is left. +void finish_search_model(SearchModel &m) { + std::vector pos = ModelPositions(m.st.models[0]); + double error2 = 0; + for (const gemmi::Chain &ch : m.st.models[0].chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + error2 += 3.0 * a.b_iso / (8.0 * PI * PI); // the B from pLDDT, back to its rmsd^2 + m.error_A = std::sqrt(error2 / static_cast(pos.size())); + gemmi::Vec3 c; + for (const gemmi::Position &p : pos) + c += p; + c *= 1.0 / static_cast(pos.size()); + m.radius_A = 0; + for (gemmi::Position &p : pos) { + p = gemmi::Position(gemmi::Vec3(p) - c); + m.radius_A = std::max(m.radius_A, p.length()); + } + SetModelPositions(m.st.models[0], pos); + m.mass_da = DALTON_PER_RESIDUE * m.residues_kept; +} + +} // namespace + SearchModel PrepareSearchModel(const std::string &path) { SearchModel out; try { @@ -491,17 +604,24 @@ SearchModel PrepareSearchModel(const std::string &path) { 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) + std::vector kept; + for (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; }); + // 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. Each kept with + // its position in the whole prediction, which its predicted aligned error is indexed by. + if (!r.atoms.empty() && r.atoms[0].b_iso >= 70.0f) { + kept.push_back(std::move(r)); + out.residue_index.push_back(out.residues_read); + } + out.residues_read++; + } + ch.residues = std::move(kept); out.residues_kept += static_cast(ch.residues.size()); } std::erase_if(model.chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); }); @@ -511,24 +631,137 @@ SearchModel PrepareSearchModel(const std::string &path) { } // 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()); + finish_search_model(out); + return out; +} + +// --- domains from the predicted aligned error --- + +std::string PaePathFor(const std::string &model_path) { + const size_t slash = model_path.find_last_of("/\\"); + const std::string dir = slash == std::string::npos ? std::string() : model_path.substr(0, slash + 1); + const std::string name = model_path.substr(dir.size()); + const size_t at = name.rfind("-model_v"); + if (at == std::string::npos) + return {}; + const size_t dot = name.find('.', at); + const std::string version = name.substr(at + 8, dot == std::string::npos ? std::string::npos : dot - at - 8); + return dir + name.substr(0, at) + "-predicted_aligned_error_v" + version + ".json"; +} + +// A prediction cut into domains by its predicted aligned error before molecular replacement, as Oeffner et +// al. (2022) Acta Cryst. D78, 1303-1314 do; here by average linkage. +std::vector> PaeDomains(const std::string &pae_path, const SearchModel &model) { + // Two domains a prediction does not place relative to each other with confidence are this far apart in + // mean PAE (within one domain it is a few angstroms); a domain is at least this many residues. + constexpr double PAE_DOMAIN_A = 10.0; + constexpr int DOMAIN_MIN_RESIDUES = 40; + // The clustering starts from runs of this many consecutive residues kept. + constexpr int PIECE = 5; + std::vector> pae; + try { + std::ifstream in(pae_path); + if (!in) + return {}; + const nlohmann::json j = nlohmann::json::parse(in); + const nlohmann::json &m = j.is_array() ? j.at(0).at("predicted_aligned_error") : j.at("predicted_aligned_error"); + for (const auto &row : m) + pae.push_back(row.get>()); + } catch (const std::exception &) { + return {}; } - SetModelPositions(model, pos); - out.mass_da = DALTON_PER_RESIDUE * out.residues_kept; + const size_t n = pae.size(); + if (n != static_cast(model.residues_read) + || std::any_of(pae.begin(), pae.end(), [n](const std::vector &r) { return r.size() != n; })) + return {}; + // The pieces: runs of consecutive kept residues, at most PIECE long. + std::vector> members; + for (size_t i = 0; i < model.residue_index.size(); i++) { + const int r = model.residue_index[i]; + if (members.empty() || static_cast(members.back().size()) == PIECE || members.back().back() != r - 1) + members.push_back({}); + members.back().push_back(r); + } + // Average linkage on the symmetrised PAE: the two clusters with the lowest mean PAE between them merged + // until that is PAE_DOMAIN_A or more. + const size_t m = members.size(); + std::vector sum(m * m, 0.0); + for (size_t a = 0; a < m; a++) + for (size_t b = 0; b < m; b++) + for (int i : members[a]) + for (int k : members[b]) + sum[a * m + b] += (pae[i][k] + pae[k][i]) / 2.0; + std::vector active(m, 1); + while (true) { + double best = PAE_DOMAIN_A; + size_t ba = 0, bb = 0; + for (size_t a = 0; a < m; a++) + for (size_t b = a + 1; b < m; b++) + if (active[a] && active[b]) { + const double mean = sum[a * m + b] / static_cast(members[a].size() * members[b].size()); + if (mean < best) { + best = mean; + ba = a; + bb = b; + } + } + if (!(best < PAE_DOMAIN_A)) + break; + for (size_t c = 0; c < m; c++) { + sum[ba * m + c] += sum[bb * m + c]; + sum[c * m + ba] += sum[c * m + bb]; + } + members[ba].insert(members[ba].end(), members[bb].begin(), members[bb].end()); + active[bb] = 0; + } + std::vector> out; + for (size_t a = 0; a < m; a++) + if (active[a] && static_cast(members[a].size()) >= DOMAIN_MIN_RESIDUES) { + std::sort(members[a].begin(), members[a].end()); + out.push_back(members[a]); + } + std::stable_sort(out.begin(), out.end(), [](const std::vector &a, const std::vector &b) { + return a.size() > b.size(); + }); + return out; +} + +SearchModel DomainModel(const SearchModel &model, const std::vector &residues) { + SearchModel out = model; + out.residue_index.clear(); + out.residues_kept = 0; + size_t i = 0; + for (gemmi::Chain &ch : out.st.models[0].chains) { + std::vector kept; + for (gemmi::Residue &r : ch.residues) { + const int index = model.residue_index[i++]; + if (std::binary_search(residues.begin(), residues.end(), index)) { + kept.push_back(std::move(r)); + out.residue_index.push_back(index); + } + } + ch.residues = std::move(kept); + out.residues_kept += static_cast(ch.residues.size()); + } + std::erase_if(out.st.models[0].chains, [](const gemmi::Chain &ch) { return ch.residues.empty(); }); + finish_search_model(out); + // The domain's residues as ranges of the prediction's numbering, from 1. + std::string ranges; + for (size_t k = 0; k < residues.size(); k++) { + if (k > 0 && residues[k] == residues[k - 1] + 1 && k + 1 < residues.size() && residues[k + 1] == residues[k] + 1) + continue; + if (k > 0 && residues[k] == residues[k - 1] + 1) + ranges += fmt::format("-{}", residues[k] + 1); + else + ranges += fmt::format("{}{}", ranges.empty() ? "" : ",", residues[k] + 1); + } + out.description = fmt::format("PAE domain of {} residues ({})", out.residues_kept, ranges); return out; } @@ -610,7 +843,13 @@ MolecularTransform::MolecularTransform(const SearchModel &model, double d_min, s 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); +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() > 0 && MapToFPhiGPUFits(dc.grid)) + grid_ = MapToFPhiGPU(dc.grid); + else +#endif + grid_ = MapToFPhi(dc.grid); + d_min_ = d_min; blur_ = dc.blur; } @@ -844,39 +1083,6 @@ std::vector> KnownStructureFactors(const MolecularTransform return out; } -double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell, - const std::vector> &known, const gemmi::Mat33 &rotation, - const gemmi::Fractional &t) { - const size_t terms = setup.Terms(), nops = setup.ops; - const gemmi::Mat33 rt = rotation.transpose(); - const gemmi::Position x = cell.orthogonalize(t); - std::vector ic(terms), bin_sum(setup.bins, 0.0); - std::vector bin_n(setup.bins, 0); - for (size_t i = 0; i < terms; i++) { - std::complex f = known[i]; - for (size_t j = 0; j < nops; j++) { - const gemmi::Vec3 &sj = setup.s[i * nops + j]; - f += ft.At(rt.multiply(sj)) * setup.phase[i * nops + j] * std::polar(1.0, 2.0 * PI * sj.dot(x)); - } - ic[i] = std::norm(f); - bin_sum[setup.bin[i]] += ic[i]; - bin_n[setup.bin[i]]++; - } - double w = 0, sx = 0, sy = 0, sxx = 0, syy = 0, sxy = 0; - for (size_t i = 0; i < terms; i++) { - const double m = setup.multiplicity[i], a = ic[i] / (bin_sum[setup.bin[i]] / bin_n[setup.bin[i]]), - b = setup.e2[i]; - w += m; - sx += m * a; - sy += m * b; - sxx += m * a * a; - syy += m * b * b; - sxy += m * a * b; - } - const double vx = sxx - sx * sx / w, vy = syy - sy * sy / w; - return vx > 0 && vy > 0 ? (sxy - sx * sy / w) / std::sqrt(vx * vy) : 0.0; -} - TranslationPeak TranslationFunction(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::Mat33 &rotation, const std::vector> *known) { TranslationPeak out; @@ -944,6 +1150,246 @@ TranslationPeak TranslationFunction(const MolecularTransform &ft, const Translat return out; } +// --- the likelihood --- +// +// The log-likelihood gain of a placement over a random structure, from the Rice (acentric) and Woolfson +// (centric) distributions of E_obs given E_calc with sigma_A, and sigma_A from the model's expected +// coordinate error, following Read (2001) Acta Cryst. D57, 1373-1382; the gain judged as a Z-score over +// the same orientation at other positions, as the likelihood-enhanced translation function's Z-score of +// McCoy et al. (2005) Acta Cryst. D61, 458-464. + +namespace { + +// ln I0(x), x >= 0. Following Abramowitz & Stegun (1964) 9.8.1-9.8.2 +double log_bessel_i0(double x) { + if (x < 3.75) { + const double t = (x / 3.75) * (x / 3.75); + return std::log(1.0 + t * (3.5156229 + t * (3.0899424 + t * (1.2067492 + t * (0.2659732 + + t * (0.0360768 + t * 0.0045813)))))); + } + const double u = 3.75 / x; + return x - 0.5 * std::log(x) + std::log(0.39894228 + u * (0.01328592 + u * (0.00225319 + u * (-0.00157565 + + u * (0.00916281 + u * (-0.02057706 + u * (0.02635537 + u * (-0.01647633 + u * 0.00392377)))))))); +} + +double log_cosh(double x) { + x = std::fabs(x); + return x + std::log1p(std::exp(-2.0 * x)) - std::log(2.0); +} + +constexpr int LIKELIHOOD_NULL_POSITIONS = 100; +constexpr unsigned LIKELIHOOD_SEED = 20261012; + +} // namespace + +LikelihoodSetup MakeLikelihoodSetup(const SearchReflections &refl, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg) { + const gemmi::GroupOps gops = sg.operations(); + const gemmi::Mat33 frac_t = cell.frac.mat.transpose(); + LikelihoodSetup out; + out.ops = static_cast(gops.sym_ops.size()); + out.bins = refl.bins; + for (size_t i = 0; i < refl.hkl.size(); i++) { + const gemmi::Miller &h = refl.hkl[i]; + if (gops.is_systematically_absent(h)) + continue; + for (const gemmi::Op &op : gops.sym_ops) { + const gemmi::Miller k = rotated_index(op, h); + out.k.push_back(k); + out.s.push_back(frac_t.multiply(miller_vec(k))); + const double ht = (h[0] * op.tran[0] + h[1] * op.tran[1] + h[2] * op.tran[2]) + / static_cast(gemmi::Op::DEN); + out.phase.push_back(std::polar(1.0, 2.0 * PI * ht)); + for (int a = 0; a < 3; a++) + out.max_k[a] = std::max(out.max_k[a], std::abs(k[a])); + } + out.eo.push_back(std::sqrt(std::max(refl.e2[i], 0.0))); + out.epsilon.push_back(gops.epsilon_factor_without_centering(h)); + out.s2.push_back(cell.calculate_1_d2(h)); + out.centric.push_back(gops.is_reflection_centric(h)); + out.bin.push_back(refl.bin[i]); + } + return out; +} + +std::vector> CopyFactors(const SearchModel &model, const LikelihoodSetup &setup, + const gemmi::UnitCell &cell, const gemmi::Mat33 &rotation, + double d_min, size_t nthreads, const MolecularTransform *ft) { + std::vector> g(setup.s.size()); + if (ft != nullptr && d_min >= ft->DMin()) { + const gemmi::Mat33 rt = rotation.transpose(); + for (size_t e = 0; e < g.size(); e++) + g[e] = ft->At(rt.multiply(setup.s[e])) * setup.phase[e]; + return g; + } + gemmi::Structure st = model.st; + std::vector pos = ModelPositions(st.models[0]); + for (gemmi::Position &p : pos) + p = gemmi::Position(rotation.multiply(p)); + SetModelPositions(st.models[0], pos); + st.cell = cell; + 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); + const gemmi::FPhiGrid f1 = MapToFPhi(dc.grid); + for (size_t e = 0; e < g.size(); e++) { + int i = setup.k[e][0], j = setup.k[e][1], l = setup.k[e][2]; + 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 %= f1.nu; if (i < 0) i += f1.nu; + j %= f1.nv; if (j < 0) j += f1.nv; + const std::complex v = f1.data[i + static_cast(f1.nu) * (j + static_cast(f1.nv) * l)]; + g[e] = (conj ? std::conj(v) : v) * std::exp(dc.blur * setup.s[e].length_sq() / 4.0) * setup.phase[e]; + } + return g; +} + +void AddCopy(const LikelihoodSetup &setup, const std::vector> &g, const gemmi::Fractional &t, + std::vector> &f) { + // exp(2 pi i m t_a) for every index m an axis needs, so that each term is three products. + std::array>, 3> table; + const double ta[3] = {t.x, t.y, t.z}; + for (int a = 0; a < 3; a++) { + const int m = setup.max_k[a]; + table[a].resize(2 * m + 1); + for (int i = -m; i <= m; i++) + table[a][i + m] = std::polar(1.0, 2.0 * PI * i * ta[a]); + } + const size_t nops = setup.ops; + for (size_t i = 0; i < setup.Terms(); i++) + for (size_t j = 0; j < nops; j++) { + const size_t e = i * nops + j; + const gemmi::Miller &k = setup.k[e]; + f[i] += g[e] * table[0][k[0] + setup.max_k[0]] * table[1][k[1] + setup.max_k[1]] + * table[2][k[2] + setup.max_k[2]]; + } +} + +namespace { + +// E_calc^2 of each term: |f|^2 / epsilon over its bin's mean (0 where the bin has none). +std::vector calc_e2(const LikelihoodSetup &setup, const std::vector> &f) { + const size_t n = setup.Terms(); + std::vector bin_sum(setup.bins, 0.0); + std::vector bin_n(setup.bins, 0); + std::vector e2(n); + for (size_t i = 0; i < n; i++) { + e2[i] = std::norm(f[i]) / setup.epsilon[i]; + bin_sum[setup.bin[i]] += e2[i]; + bin_n[setup.bin[i]]++; + } + for (size_t i = 0; i < n; i++) { + const double mean = bin_sum[setup.bin[i]] / bin_n[setup.bin[i]]; + e2[i] = mean > 0 ? e2[i] / mean : 0.0; + } + return e2; +} + +// sigma_A's fall-off with resolution for the model's coordinate error, Luzzati's D, at each term. +std::vector error_factor(const LikelihoodSetup &setup, double error_A) { + std::vector out(setup.Terms()); + for (size_t i = 0; i < out.size(); i++) + out[i] = std::exp(-2.0 * PI * PI * error_A * error_A * setup.s2[i] / 3.0); + return out; +} + +double llg_sum(const LikelihoodSetup &setup, const std::vector &ec2, const std::vector &d, + double scale) { + double llg = 0; + for (size_t i = 0; i < setup.Terms(); i++) { + const double eo = setup.eo[i]; + const double sa = scale * d[i]; + const double v = 1.0 - sa * sa; + const double x = 2.0 * sa * eo * std::sqrt(ec2[i]) / v; + if (setup.centric[i]) + llg += -0.5 * std::log(v) - sa * sa * (eo * eo + ec2[i]) / (2.0 * v) + log_cosh(x / 2.0); + else + llg += -std::log(v) - sa * sa * (eo * eo + ec2[i]) / v + log_bessel_i0(x); + } + return llg; +} + +} // namespace + +double LogLikelihoodGain(const LikelihoodSetup &setup, const std::vector> &f, double error_A, + double scale) { + return llg_sum(setup, calc_e2(setup, f), error_factor(setup, error_A), scale); +} + +double RefinedLogLikelihoodGain(const LikelihoodSetup &setup, const std::vector> &f, + double error_A, double *scale) { + // Golden-section search for the scale in (0, 0.99), to a thousandth: the gain is 0 at 0 and has one + // maximum. + const std::vector ec2 = calc_e2(setup, f), d_error = error_factor(setup, error_A); + const double golden = (std::sqrt(5.0) - 1.0) / 2.0; + double a = 0.0, b = 0.99; + double c = b - golden * (b - a), d = a + golden * (b - a); + double fc = llg_sum(setup, ec2, d_error, c), fd = llg_sum(setup, ec2, d_error, d); + for (int iter = 0; iter < 14; iter++) { + if (fc > fd) { + b = d; d = c; fd = fc; + c = b - golden * (b - a); + fc = llg_sum(setup, ec2, d_error, c); + } else { + a = c; c = d; fc = fd; + d = a + golden * (b - a); + fd = llg_sum(setup, ec2, d_error, d); + } + } + if (scale != nullptr) + *scale = fc > fd ? c : d; + return std::max(fc, fd); +} + +double LikelihoodGain(const LikelihoodSetup &setup, const std::vector> &known, double known_llg, + const std::vector> &g, const gemmi::Fractional &t, double error_A, + double *llg, double *scale) { + std::vector> f = known; + AddCopy(setup, g, t, f); + const double all = RefinedLogLikelihoodGain(setup, f, error_A, scale); + if (llg != nullptr) + *llg = all; + return all - known_llg; +} + +LikelihoodPeak LikelihoodTranslationZ(const LikelihoodSetup &setup, const std::vector> &known, + double known_llg, const std::vector> &g, + const gemmi::Fractional &t, double error_A, size_t nthreads) { + LikelihoodPeak out; + out.gain = LikelihoodGain(setup, known, known_llg, g, t, error_A, &out.llg, &out.scale); + // The positions drawn from the raw 32-bit stream, which every standard library makes the same. + std::mt19937 rng(LIKELIHOOD_SEED); + std::vector where(LIKELIHOOD_NULL_POSITIONS); + for (gemmi::Fractional &w : where) { + const double x = rng() / 4294967296.0, y = rng() / 4294967296.0, z = rng() / 4294967296.0; + w = gemmi::Fractional(x, y, z); + } + std::vector null(where.size()); + const std::vector d_error = error_factor(setup, error_A); + ParallelFor(static_cast(where.size()), nthreads, [&](int i) { + std::vector> fr = known; + AddCopy(setup, g, where[i], fr); + null[i] = llg_sum(setup, calc_e2(setup, fr), d_error, out.scale); + }); + double mean = 0, s2 = 0; + for (double x : null) + mean += x; + mean /= static_cast(null.size()); + for (double x : null) + s2 += (x - mean) * (x - mean); + const double sd = std::sqrt(s2 / static_cast(null.size() - 1)); + out.z = sd > 0 ? (out.llg - mean) / sd : 0.0; + return out; +} + // --- translational NCS --- TranslationalNcs DetectTranslationalNcs(const std::vector &merged, const gemmi::UnitCell &cell, @@ -1186,6 +1632,168 @@ gemmi::Structure AssembleCopies(const std::vector &copies) { return out; } +std::vector SplitCopies(const gemmi::Structure &assembly, int copies) { + std::vector body(copies, assembly); + const size_t chains = assembly.models[0].chains.size(); + for (int k = 0; k < copies; k++) { + body[k].models[0].chains.clear(); + for (size_t i = 0; i < chains; i++) + if (static_cast(i * copies / chains) == k) + body[k].models[0].chains.push_back(assembly.models[0].chains[i]); + } + return body; +} + +std::optional SubgroupBasis(const gemmi::UnitCell &g, const gemmi::SpaceGroup &g_sg, + const gemmi::UnitCell &h, const gemmi::SpaceGroup &h_sg) { + constexpr double TOLERANCE = 0.01; // relative, on lengths, and on the cosines of the angles + const gemmi::GroupOps gops = g_sg.operations(); + const auto lattice = [&](const gemmi::Vec3 &v) { + for (const gemmi::Op::Tran &c : gops.cen_ops) { + const double d[3] = {v.x - static_cast(c[0]) / gemmi::Op::DEN, + v.y - static_cast(c[1]) / gemmi::Op::DEN, + v.z - static_cast(c[2]) / gemmi::Op::DEN}; + if (std::all_of(d, d + 3, [](double x) { return std::fabs(x - std::round(x)) < 1e-9; })) + return true; + } + return false; + }; + const double length[3] = {h.a, h.b, h.c}; + // Lattice vectors of g, in g's fractional coordinates, as long as each axis of h. + std::array, 3> axes; + for (int i = -4; i <= 4; i++) + for (int j = -4; j <= 4; j++) + for (int k = -4; k <= 4; k++) { + const gemmi::Vec3 v(i / 2.0, j / 2.0, k / 2.0); + if (!lattice(v)) + continue; + const double len = g.orthogonalize_difference(gemmi::Fractional(v)).length(); + for (int a = 0; a < 3; a++) + if (std::fabs(len - length[a]) < TOLERANCE * length[a]) + axes[a].push_back(v); + } + const double cosine[3] = {std::cos(h.alpha * PI / 180), std::cos(h.beta * PI / 180), std::cos(h.gamma * PI / 180)}; + const auto cos_between = [&](const gemmi::Vec3 &u, const gemmi::Vec3 &v) { + const gemmi::Position x = g.orthogonalize_difference(gemmi::Fractional(u)); + const gemmi::Position y = g.orthogonalize_difference(gemmi::Fractional(v)); + return x.dot(y) / (x.length() * y.length()); + }; + // The origin of h's setting in it: g's operators carried into h's axes, those of them h has shifted by o + // onto h's own - each translation t' + (1 - R) o the same as h's, modulo the lattice. + const gemmi::GroupOps hops = h_sg.operations(); + const auto origin = [&](const gemmi::Mat33 &m) -> std::optional { + const gemmi::Mat33 inv = m.inverse(); + struct Pair { gemmi::Mat33 r; gemmi::Vec3 t, t_h; }; + std::vector pairs; + for (const gemmi::Op &op : gops.all_ops_sorted()) { + gemmi::Mat33 r; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + r.a[i][j] = static_cast(op.rot[i][j]) / gemmi::Op::DEN; + const gemmi::Mat33 rh = inv.multiply(r).multiply(m); + const gemmi::Vec3 th = inv.multiply(gemmi::Vec3(op.tran[0], op.tran[1], op.tran[2]) / gemmi::Op::DEN); + for (const gemmi::Op &hop : hops.all_ops_sorted()) { + bool same = true; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + same = same && std::fabs(rh.a[i][j] - static_cast(hop.rot[i][j]) / gemmi::Op::DEN) < 1e-6; + if (same) + pairs.push_back({rh, th, gemmi::Vec3(hop.tran[0], hop.tran[1], hop.tran[2]) / gemmi::Op::DEN}); + } + } + for (int i = 0; i < 24; i++) + for (int j = 0; j < 24; j++) + for (int k = 0; k < 24; k++) { + const gemmi::Vec3 o(i / 24.0, j / 24.0, k / 24.0); + const bool all = std::all_of(pairs.begin(), pairs.end(), [&](const Pair &p) { + // Each of g's operators that lands on h's rotations has to land on one of h's operators. + const gemmi::Vec3 t = p.t + o - p.r.multiply(o); + for (const gemmi::Op &hop : hops.all_ops_sorted()) { + gemmi::Mat33 hr; + for (int a = 0; a < 3; a++) + for (int b = 0; b < 3; b++) + hr.a[a][b] = static_cast(hop.rot[a][b]) / gemmi::Op::DEN; + if (!hr.approx(p.r, 1e-6)) + continue; + const gemmi::Vec3 d = t - gemmi::Vec3(hop.tran[0], hop.tran[1], hop.tran[2]) / gemmi::Op::DEN; + if (std::fabs(d.x - std::round(d.x)) < 1e-6 && std::fabs(d.y - std::round(d.y)) < 1e-6 + && std::fabs(d.z - std::round(d.z)) < 1e-6) + return true; + } + return false; + }); + if (all) + return o; + } + return std::nullopt; + }; + std::optional best; + double best_error = INFINITY; + for (const gemmi::Vec3 &a : axes[0]) + for (const gemmi::Vec3 &b : axes[1]) + for (const gemmi::Vec3 &c : axes[2]) { + const double error = std::fabs(cos_between(b, c) - cosine[0]) + std::fabs(cos_between(a, c) - cosine[1]) + + std::fabs(cos_between(a, b) - cosine[2]); + const gemmi::Mat33 m(a.x, b.x, c.x, a.y, b.y, c.y, a.z, b.z, c.z); + if (error < 3 * TOLERANCE && m.determinant() > 0 && error < best_error) + if (const auto o = origin(m)) { + best_error = error; + best = SubgroupSetting{m, *o}; + } + } + return best; +} + +gemmi::Structure ExpandIntoSubgroup(const gemmi::Structure &placed, int copies, const gemmi::SpaceGroup &g_sg, + const SubgroupSetting &setting, const gemmi::UnitCell &h, + const gemmi::SpaceGroup &h_sg, int *copies_out) { + constexpr double SAME = 0.01; // fractional, between the centroids of two images + const gemmi::Mat33 to_h = setting.basis.inverse(); + const gemmi::GroupOps hops = h_sg.operations(); + std::vector out; + std::vector centroids; + for (const gemmi::Structure © : SplitCopies(placed, copies)) + for (const gemmi::Op &op : g_sg.operations().all_ops_sorted()) { + gemmi::Structure image = copy; + std::vector pos = ModelPositions(image.models[0]); + std::vector f; + gemmi::Vec3 c; + for (const gemmi::Position &p : pos) { + const gemmi::Fractional x = placed.cell.fractionalize(p); + const std::array y = op.apply_to_xyz({{x.x, x.y, x.z}}); + f.push_back(to_h.multiply(gemmi::Vec3(y[0], y[1], y[2])) + setting.origin); + c += f.back(); + } + c *= 1.0 / static_cast(f.size()); + const gemmi::Vec3 shift(-std::floor(c.x), -std::floor(c.y), -std::floor(c.z)); + const gemmi::Fractional centre(c + shift); + // An image the subgroup's own symmetry makes of one already kept is the same copy. + const bool seen = std::any_of(centroids.begin(), centroids.end(), [&](const gemmi::Fractional &k) { + for (const gemmi::Op &hop : hops.all_ops_sorted()) { + const std::array y = hop.apply_to_xyz({{k.x, k.y, k.z}}); + const double d[3] = {y[0] - centre.x, y[1] - centre.y, y[2] - centre.z}; + if (std::all_of(d, d + 3, [](double x) { return std::fabs(x - std::round(x)) < SAME; })) + return true; + } + return false; + }); + if (seen) + continue; + centroids.push_back(centre); + for (size_t i = 0; i < pos.size(); i++) + pos[i] = h.orthogonalize(gemmi::Fractional(f[i] + shift)); + SetModelPositions(image.models[0], pos); + out.push_back(std::move(image)); + } + gemmi::Structure st = AssembleCopies(out); + st.cell = h; + st.spacegroup_hm = h_sg.xhm(); + st.setup_cell_images(); + if (copies_out != nullptr) + *copies_out = static_cast(out.size()); + return st; +} + // The density of the Matthews coefficient over protein crystals: a Gumbel (extreme-value) distribution // - the family Kantardjieff & Rupp fit to this highly skewed peak - with their mode (2.34 A^3/Da) and // mean (2.69 A^3/Da) for 10471 protein crystal forms, which make its scale (mean - mode) / Euler's @@ -1222,19 +1830,17 @@ std::vector SameRotationGroups(const gemmi::SpaceGroup &sg) { // --- the search and the decision --- -MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, - Logger &logger, size_t nthreads) { +namespace { + +// The search and the decision for one search model. +MolecularReplacementResult search_model(const SearchModel &model, 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; - } + result.search_model = model.description; + auto t0 = std::chrono::steady_clock::now(); if (merges.empty() || merges[0].merged == nullptr) { result.failure_reason = "no merge to search"; return result; @@ -1246,12 +1852,18 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path 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); + logger.Info("Model prior: {}: {} of {} residues kept at pLDDT >= 70 ({:.1f} kDa, radius {:.1f} A)", + model.description, 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); + std::optional ft_made; + try { + ft_made.emplace(model, result.d_min_search, nthreads); + } catch (const std::exception &e) { // the GPU's FFT + result.failure_reason = std::string("the molecular transform: ") + e.what(); + return result; + } + const MolecularTransform &ft = *ft_made; #ifdef JFJOCH_USE_CUDA // The rotation and translation functions on the GPU where there is one; the rest of the search is // the same either way. @@ -1287,6 +1899,18 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path result.d_min_search, result.d_max_search); const std::vector pg = CartesianPointGroup(merge.space_group, cell); const std::vector rotations = RotationAsu(sphere, pg); + // The likelihood's reflections: down to LLG_D_MIN, or as far as the lowest-resolution LLG_TERMS of + // them reach where there are more - at least the search's own. + std::vector ds; + for (const MergedReflection &r : *merge.merged) + if (std::isfinite(r.I) && r.d <= result.d_max_search) + ds.push_back(r.d); + std::sort(ds.begin(), ds.end(), std::greater<>()); + double d_llg = std::max(LLG_D_MIN, data_d_min); + if (ds.size() > LLG_TERMS) + d_llg = std::max(d_llg, std::min(result.d_min_search, ds[LLG_TERMS - 1])); + const SearchReflections lrefl = MakeSearchReflections(*merge.merged, cell, merge.space_group, d_llg, + result.d_max_search); #ifdef JFJOCH_USE_CUDA const std::vector rf = gpu ? gpu->RotationFunction(refl, cell, merge.space_group, rotations) : RotationFunction(ft, refl, cell, merge.space_group, rotations, nthreads); @@ -1438,6 +2062,11 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path by_tf.insert(by_tf.begin(), first.begin(), first.end()); } const TranslationPeak *tf = tfs.data(); + const LikelihoodSetup lsetup = MakeLikelihoodSetup(lrefl, cell, sg); + const std::vector> no_known(lsetup.Terms(), 0.0); + // The further copies' peaks are ranked by the likelihood on the search's own reflections, read off + // the molecular transform. + const LikelihoodSetup rsetup = MakeLikelihoodSetup(refl, cell, sg); // --- the first copy: the best few through packing, the rigid body and R-free --- RigidBodyEngines engines; @@ -1450,6 +2079,10 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path std::vector copies; std::vector copy_rotation; std::vector copy_position; + // CopyFactors of each, for the likelihood and for ranking (rsetup), as they are needed. + std::vector>> copy_g, copy_rg; + // CopyFactors of the orientations searched, as they are needed. + std::vector>> rot_g(rots.size()); // Down the list until RESCORED_PER_GROUP placements that pack have been refined: a placement // clashing with its own mates over twice the limit is not refined - the rigid body moves a // search placement by a few degrees and a fraction of an angstrom, and the clash fraction moves @@ -1475,6 +2108,7 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path } else if (!fits) { refined++; } + // The refined placement's likelihood, judged against the same orientation elsewhere. // The refined placement with the lowest R-free; the best translation-function one where // none was refined. if (!have || (std::isfinite(score.r_free) && !(cand.r_free <= score.r_free))) { @@ -1493,11 +2127,25 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path copies = {*placed}; copy_rotation = {moved.multiply(rots[i])}; copy_position = {centroid(*placed)}; + copy_g = {{}}; + copy_rg = {{}}; } glog.Info("Model prior: {:<12} RF Z {:5.1f} TF Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", sg.xhm(), rf_sd > 0 ? (rfs[i] - rf_mean) / rf_sd : 0.0, tf[i].z, 100 * clash, score.r_work, score.r_free); } + // The likelihood of the placement kept, judged against the same orientation elsewhere - only where + // its translation function's Z-score does not already stand out. + if (have && std::isfinite(cand.r_free) && cand.tf_z < TF_Z_PLACED) { + copy_g[0] = CopyFactors(model, lsetup, cell, copy_rotation[0], d_llg, nthreads, &ft); + const LikelihoodPeak lk = LikelihoodTranslationZ(lsetup, no_known, 0.0, copy_g[0], copy_position[0], + model.error_A, nthreads); + cand.llg = lk.llg; + cand.llg_z = lk.z; + cand.copy_llg_z = {lk.z}; + glog.Info("Model prior: {:<12} the placement kept: LLG {:.1f}, Z {:.1f} over its orientation elsewhere", + sg.xhm(), lk.llg, lk.z); + } // --- further copies, one at a time, with the copies placed as a known partial structure --- // Each copy is refined on its own and all of them scored together. @@ -1515,6 +2163,8 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path &fixed.models[0]); copy_rotation[c] = moved.multiply(copy_rotation[c]); copy_position[c] = centroid(copies[c]); + copy_g[c].clear(); + copy_rg[c].clear(); } const gemmi::Structure assembly = AssembleCopies(copies); const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol}; @@ -1555,19 +2205,43 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path struct Trial { size_t i = 0; double clash = 1; - double cc = 0; + double gain = 0; + LikelihoodPeak refined; + std::vector> g; PlacementScore score{NAN, NAN, NAN, NAN}; gemmi::Structure copy, assembly; gemmi::Mat33 rotation; }; std::optional best; const int n = cand.copies + 1; - // The highest peaks that do not sit on a copy already placed (or a mate of one), ranked by - // the correlation of the assembly with E^2, and the best few of those through the rigid body. + // The highest peaks that do not sit on a copy already placed (or a mate of one), ranked by the + // likelihood gain of the copy beside the known ones, and the best few of those through the + // rigid body. A copy that only repeats a known one - moved by one of the group's origin shifts, + // it doubles F on half the reflections and cancels it on the other half - gains nothing. + std::vector> rknown(rsetup.Terms(), 0.0); + for (size_t c = 0; c < copies.size(); c++) { + if (copy_rg[c].empty()) + copy_rg[c] = CopyFactors(model, rsetup, cell, copy_rotation[c], result.d_min_search, nthreads, &ft); + AddCopy(rsetup, copy_rg[c], copy_position[c], rknown); + } + double known_scale = 0; + const double rknown_llg = RefinedLogLikelihoodGain(rsetup, rknown, model.error_A, &known_scale); + // Down the highest peaks until COPY_TRIALS that gain likelihood - a copy that lowers it is not + // one, and a peak on a copy already placed (or a mate of one) does - and pack. The gain that + // ranks them is taken at the known copies' sigma_A. std::vector trials; - for (size_t k = 0; k < by_tf2.size() && trials.size() < COPY_TRIALS; k++) { + for (size_t k = 0; k < by_tf2.size() && k < 2 * COPY_TRIALS && trials.size() < COPY_TRIALS; k++) { Trial trial; trial.i = by_tf2[k]; + std::vector> &g = trial.i < rots.size() ? rot_g[trial.i] + : copy_rg[trial.i - rots.size()]; + if (g.empty()) + g = CopyFactors(model, rsetup, cell, rots2[trial.i], result.d_min_search, nthreads, &ft); + std::vector> f = rknown; + AddCopy(rsetup, g, tf2[trial.i].t, f); + trial.gain = LogLikelihoodGain(rsetup, f, model.error_A, known_scale) - rknown_llg; + if (!(trial.gain > 0)) + continue; trial.copy = PlaceSearchModel(model, rots2[trial.i], tf2[trial.i].t, cell, sg); std::vector all = copies; all.push_back(trial.copy); @@ -1576,14 +2250,14 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path continue; trials.push_back(std::move(trial)); } - ParallelFor(static_cast(trials.size()), nthreads, [&](int k) { - trials[k].cc = AssemblyCorrelation(ft, *setup2, cell, known, rots2[trials[k].i], - tf2[trials[k].i].t); - }); - std::stable_sort(trials.begin(), trials.end(), [](const Trial &a, const Trial &b) { return a.cc > b.cc; }); + std::stable_sort(trials.begin(), trials.end(), + [](const Trial &a, const Trial &b) { return a.gain > b.gain; }); if (trials.size() > RESCORED_PER_GROUP) trials.resize(RESCORED_PER_GROUP); const gemmi::Structure others = AssembleCopies(copies); + // The known copies on the likelihood's own reflections, made where a trial needs its Z-score. + std::optional>> lknown; + double known_llg = 0; for (Trial &trial : trials) { std::vector all = copies; gemmi::Mat33 moved; @@ -1595,11 +2269,26 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path const PlacementScore held{NAN, NAN, cand.k_sol, cand.b_sol}; trial.score = score_placement(trial.assembly, sg, *merge.merged, d_min, nthreads, &scorers, &held); trial.rotation = moved.multiply(rots2[trial.i]); - glog.Info("Model prior: {:<12} copy {} TF Z {:5.1f} CC {:.3f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", - sg.xhm(), n, tf2[trial.i].z, trial.cc, 100 * trial.clash, trial.score.r_work, - trial.score.r_free); + trial.refined.gain = trial.refined.llg = trial.refined.z = NAN; + if (tf2[trial.i].z < TF_Z_PLACED) { + if (!lknown) { + lknown.emplace(lsetup.Terms(), 0.0); + for (size_t c = 0; c < copies.size(); c++) { + if (copy_g[c].empty()) + copy_g[c] = CopyFactors(model, lsetup, cell, copy_rotation[c], d_llg, nthreads, &ft); + AddCopy(lsetup, copy_g[c], copy_position[c], *lknown); + } + known_llg = RefinedLogLikelihoodGain(lsetup, *lknown, model.error_A); + } + trial.g = CopyFactors(model, lsetup, cell, trial.rotation, d_llg, nthreads, &ft); + trial.refined = LikelihoodTranslationZ(lsetup, *lknown, known_llg, trial.g, centroid(trial.copy), + model.error_A, nthreads); + } + glog.Info("Model prior: {:<12} copy {} TF Z {:5.1f} LLG gain {:7.1f}, refined {:7.1f} Z {:5.1f} clash {:4.1f}% R-work {:.3f} R-free {:.3f}", + sg.xhm(), n, tf2[trial.i].z, trial.gain, trial.refined.gain, + trial.refined.z, 100 * trial.clash, trial.score.r_work, trial.score.r_free); // The first, in that order, that stands out, packs and lowers R-free. - if (tf2[trial.i].z >= TF_Z_PLACED && trial.clash <= CLASH_MAX + if ((tf2[trial.i].z >= TF_Z_PLACED || trial.refined.z >= TF_Z_PLACED) && trial.clash <= CLASH_MAX && trial.score.r_free <= cand.r_free - COPY_R_FREE_GAIN) { best = trial; break; @@ -1610,9 +2299,14 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path copies.push_back(best->copy); copy_rotation.push_back(best->rotation); copy_position.push_back(centroid(best->copy)); + copy_g.push_back(best->g); + copy_rg.push_back({}); cand.copies = n; cand.copy_tf_z.push_back(tf2[best->i].z); cand.tf_z = std::max(cand.tf_z, tf2[best->i].z); + cand.copy_llg_z.push_back(best->refined.z); + cand.llg_z = std::fmax(cand.llg_z, best->refined.z); + cand.llg = best->refined.llg; // NaN where the translation function stood out cand.clash_fraction = best->clash; cand.r_work = best->score.r_work; cand.r_free = best->score.r_free; @@ -1642,7 +2336,10 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path cand.verdict = "MATTHEWS"; else if (cand.clash_fraction > CLASH_MAX) cand.verdict = "CLASH"; - else if (!(cand.tf_z >= TF_Z_PLACED)) + // The likelihood's Z-score stands in for the translation function's only where the copies placed + // are as many as the cell most probably holds: beside content the search has not found, a copy's + // likelihood over random positions of its own is no measure of the placement. + else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely))) cand.verdict = "NOT_PLACED"; else cand.verdict = "PLACED"; @@ -1660,9 +2357,9 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path std::vector found; for (auto &r : running) found.push_back(r.get()); - // Every group the best one leads searched again with the best one's copies' orientations tried - // first: the groups share the rotations, so a copy the search placed in one is in the same - // orientation in the others, and one group's search failing - next to copies related by a + // With tNCS, every group that placed fewer copies than the best one searched again with the best one's + // copies' orientations tried first: the groups share the rotations, so a copy the search placed in one is in + // the same orientation in the others, and one group's search failing - next to copies related by a // translation, the first copy alone fits no better than at random, and the copy the search keeps // can be the wrong one - is then not taken for evidence against it. The better of the two searches // is kept. @@ -1670,10 +2367,11 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path for (size_t g = 0; g < found.size(); g++) if (found[g].cand.verdict == "PLACED" && (lead < 0 || found[g].cand.r_free < found[lead].cand.r_free)) lead = static_cast(g); - if (lead >= 0) { + if (lead >= 0 && result.tncs.present) { std::vector> again(found.size()); for (size_t g = 0; g < found.size(); g++) - if (static_cast(g) != lead && found[g].cand.vm >= VM_MIN) + if (static_cast(g) != lead && found[g].cand.vm >= VM_MIN + && found[g].cand.copies < found[lead].cand.copies) again[g] = std::async(std::launch::async, search_group, g, found[lead].copy_rotation); for (size_t g = 0; g < found.size(); g++) if (again[g].valid()) { @@ -1754,12 +2452,24 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path if (static_cast(i) == best || c.merge != merge || c.clash_fraction > CLASH_MAX) continue; // With tNCS one copy alone fits no better than at random, so a group whose search stopped at - // fewer copies than the chosen one's has not been shown unable to hold them. + // fewer copies than the chosen one's has not been shown unable to hold them - where the tNCS + // vector makes up the groups' screw translations, or where they are each other's enantiomorph. if (result.tncs.present && c.copies < result.candidates[best].copies - && TncsRelatesGroups(c.space_group, result.candidates[best].space_group, result.tncs.vector)) { - result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS", - c.space_group.xhm(), result.candidates[best].space_group.xhm()); - return -2; + && (TncsRelatesGroups(c.space_group, result.candidates[best].space_group, result.tncs.vector) + || enantiomorphs(c.space_group, result.candidates[best].space_group))) { + // Unless the chosen copies, carried into this group as they are and refined, build a crystal + // that fits clearly worse: where the tNCS vector does make up the difference between the groups, + // they build the same one. + const MRCandidate &b = result.candidates[best]; + const PlacementScore carried = carry_copies(*b.placed, b.copies, c.space_group, *merges[merge].merged, + result.d_min_score, nthreads, &scorers, logger); + logger.Info("Model prior: {}'s {} copies carried into {}: R-free {:.3f} against {:.3f}", + b.space_group.xhm(), b.copies, c.space_group.xhm(), carried.r_free, b.r_free); + if (!(carried.r_free - b.r_free >= R_FREE_MARGIN)) { + result.decision = fmt::format("UNDECIDED: {} placed fewer copies than {} beside tNCS", + c.space_group.xhm(), b.space_group.xhm()); + return -2; + } } if (!(c.r_free - result.candidates[best].r_free >= R_FREE_MARGIN)) { result.decision = fmt::format("UNDECIDED: {} and {} within {:.2f} in R-free", @@ -1798,6 +2508,48 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path } else if (answer == -1) { result.decision = "NOT_PLACED: no group of the answer's merge holds a placed model"; } + // The twin comparison, where the answer's merge chose and the twin hypothesis' merge is there: on that + // merge, against its data as twinned by its law, the chosen placement with all its symmetry images - the + // group's crystal - and the model searched for afresh in the subgroup, as a crystal of the subgroup. A + // crystal of the group twinned by the law reads like one of the subgroup, and its placement read in the + // subgroup fits as well as it fits in the group; a crystal of the subgroup has copies the group's + // symmetry does not relate, which only a search in the subgroup places. The subgroup is chosen where its + // own placement fits the twinned data better by R_FREE_MARGIN in R-free - on reflections none of the + // copies, nor the fraction, was fitted to - and passes the same choice as any merge's (its null, and a + // lead over the subgroup's other groups). + if (result.chosen >= 0 && result.candidates[result.chosen].merge == 0 && merges.size() > 1 + && merges[1].twin_law) { + const MRCandidate g = result.candidates[result.chosen]; + const MRMerge &tm = merges[1]; + const gemmi::UnitCell h_cell(tm.cell); + if (const auto setting = SubgroupBasis(gemmi::UnitCell(merges[0].cell), g.space_group, h_cell, tm.space_group)) { + const gemmi::Op &law = *tm.twin_law; + int n = 0; + const gemmi::Structure expanded = ExpandIntoSubgroup(*g.placed, g.copies, g.space_group, *setting, h_cell, + tm.space_group, &n); + const PlacementScore as_group = score_placement(expanded, tm.space_group, *tm.merged, result.d_min_score, + nthreads, &scorers, nullptr, &law); + result.candidates[result.chosen].r_free_twin = as_group.r_free_twin; + const std::string decision = result.decision; + search_merge(1); + const int h = choose_in(1); + if (h >= 0) { + MRCandidate &c = result.candidates[h]; + const PlacementScore t = score_placement(*c.placed, c.space_group, *tm.merged, result.d_min_score, + nthreads, &scorers, nullptr, &law); + c.r_work_twin = t.r_work_twin; + c.r_free_twin = t.r_free_twin; + c.twin_fraction = t.twin_fraction; + logger.Info("Model prior: twin comparison on the {} merge, twinned by {}: {} read as its {} copies " + "there R-free {:.3f}; searched for in {}, {} copies, R-free {:.3f} at a fraction {:.2f}", + tm.space_group.xhm(), law.as_hkl().triplet('h'), g.space_group.xhm(), n, + as_group.r_free_twin, c.space_group.xhm(), c.copies, t.r_free_twin, t.twin_fraction); + if (t.r_free_twin <= as_group.r_free_twin - R_FREE_MARGIN) + result.chosen = h; + } + result.decision = decision; + } + } if (result.chosen >= 0) { const MRCandidate &c = result.candidates[result.chosen]; const gemmi::SpaceGroup &data = merges[0].space_group; @@ -1820,6 +2572,95 @@ MolecularReplacementResult RunMolecularReplacement(const std::string &model_path return result; } +} // namespace + +MolecularReplacementResult RunMolecularReplacement(const std::string &model_path, const std::vector &merges, + Logger &logger, size_t nthreads) { + const auto t0 = std::chrono::steady_clock::now(); + SearchModel model = PrepareSearchModel(model_path); + if (!model.failure_reason.empty()) { + MolecularReplacementResult result; + result.model_path = model_path; + result.residues_read = model.residues_read; + result.failure_reason = model.failure_reason; + return result; + } + model.description = "the whole model"; + const std::string pae_path = PaePathFor(model_path); + const std::vector> domains = pae_path.empty() ? std::vector>() + : PaeDomains(pae_path, model); + const double t_model = seconds_since(t0); + if (!pae_path.empty()) + logger.Info("Model prior: {} rigid domain(s) in the predicted aligned error {}", domains.size(), pae_path); + + MolecularReplacementResult result = search_model(model, merges, logger, nthreads); + result.model_path = model_path; + result.pae_domains = static_cast(domains.size()); + result.t_model += t_model; + // 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. + std::vector logs; + for (size_t i = 0; i < domains.size(); i++) + logs.push_back(Logger::Buffered()); + std::vector> 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); + })); + std::vector found; + for (size_t i = 0; i < domains.size(); i++) { + MolecularReplacementResult r = running[i].get(); + logs[i].ReplayInto(logger); + 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; + for (size_t i = 0; i < found.size(); i++) { + if (found[i].chosen < 0) { + agree = false; + continue; + } + if (best >= 0 && found[i].candidates[found[i].chosen].space_group.xhm() + != found[best].candidates[found[best].chosen].space_group.xhm()) + agree = false; + if (best < 0 || found[i].candidates[found[i].chosen].r_free < found[best].candidates[found[best].chosen].r_free) + best = static_cast(i); + } + // A domain is part of the model: it may confirm the answer's group, not change it. + if (best >= 0 && agree && found[best].decision != "AGREES") { + 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_path = model_path; + 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; + } else { + result.t_domains = domain_seconds; + if (!found.empty()) + result.decision = "UNDECIDED: the PAE domains placed disagree"; + } + logger.Info("Model prior: with the PAE domains searched by themselves: {}{}", result.decision, + result.chosen >= 0 ? " - " + result.candidates[result.chosen].space_group.xhm() + ", " + + result.search_model : ""); + } + return result; +} + void CheckEnantiomorphHand(MolecularReplacementResult &result, const std::vector &merged, const gemmi::SpaceGroup &data_space_group, Logger &logger, size_t nthreads) { const MRCandidate &c = result.candidates[result.chosen]; diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index 30db4f53e..d775e8b09 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -23,6 +23,7 @@ #include #include +#include #include #include @@ -43,10 +44,25 @@ struct SearchModel { 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 + double error_A = 0; // rms of the coordinate errors the pLDDTs predict, over the atoms kept + std::vector residue_index; // of each residue kept, in model order: its position in the prediction + std::string description; // "the whole model", or which domain std::string failure_reason; }; SearchModel PrepareSearchModel(const std::string &path); +// The rigid domains of a prediction from its predicted aligned error (PAE, the AlphaFold DB's JSON: a matrix +// "predicted_aligned_error" over the residues of the whole prediction): the kept residues clustered by average +// linkage on the symmetrised PAE until every two clusters are a mean PAE_DOMAIN_A apart; clusters of at least +// DOMAIN_MIN_RESIDUES, largest first, each a list of positions in the prediction. Empty where the file cannot +// be read or does not match the model. +std::vector> PaeDomains(const std::string &pae_path, const SearchModel &model); +// The search model of one domain: the residues of `model` at those positions, recentred. +SearchModel DomainModel(const SearchModel &model, const std::vector &residues); +// The PAE file the AlphaFold DB puts beside a model: AF--F1-model_v.cif -> +// AF--F1-predicted_aligned_error_v.json in the same directory; empty where the name does not follow it. +std::string PaePathFor(const std::string &model_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. @@ -75,11 +91,13 @@ public: const gemmi::FPhiGrid &Grid() const { return grid_; } double Box() const { return box_; } double Blur() const { return blur_; } + double DMin() const { return d_min_; } private: gemmi::FPhiGrid grid_; // half-l, on a P1 box `box_` A wide double box_ = 0; double blur_ = 0; + double d_min_ = 0; }; // The observed reflections the search scores against: E^2 in resolution bins, with each unique @@ -156,15 +174,62 @@ TranslationPeak TranslationFunction(const MolecularTransform &ft, const Translat const gemmi::Mat33 &rotation, const std::vector> *known = nullptr); -// The correlation of E^2 with |F|^2 of the known copies (KnownStructureFactors) and one more copy turned by -// `rotation` with its centroid at `t`, over the terms of `setup`, each weighted by its multiplicity and each -// bin's |F|^2 over its mean: the score the translation function's peaks for a further copy are ranked by. -// Unlike the translation function, which is linear in |F|^2, it does not reward a copy that only repeats -// the known ones - their own copy moved by one of the group's origin shifts doubles F on half the -// reflections and cancels it on the other half, which the overlap with E^2 counts as a gain. -double AssemblyCorrelation(const MolecularTransform &ft, const TranslationSetup &setup, const gemmi::UnitCell &cell, - const std::vector> &known, const gemmi::Mat33 &rotation, - const gemmi::Fractional &t); +// --- The likelihood a placement is judged by --- +// The terms it is summed over: the reflections of a SearchReflections not absent in a group, each with +// its E_obs, whether it is centric, its epsilon and 1/d^2, and for every operator j of the group h R_j (as +// a Cartesian vector and as indices) and exp(2 pi i h.t_j), as in TranslationSetup. +struct LikelihoodSetup { + int ops = 0; + std::vector s; // terms x ops + std::vector k; // terms x ops + std::vector> phase; // terms x ops + std::vector eo, epsilon, s2; // terms + std::vector centric; // terms + std::vector bin; // terms + int bins = 0; + int max_k[3] = {0, 0, 0}; // the largest |index| of h R_j along each axis + size_t Terms() const { return eo.size(); } +}; +LikelihoodSetup MakeLikelihoodSetup(const SearchReflections &refl, const gemmi::UnitCell &cell, + const gemmi::SpaceGroup &sg); + +// g_j = G_j(h) exp(2 pi i h.t_j) of one copy in orientation `rotation`, terms x ops: the copy with its +// centroid at t adds sum_j g_j exp(2 pi i (h R_j).t) to F(h). G_j(h) = F1(h R_j), F1 the transform of +// the copy in `cell` with its centroid at the origin - one FFT of the cell's grid to `d_min`, read off at +// the lattice points themselves, so exact at any resolution. Where `ft` is given and reaches `d_min`, read +// off it instead, interpolated as the rotation and translation functions read it. +std::vector> CopyFactors(const SearchModel &model, const LikelihoodSetup &setup, + const gemmi::UnitCell &cell, const gemmi::Mat33 &rotation, + double d_min, size_t nthreads, const MolecularTransform *ft = nullptr); +// `f` += that copy's F with its centroid at the fractional position t. +void AddCopy(const LikelihoodSetup &setup, const std::vector> &g, const gemmi::Fractional &t, + std::vector> &f); + +// The log-likelihood gain of the model structure factors `f` (one per term) over a random structure, +// with sigma_A = scale exp(-2 pi^2 error^2 / (3 d^2)): E_calc are |f| normalised in the resolution bins. +double LogLikelihoodGain(const LikelihoodSetup &setup, const std::vector> &f, double error_A, + double scale); +// The same with the scale refined to its maximum (returned in `scale`). +double RefinedLogLikelihoodGain(const LikelihoodSetup &setup, const std::vector> &f, + double error_A, double *scale = nullptr); + +// How the likelihood of one more copy in orientation `rotation` at `t`, beside the `known` structure factors +// (zero for none), stands out over the same copy in the same orientation at random positions: the gain +// with the scale refined at t, and its Z-score over the gains at those positions with that scale held. +struct LikelihoodPeak { + double llg = 0; // of all the copies, the known and this one + double gain = 0; // over the known ones alone, at their own refined scale + double z = 0; + double scale = 0; +}; +// `g`: CopyFactors of the copy in that orientation; `known_llg`: RefinedLogLikelihoodGain of `known` alone. +LikelihoodPeak LikelihoodTranslationZ(const LikelihoodSetup &setup, const std::vector> &known, + double known_llg, const std::vector> &g, + const gemmi::Fractional &t, double error_A, size_t nthreads); +// The gain alone: of the known structure factors and the copy over the known ones (`llg`: of all of them). +double LikelihoodGain(const LikelihoodSetup &setup, const std::vector> &known, double known_llg, + const std::vector> &g, const gemmi::Fractional &t, double error_A, + double *llg = nullptr, double *scale = nullptr); // 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, @@ -177,6 +242,28 @@ double ClashFraction(const gemmi::Structure &placed, size_t copies = 1); // Copies of the model as one structure: their chains in order, renamed so that no two share a name. gemmi::Structure AssembleCopies(const std::vector &copies); +// The copies of an assembly (AssembleCopies), its chains split evenly in order. +std::vector SplitCopies(const gemmi::Structure &assembly, int copies); + +// The setting of a subgroup's merge (cell `h`, group `h_sg`) in the cell `g` of a group `g_sg`: `basis` has h's +// axes as columns, lattice vectors of g (with g's centring) in g's fractional coordinates, right-handed and +// matching h's lengths and angles to 1%; `origin` is where g's origin lies in h's fractional coordinates, so +// that the operators of g with h's rotations are h's own: x_h = basis^-1 x_g + origin. Empty where there is no +// such setting. +struct SubgroupSetting { + gemmi::Mat33 basis; + gemmi::Vec3 origin; +}; +std::optional SubgroupBasis(const gemmi::UnitCell &g, const gemmi::SpaceGroup &g_sg, + const gemmi::UnitCell &h, const gemmi::SpaceGroup &h_sg); + +// The `copies` of a placement in group `g_sg` and every symmetry image of them that the subgroup `h_sg` does +// not make itself, as independent copies in the subgroup's cell `h`, in its setting (SubgroupBasis); their +// number in `copies_out`. +gemmi::Structure ExpandIntoSubgroup(const gemmi::Structure &placed, int copies, const gemmi::SpaceGroup &g_sg, + const SubgroupSetting &setting, const gemmi::UnitCell &h, + const gemmi::SpaceGroup &h_sg, int *copies_out); + // The probabilities of 1, 2, ... copies of a model of `mass_da` in the asymmetric unit of a cell of // `asu_volume` A^3, from the distribution of the Matthews coefficient over protein crystals; one entry // per copy number whose Vm is at least 1.5 A^3/Da. @@ -209,6 +296,7 @@ struct MRMerge { const std::vector *merged = nullptr; UnitCell cell; gemmi::SpaceGroup space_group; // the group it was merged in + std::optional twin_law; // the twin hypothesis' merge: the law it carries }; // One space group's best placement: its first copy and every further one the search added. @@ -222,10 +310,18 @@ struct MRCandidate { double rf_z = 0; // of the first copy's orientation, over the rotation function double tf_z = 0; // the highest of the copies' std::vector copy_tf_z; // each copy's, the later ones beside the earlier + // The likelihood (LikelihoodTranslationZ), after the rigid body, taken only for a copy whose translation + // function does not stand out (NaN otherwise): of all the copies with the last such one, and the highest + // of their Z-scores, and each copy's. + double llg = NAN; + double llg_z = NAN; + std::vector copy_llg_z; double clash_fraction = 0; // of all the copies double r_work = NAN, r_free = NAN; // after the rigid body, to d_min_score double k_sol = NAN, b_sol = NAN; // the bulk solvent R was computed with double null_sigma = NAN; // R-work above its random-orientation null, in sd + // On the twin hypothesis' merge, against the data as twinned by its law (the twin comparison). + double r_work_twin = NAN, r_free_twin = NAN, twin_fraction = NAN; std::string verdict; // PLACED, or why not std::shared_ptr placed; // every copy }; @@ -236,6 +332,8 @@ struct MolecularReplacementResult { std::string model_path; 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 + 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; double d_min_score = 0; // the limit R-work / R-free are computed to std::vector candidates; @@ -261,6 +359,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) }; // Search `merges` (the answer's first) for the predicted model at `model_path`, in every space group with diff --git a/image_analysis/structure_refinement/MolecularReplacementGPU.cu b/image_analysis/structure_refinement/MolecularReplacementGPU.cu index 8a062fee8..82a44f776 100644 --- a/image_analysis/structure_refinement/MolecularReplacementGPU.cu +++ b/image_analysis/structure_refinement/MolecularReplacementGPU.cu @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -260,6 +261,25 @@ __global__ void translation_stats(const float *map, size_t n, double *sum, doubl } } +// cuFFT's r2c output of a map stored u fastest ([w][v][u/2+1], u halved) taken to MapToFPhi's half-l +// layout (u fastest, w halved), conjugated and scaled as gemmi does. +__global__ void half_u_to_half_w(const cufftComplex *in, int nu, int nv, int nw, float norm, float2 *out) { + const int nwh = nw / 2 + 1, nuh = nu / 2 + 1; + const size_t n = static_cast(nu) * nv * nwh; + const size_t i = static_cast(blockIdx.x) * blockDim.x + threadIdx.x; + if (i >= n) + return; + const int u = static_cast(i % nu), v = static_cast((i / nu) % nv), l = static_cast(i / (static_cast(nu) * nv)); + cufftComplex x; + if (u < nuh) { + x = in[u + static_cast(nuh) * (v + static_cast(nv) * l)]; + } else { // F(-h) = conj F(h) + x = in[(nu - u) + static_cast(nuh) * ((nv - v) % nv + static_cast(nv) * ((nw - l) % nw))]; + x.y = -x.y; + } + out[i] = make_float2(x.x * norm, -x.y * norm); +} + Rot transposed(const gemmi::Mat33 &m) { Rot r; for (int i = 0; i < 3; i++) @@ -292,6 +312,51 @@ std::vector bin_counts(const std::vector &bin, int bins) { } // namespace +bool MapToFPhiGPUFits(const gemmi::Grid &map) { + const size_t points = static_cast(map.nu) * map.nv * map.nw; + const size_t half_u = static_cast(map.nu / 2 + 1) * map.nv * map.nw; + const size_t half_w = static_cast(map.nu) * map.nv * (map.nw / 2 + 1); + const size_t bytes = points * sizeof(float) + (half_u + half_w) * sizeof(cufftComplex); + size_t free_bytes = 0, total = 0; + cuda_err(cudaMemGetInfo(&free_bytes, &total)); + return bytes <= total / 4; +} + +gemmi::FPhiGrid MapToFPhiGPU(const gemmi::Grid &map) { + // One at a time: concurrent searches (the PAE domains') would each hold a large grid on the card. + static std::mutex m; + std::lock_guard lock(m); + gemmi::FPhiGrid hkl; + hkl.unit_cell = map.unit_cell; + hkl.spacegroup = map.spacegroup; + hkl.axis_order = map.axis_order; + hkl.half_l = true; + hkl.set_size_without_checking(map.nu, map.nv, map.nw / 2 + 1); + const float norm = static_cast(map.unit_cell.volume / map.point_count()); + CudaStream stream; + auto d_map = upload(map.data, stream); + CudaDevicePtr d_half(static_cast(map.nu / 2 + 1) * map.nv * map.nw); + CudaDevicePtr d_out(hkl.data.size()); + cufftHandle plan; + cufft_err(cufftPlan3d(&plan, map.nw, map.nv, map.nu, CUFFT_R2C)); + try { + cufft_err(cufftSetStream(plan, stream)); + cufft_err(cufftExecR2C(plan, d_map.get(), d_half.get())); + half_u_to_half_w<<((hkl.data.size() + THREADS - 1) / THREADS), THREADS, 0, stream>>>( + d_half, map.nu, map.nv, map.nw, norm, d_out); + cuda_err(cudaGetLastError()); + static_assert(sizeof(std::complex) == sizeof(float2)); + cuda_err(cudaMemcpyAsync(hkl.data.data(), d_out.get(), hkl.data.size() * sizeof(float2), + cudaMemcpyDeviceToHost, stream)); + cuda_err(cudaStreamSynchronize(stream)); + } catch (...) { + cufftDestroy(plan); + throw; + } + cufftDestroy(plan); + return hkl; +} + struct MolecularReplacementGPUState { CudaStream stream; CudaDevicePtr grid; diff --git a/image_analysis/structure_refinement/MolecularReplacementGPU.h b/image_analysis/structure_refinement/MolecularReplacementGPU.h index 4efb7f170..7924b5148 100644 --- a/image_analysis/structure_refinement/MolecularReplacementGPU.h +++ b/image_analysis/structure_refinement/MolecularReplacementGPU.h @@ -15,6 +15,13 @@ #include "MolecularReplacement.h" +// MapToFPhi (ModelFFT.h) with cuFFT on the calling thread's current device: the same transform, the same +// half-l layout and conventions, to float rounding. For the molecular transform of a large model, whose box is +// hundreds of millions of points. MapToFPhiGPUFits: whether its device buffers take at most a quarter of the +// card's total memory - decided on the total, so the same input always takes the same path. +bool MapToFPhiGPUFits(const gemmi::Grid &map); +gemmi::FPhiGrid MapToFPhiGPU(const gemmi::Grid &map); + struct MolecularReplacementGPUState; struct MolecularReplacementGPUGroup; diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index a71e3a743..96a41fafb 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2587,6 +2587,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, 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, KeyInt("AF_PAE_DOMAINS", mp.pae_domains)); + Add(s, KeyText("MODEL_PRIOR_SEARCH_MODEL", mp.search_model)); 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)); @@ -2612,22 +2614,23 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyText("MODEL_PRIOR_FILE_PREFIX", result.model_prior_file_prefix)); Add(s, KeyEnum("MODEL_PRIOR_SEARCH_DEVICE", mp.gpu_search ? "GPU" : "CPU", {"GPU", "CPU"})); Add(s, KeyReal("MODEL_PRIOR_SECONDS", mp.t_model + mp.t_transform + mp.t_rotation + mp.t_translation - + mp.t_rescore + mp.t_null + mp.hand_check.seconds, "{:.2f}")); + + mp.t_rescore + mp.t_null + mp.t_domains + mp.hand_check.seconds, "{:.2f}")); Add(s, Blank()); ReportEntry e; e.kind = ReportEntry::Kind::Table; - e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "CLASH", - "R_WORK", "R_FREE", "NULL_SIGMA", "VERDICT"}; - e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z CLASH R_WORK R_FREE NULL_SIGMA VERDICT\n" - " ------------------- ----------- ----- ------ ------ ----- ----- ----- ------ ------ ---------- ----------"; + e.table.columns = {"MERGE", "SPACE_GROUP", "VM", "COPIES", "COPIES_LIKELY", "RF_Z", "TF_Z", "LLG_Z", "CLASH", + "R_WORK", "R_FREE", "R_FREE_TWIN", "NULL_SIGMA", "VERDICT"}; + e.table.text_header = " MERGE SPACE_GROUP VM COPIES LIKELY RF_Z TF_Z LLG_Z CLASH R_WORK R_FREE R_FREE_TWIN 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} {:6d} {:6d} {:5.1f} {:5.1f} {:4.1f}% {:>6s} {:>6s} {:>10s} {}", + " {:<19s} {:<11s} {:5.2f} {:6d} {:6d} {:5.1f} {:5.1f} {:>5s} {:4.1f}% {:>6s} {:>6s} {:>11s} {:>10s} {}", merge, c.space_group.xhm(), c.vm, c.copies, c.copies_likely, c.rf_z, c.tf_z, + std::isfinite(c.llg_z) ? fmt::format("{:.1f}", c.llg_z) : std::string("-"), 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("-"), + r(c.r_free), r(c.r_free_twin), std::isfinite(c.null_sigma) ? fmt::format("{:+.1f}", c.null_sigma) : std::string("-"), c.verdict)); std::vector row; row.push_back(KeyText("", merge).value); @@ -2637,9 +2640,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, row.push_back(KeyInt("", c.copies_likely).value); row.push_back(KeyReal("", c.rf_z, "{:.1f}").value); row.push_back(KeyReal("", c.tf_z, "{:.1f}").value); + row.push_back(KeyReal("", c.llg_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.r_free_twin, "{:.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)); @@ -2651,16 +2656,24 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " cannot tell apart (rotation function on a 6 deg sampling, translation function, {:.1f}-{:.1f} A),\n" " and the best placements taken through the rigid body; then further copies, one at a time, up to\n" " the largest copy number with a Matthews probability of 1% (COPIES_LIKELY: the most probable),\n" - " each with the translation function beside the copies already placed, kept where its TF_Z is 8\n" - " or more, it packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES, to\n" - " SCORE_D_MIN; TF_Z is the highest of the copies'. PLACED needs TF_Z >= 8, at most 5% of the\n" + " each with the translation function beside the copies already placed, its peaks ranked by the\n" + " likelihood gain of the copy (sigma_A from the pLDDT's coordinate error, its scale refined), kept\n" + " where the copy stands out, packs and R-free drops by 0.01. R_WORK and R_FREE are of all COPIES,\n" + " to SCORE_D_MIN; TF_Z is the highest of the copies'; LLG_Z the highest Z-score of a copy's refined\n" + " likelihood gain over the same copy at random positions, where its TF_Z was below 8. A copy stands\n" + " out where TF_Z or LLG_Z is 8 or more. PLACED needs a copy that stands out, at most 5% of the\n" " C-alpha atoms clashing with another copy or a symmetry mate, one copy fitting the cell\n" " (MATTHEWS: Vm below 1.5) and, for the candidate the choice falls on, an R-work at least 3 sd\n" " better than the same copies refitted from random orientations (NULL where it is not; the next\n" " candidate is then taken). The model chooses within a merge where its placed group leads every\n" - " other that packs by 0.03 in R-free, and the twin hypothesis' merge - searched only then - where\n" - " no group of the answer's merge can hold it at all. A group one copy cannot fit is checked for\n" - " packing only, not refined.\n" + " other that packs by 0.03 in R-free, and the twin hypothesis' merge where no group of the\n" + " answer's merge can hold it at all, or where, on that merge and against its data as twinned by\n" + " its law, the model searched for in the subgroup leads the chosen placement read as a crystal\n" + " there by 0.03 in R-free (R_FREE_TWIN). A group one copy cannot fit is checked for\n" + " packing only, not refined. Where the AlphaFold DB's predicted aligned error sits beside the\n" + " model, the prediction is cut into rigid domains (AF_PAE_DOMAINS), residues in none left out;\n" + " where the whole model is placed nowhere, each domain is searched by itself, and the domains\n" + " placed have to agree (MODEL_PRIOR_SEARCH_MODEL: what the decision was made with).\n" " The answer's files take the enantiomorph from it - only where model validation of the placed\n" " model accepts it as it accepts a --model's: a fit over its null and an anomalous map at its\n" " atoms that is not inverted (ENANTIOMORPH_VETOED otherwise) - and nothing else; another group it\n" @@ -2669,8 +2682,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, 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}, hand check {:.3f}.", mp.t_model, mp.t_transform, mp.t_rotation, - mp.t_translation, mp.t_rescore, mp.t_null, mp.hand_check.seconds), true)); + " R-free {:.2f}, null {:.2f}, PAE domains {:.2f}, hand check {:.3f}.", mp.t_model, mp.t_transform, + mp.t_rotation, mp.t_translation, mp.t_rescore, mp.t_null, mp.t_domains, mp.hand_check.seconds), true)); } doc.sections.push_back(std::move(s)); } diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 07c9095b4..ce840cce9 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -2874,16 +2874,22 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // too, with the law - so a refinement program can take the twin hypothesis without a re-run. // The answer itself is not touched by any of this. const bool pre_twin_like = result.pre_promotion_twinning && result.pre_promotion_twinning->twinning_suspected; + // A predicted model can test the twin hypothesis itself (--af-model's twin comparison), so with one the + // merge in the indexed lattice's group is made wherever a metric re-ask promoted from it. + const bool model_tests_twin = !config_.af_model_path.empty() && twin_metric_subgroup.has_value(); if (!geometry_prepass && result.space_group && sg_search.point_group_order > 1 - && (twin_overruled_subgroup || pre_twin_like)) { + && (twin_overruled_subgroup || pre_twin_like || model_tests_twin)) { twin_choice = ChooseTwinSubgroup(*result.space_group, sg_search, twin_overruled_subgroup ? twin_overruled_subgroup : twin_metric_subgroup); twin_choice_group_number = result.space_group->number; twin_choice_reason = twin_overruled_subgroup ? fmt::format("a twin gate refused it ({}) and was overruled", twin_overruled_reason) - : fmt::format("the merge read twin-like before the space-group search ({})", - TwinningVerdictLine(*result.pre_promotion_twinning)); + : pre_twin_like + ? fmt::format("the merge read twin-like before the space-group search ({})", + TwinningVerdictLine(*result.pre_promotion_twinning)) + : std::string("it was promoted from the indexed lattice's group, which the predicted model " + "can test"); if (twin_choice) logger.Info("Twin hypothesis: {} was adopted, but {} - the merge in {} is written too, " "with twin law {}", result.space_group->xhm(), twin_choice_reason, @@ -3992,7 +3998,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { 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}); + merges.push_back({"TWIN_SUBGROUP_MERGE", &twin_file->merged, twin_file->cell, twin_file->space_group, + twin_file->law}); MolecularReplacementResult mr = RunMolecularReplacement(config_.af_model_path, merges, logger, static_cast(config_.nthreads)); if (!mr.ok) diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index fb6fc0b02..994508557 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -202,7 +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