From c6912e2a066008fc15a020ab6af12ee8030c008e Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 12 Oct 2026 16:49:02 +0200 Subject: [PATCH] rugnux --af-model: tNCS half-vector moves checked, the merge's own group not ruled out beside tNCS, P1's first copy judged by its orientation (WIP, validation pending) Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3 --- docs/ACKNOWLEDGEMENT.md | 7 +- .../MolecularReplacement.cpp | 172 +++++++++++++++++- 2 files changed, 168 insertions(+), 11 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index d6264b1b6..856cba036 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -481,7 +481,12 @@ boundaries of molecular replacement with maximum likelihood" (2001), Acta Cryst. 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 +[doi:10.1107/S0907444905001617](https://doi.org/10.1107/S0907444905001617). A model whose expected +log-likelihood gain - the sum of sigma_A^4 / 2 over the reflections - is too small for a placement to stand +out is not searched for: A. J. McCoy, R. D. Oeffner, A. G. Wrobel, J. R. M. Ojala, K. Tryggvason, +B. Lohkamp and R. J. Read, "Ab initio solution of macromolecular crystal structures without direct +methods" (2017), Proc. Natl. Acad. Sci. USA 114, 3637-3641 +[doi:10.1073/pnas.1701640114](https://doi.org/10.1073/pnas.1701640114). 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 diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 9c11693a9..450225a19 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -102,6 +102,14 @@ constexpr double NULL_SIGMA_MIN = 4.51; // A model whose expected log-likelihood gain (ExpectedLogLikelihoodGain) is below this is not searched for. The // lowest of a model the search did place, over 193 survey sets: 26; every PAE domain below 24 placed nowhere. constexpr double ELLG_MIN = 10.0; +// tNCS by half a lattice vector v: copies v apart make the reflections with h.v half-integral weak, and a +// placement moved by (v + o) / 2, o an origin shift of the group, fits the strong ones as well as the right one +// does; only the weak ones tell them apart. A placement is compared with every such move of itself on those +// reflections (R_tncs_weak), and one of them that fits them better by TNCS_SHIFT_MARGIN is taken instead; where the +// two best are closer than that, the placement is not decided on. A vector is taken as half-integral within +// TNCS_HALF_TOLERANCE. +constexpr double TNCS_SHIFT_MARGIN = 0.01; +constexpr double TNCS_HALF_TOLERANCE = 0.05; double seconds_since(std::chrono::steady_clock::time_point t0) { return std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); @@ -161,6 +169,9 @@ struct PlacementScore { // 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; + // Over the reflections a tNCS vector v weakens - h.v half-integral - work and free together (score_placement's + // `tncs`). + double r_tncs_weak = NAN; }; gemmi::AsuData> working_amplitudes(const std::vector &merged, @@ -205,37 +216,54 @@ void model_factors(const gemmi::Structure &st, double d_min, gemmi::AsuData> &fmodel, const gemmi::Scaling &scaling, - const gemmi::SpaceGroup &sg, const std::vector &merged, double d_min) { + const gemmi::SpaceGroup &sg, const std::vector &merged, double d_min, + const gemmi::Fractional *tncs = nullptr) { const gemmi::GroupOps gops = sg.operations(); const gemmi::ReciprocalAsu asu(&sg); - std::unordered_map> obs; + struct Obs { + double f; + bool free; + bool weak; // h.v half-integral + }; + std::unordered_map obs; for (const MergedReflection &r : merged) { if (std::isnan(r.F) || r.d < d_min) continue; gemmi::Miller h{{r.h, r.k, r.l}}; if (!asu.is_in(h)) h = asu.to_asu(h, gops).first; - obs[hkl_key(h)] = {r.F, r.rfree_flag}; + bool weak = false; + if (tncs != nullptr) { + const double hv = h[0] * tncs->x + h[1] * tncs->y + h[2] * tncs->z; + weak = std::fabs(hv - std::round(hv)) > 0.25; + } + obs[hkl_key(h)] = {r.F, r.rfree_flag, weak}; } - double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + double num_w = 0, den_w = 0, num_f = 0, den_f = 0, num_t = 0, den_t = 0; for (const auto &hv : fmodel.v) { const auto it = obs.find(hkl_key(hv.hkl)); if (it == obs.end()) continue; - const double fo = it->second.first, fc = std::abs(hv.value); - if (it->second.second) { + const double fo = it->second.f, fc = std::abs(hv.value); + if (it->second.free) { num_f += std::fabs(fo - fc); den_f += fo; } else { num_w += std::fabs(fo - fc); den_w += fo; } + if (it->second.weak) { + num_t += std::fabs(fo - fc); + den_t += fo; + } } PlacementScore out; out.k_sol = scaling.k_sol; out.b_sol = scaling.b_sol; out.r_work = den_w > 0 ? num_w / den_w : 1; out.r_free = den_f > 0 ? num_f / den_f : 1; + if (tncs != nullptr) + out.r_tncs_weak = den_t > 0 ? num_t / den_t : 1; return out; } @@ -423,7 +451,7 @@ void twin_r_factors(const gemmi::AsuData> &fmodel, const gem 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, - const gemmi::Op *law = nullptr) { + const gemmi::Op *law = nullptr, const gemmi::Fractional *tncs = nullptr) { gemmi::AsuData> fmask, fmodel; #ifdef JFJOCH_USE_CUDA ScoreEngines::Engine *gpu = engines != nullptr ? engines->get(st, sg, d_min) : nullptr; @@ -454,7 +482,7 @@ PlacementScore score_placement(const gemmi::Structure &st, const gemmi::SpaceGro FitModelScale(scaling, {.b_hi = SolventBMax(st.models[0])}, nthreads); } scaling.scale_data(fmodel, &fmask); - PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min); + PlacementScore out = r_factors(fmodel, scaling, sg, merged, d_min, tncs); if (law != nullptr) twin_r_factors(fmodel, sg, merged, d_min, *law, out); return out; @@ -651,6 +679,33 @@ bool enantiomorphs(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) { return other != nullptr && other->number == b.number; } +// The tNCS vector where it is half a lattice vector - every component within TNCS_HALF_TOLERANCE of 0 or 1/2 - +// on those halves; empty otherwise. +std::optional half_tncs_vector(const TranslationalNcs &tncs) { + if (!tncs.present) + return std::nullopt; + double v[3] = {tncs.vector.x, tncs.vector.y, tncs.vector.z}; + for (double &x : v) { + const double h = std::round(2 * x) / 2; + if (std::fabs(x - h) > TNCS_HALF_TOLERANCE) + return std::nullopt; + x = h - std::floor(h); + } + return gemmi::Fractional(v[0], v[1], v[2]); +} + +// Whether moving the origin by `s` (a multiple of 1/24 in each component) leaves `sg`'s operators as they are: +// a placement moved by s is then the same crystal. +bool is_origin_shift(const gemmi::SpaceGroup &sg, const gemmi::Fractional &s) { + gemmi::Op shift = gemmi::Op::identity(); + const double x[3] = {s.x, s.y, s.z}; + for (int a = 0; a < 3; a++) + shift.tran[a] = static_cast(std::lround(x[a] * gemmi::Op::DEN)) % gemmi::Op::DEN; + gemmi::GroupOps ops = sg.operations(); + ops.change_basis_forward(shift); + return ops.is_same_as(sg.operations()); +} + } // namespace // --- the model --- @@ -2531,8 +2586,11 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec cand.verdict = "CLASH"; // 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))) + // likelihood over random positions of its own is no measure of the placement. In P1 the first copy sits + // at the origin, which is arbitrary there and which no function of its position can single out: one copy + // is judged by its orientation alone, against the null. + else if (!(cand.tf_z >= TF_Z_PLACED || (cand.llg_z >= TF_Z_PLACED && cand.copies >= cand.copies_likely) + || (gops.order() == 1 && cand.copies == 1))) cand.verdict = "NOT_PLACED"; else cand.verdict = "PLACED"; @@ -2630,6 +2688,79 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec // by side, each taken where the choice comes to it. std::vector nulls(result.candidates.size(), NAN); bool first = true; + // With tNCS by half a lattice vector, candidate i against its moves by half that vector (TNCS_SHIFT_MARGIN): + // the better taken into it, and whether the two best are too close to tell apart. + // The tNCS vector is the answer's merge's, in its cell. + const std::optional half = merge == 0 ? half_tncs_vector(result.tncs) : std::nullopt; + std::vector shift_checked(result.candidates.size(), 0); + const auto check_shifts = [&](int i) -> bool { + MRCandidate &c = result.candidates[i]; + const gemmi::Fractional v = *half; + const std::vector &merged = *merges[merge].merged; + const double d_min = result.d_min_score; + struct Option { + gemmi::Fractional shift; + std::shared_ptr placed; + PlacementScore score; + double clash = 0; + }; + std::vector