From 613b3eda2cdfa282b701f788d40a51be8af0c030 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 27 Oct 2026 16:09:40 +0100 Subject: [PATCH 1/2] RESULT-CHANGING: Refinement: copies of a molecule found by sequence alignment and restrained locally, decided by R-free Polymer chains whose sequences align with more than 15 residues at 80% identity (REFMAC5's criterion) are copies; atoms of the same name in aligned residues of the same name correspond. Each copy is paired with the next of its molecule and restrained on its interatomic distances below 4.2 A (bonds and angles excepted) against the other copy's, under a Geman-McClure loss that saturates past 0.5 A, so regions that really differ release themselves and rigid domain motion costs nothing. No superposition, no B-factor NCS, waters and ligands free. Use and strength are decided late, by the weight trials: every rung of the first ladder runs without NCS and then with it at sigma 0.1 A, and every trajectory step also tries it switched on / halved / doubled (off beyond 0.4 A); listed after its twin, NCS is kept only where it lowers R-free by more than the tie. A model without copies runs exactly the trials it ran before (bit-identical refined model on a single-copy control). New kind RestraintKind::NcsDistance evaluated by the CPU and the GPU restraint kernels; REFINE_NCS_ report keys where copies exist. On MD-perturbed starts of the refine-battery sets (refine mode, depositor's free set): 3.1-3.5 A sets R-free -0.012 to -0.073 (one start +0.001), 1.95 A -0.001 to -0.004. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/ACKNOWLEDGEMENT.md | 16 ++ docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 31 ++- docs/RUGNUX_REPORT.md | 8 + .../structure_refinement/CMakeLists.txt | 1 + .../structure_refinement/Refinement.cpp | 115 +++++--- .../structure_refinement/Refinement.h | 14 + .../structure_refinement/RefinementGPU.cu | 5 + .../structure_refinement/RefinementMath.h | 26 +- .../structure_refinement/RefinementNcs.cpp | 249 ++++++++++++++++++ .../structure_refinement/RefinementNcs.h | 48 ++++ .../RefinementRestraints.cpp | 8 + rugnux/ResultReport.cpp | 10 + tests/RefinementTest.cpp | 106 ++++++++ 13 files changed, 604 insertions(+), 33 deletions(-) create mode 100644 image_analysis/structure_refinement/RefinementNcs.cpp create mode 100644 image_analysis/structure_refinement/RefinementNcs.h diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index ac8c70426..10332fba5 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -596,6 +596,22 @@ predicted distances around metal sites in proteins" (2006), Acta Cryst. D62, 678 nucleic acids, F. Leonarski, L. D'Ascenzo and P. Auffinger, "Mg2+ ions: do they bind to nucleobase nitrogens?" (2017), Nucleic Acids Res. 45, 987-1004 [doi:10.1093/nar/gkw1175](https://doi.org/10.1093/nar/gkw1175). +**NCS restraints** — copies of a molecule in the model are found by sequence alignment with REFMAC5's criterion +(more than 15 residues aligned, 80% identity; Murshudov et al. 2011, above) and restrained locally, on pairs of +corresponding interatomic distances under the Geman-McClure robust loss, as REFMAC5's local NCS restraints are, with +the distance range of R. A. Nicholls, F. Long and G. N. Murshudov, "Low-resolution refinement tools in REFMAC5" +(2012), Acta Cryst. D68, 404-417 [doi:10.1107/S090744491105606X](https://doi.org/10.1107/S090744491105606X), and +the bonds and angles left out as they leave them out. The case for local restraints over a superposition, and for +deciding them by R-free, is G. J. Kleywegt's, "Use of non-crystallographic symmetry in protein structure refinement" +(1996), Acta Cryst. D52, 842-857 [doi:10.1107/S0907444995016477](https://doi.org/10.1107/S0907444995016477); I. Usón, +E. Pohl, T. R. Schneider, Z. Dauter, A. Schmidt, H.-J. Fritz and G. M. Sheldrick, "1.7 Å structure of the stabilized +REIv mutant T39K. Application of local NCS restraints" (1999), Acta Cryst. D55, 1158-1167 +[doi:10.1107/S0907444999003972](https://doi.org/10.1107/S0907444999003972); and J. J. Headd, N. Echols, +P. V. Afonine, N. W. Moriarty, R. J. Gildea and P. D. Adams, "Flexible torsion-angle noncrystallographic symmetry +restraints for improved macromolecular structure refinement" (2014), Acta Cryst. D70, 1346-1356 +[doi:10.1107/S1399004714003277](https://doi.org/10.1107/S1399004714003277), whose comparison of global, torsion and +no NCS is why the copies' B factors are refined independently here. + **Unmodelled density** — the blobs the refinement lists are the connected regions of the difference map above a cutoff, found by flood fill as Coot's "unmodelled blobs" finds them, in GEMMI's implementation (`find_blobs_by_flood_fill`): P. Emsley, B. Lohkamp, W. G. Scott and K. Cowtan, "Features and development diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index 0c9d74b4f..d7168120d 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -197,6 +197,26 @@ common atom, 0.2 and 0.3 for the side chains and the rest, and 0.5 for atoms wit among them, symmetry mates included) that no bond or angle relates. No atom's B goes below 2 Ų or above the ceiling of §15.1. +**NCS restraints.** Where the model holds several copies of one molecule, they are found from the model alone: two +polymer chains of one kind (protein, or nucleic acid) are copies where their sequences align with more than 15 +residues at 80% identity or more, REFMAC5's criterion. In the residues aligned to a residue of the same name, the +atoms of the same name correspond - not a hydrogen, not an alternative conformation, both refinable. Each copy is +paired with the next of its molecule (and the last with the first, where there are three or more), and each pair is +restrained **locally**, on its interatomic distances rather than on a superposition: every distance within a chain +shorter than 4.2 Å in either copy, except those a bond or an angle of the dictionary already fixes, against the same +distance in the other copy, with the robust loss of Geman and McClure, + +$$ +E = rac{(\Delta/\sigma)^2}{1 + (\Delta/0.5\, ext{Å})^2},\qquad \Delta = d_\mathrm{copy\,1} - d_\mathrm{copy\,2}, +$$ + +the squared normalised difference where the two distances agree and never more than $(0.5\, ext{Å}/\sigma)^2$ where +they do not. A region in which the copies really differ - a loop in a crystal contact, a side chain in another +rotamer, a domain turned on a hinge - lets itself go as the refinement moves it, without a rule deciding in advance +which residues to exclude, and rigid motions of whole regions cost nothing, since only distances within 4.2 Å are +compared. Waters and ligands are not restrained, nor are B factors: the copies' B factors are refined independently. +Whether the restraints are used at all, and how tightly ($\sigma$), the weight trials decide (§15.7). + ### 15.6 The minimiser Positions and B factors are minimised in turn, each by an L-BFGS of rugnux's own (`RefinementMinimiser`), @@ -248,6 +268,12 @@ are: trajectory ends where it began. 3. **The final ladder**: the trajectory's weights times $\{0.5, 1, 2\}$, kept where one displaces it. +Where the model has copies (§15.5), the NCS restraints are one more weight: every rung of the first ladder is run +without them and then with them at $\sigma = 0.1$ Å, and every step of the trajectory also tries them changed - +switched on where they are off, and in turn tightened ($\sigma$ halved) and loosened ($\sigma$ doubled, off beyond +0.4 Å) where they are on. Listed after their twin without, the restraints are kept only where they lower R-free by +more than the tie. A model without copies runs exactly the trials it ran before. + $w$ never exceeds $w_0/8$. The trials of one set run side by side - on one GPU, on several, or on the CPU's threads - each from the same start and taking nothing from the one before it, so the result does not depend on how many run at once. @@ -387,7 +413,9 @@ What the refinement does not do yet - each a known gap, not a design choice: - **TLS** refinement. - **Anisotropic ADPs**: an anisotropic B is replaced by its isotropic equivalent, and only isotropic B is refined. -- **NCS restraints**: copies are refined independently. +- **Global NCS and NCS of B factors**: copies are restrained locally, on their interatomic distances (§15.5), + never by a superposition, and their B factors are refined independently; ligands and waters are not restrained + between copies. - **Twin refinement**: the target is not a twinned one, also where the twin hypothesis' merge was written. - **Ion identification**: waters are built and removed (§15.10), but no ion is placed or relabelled; the solvent audit scores the waters and ions there are. @@ -417,6 +445,7 @@ What the refinement does not do yet - each a known gap, not a design choice: - A. T. Brünger, M. Karplus & G. A. Petsko, "Crystallographic refinement by simulated annealing: application to crambin", *Acta Cryst.* **A45** (1989), 50-61 (the starting weight from the ratio of the gradients, §15.7). - P. V. Afonine et al., "Towards automated crystallographic structure refinement with phenix.refine", *Acta Cryst.* **D68** (2012), 352-367 (weights chosen among refinements by R-free within geometry bounds, and the overall B moved into the atoms, §15.3, §15.7). - A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free). +- G. J. Kleywegt, "Use of non-crystallographic symmetry in protein structure refinement", *Acta Cryst.* **D52** (1996), 842-857; I. Usón, E. Pohl, T. R. Schneider, Z. Dauter, A. Schmidt, H.-J. Fritz & G. M. Sheldrick, "1.7 Å structure of the stabilized REIv mutant T39K. Application of local NCS restraints", *Acta Cryst.* **D55** (1999), 1158-1167; R. A. Nicholls, F. Long & G. N. Murshudov, "Low-resolution refinement tools in REFMAC5", *Acta Cryst.* **D68** (2012), 404-417; J. J. Headd et al., "Flexible torsion-angle noncrystallographic symmetry restraints for improved macromolecular structure refinement", *Acta Cryst.* **D70** (2014), 1346-1356 (local NCS restraints with a robust loss, decided by R-free, and Murshudov et al. 2011 above for the copies found by alignment, §15.5, §15.7). - P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (unmodelled blobs, §15.9). - N. Naleem, A. Henning-Knechtel, S. Kirmizialtin & P. Auffinger, "Cat_Wiz: a stereochemistry-guided toolkit for locating, diagnosing, and annotating Mg2+ ions in RNA structures", *Nucleic Acids Res.* **54** (2026), gkag120; F. Leonarski, L. D'Ascenzo & P. Auffinger, *Nucleic Acids Res.* **45** (2017), 987-1004 and *RNA* **25** (2019), 173-192; F. Leonarski, A. Henning-Knechtel, S. Kirmizialtin, E. Ennifar & P. Auffinger, *Nucleic Acids Res.* **53** (2025), gkae1148 (the coordination bands, distances and exclusion zone of the solvent audit, §15.10). - P. V. Afonine et al., "Towards automated crystallographic structure refinement with phenix.refine", *Acta Cryst.* **D68** (2012), 352-367, and P. Emsley, B. Lohkamp, W. G. Scott & K. Cowtan, "Features and development of Coot", *Acta Cryst.* **D66** (2010), 486-501 (the peak cutoff, the partner and clash distances and the resolution limit of the water build, §15.10). diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 8865f17d4..efbf09e0d 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1197,6 +1197,14 @@ Where a model was placed and refined ([Refinement](RUGNUX_MODEL_AND_REFERENCE.md | `REFINE_RMSD_BONDS_START`, `REFINE_RMSD_ANGLES_START` | *dev*. The same for the model as given | | `REFINE_CLASHES` | Pairs of atoms in the refined model more than 0.4 Å closer than their contact distance (symmetry mates included) | | `REFINE_XRAY_WEIGHT` | *dev*. The X-ray weight of the positions that was kept | +| `REFINE_NCS_GROUPS` | Molecules present more than once in the model (copies found by sequence alignment, §15.5 of the refinement reference). Present only where there are any, as are the other `REFINE_NCS_` keys | +| `REFINE_NCS_COPIES` | The chains of those molecules | +| `REFINE_NCS_CHAINS` | The copies by chain name, `=` between copies of one molecule, e.g. `A=B=C D=E` | +| `REFINE_NCS_RESTRAINTS` | *dev*. The local distance restraints between pairs of copies | +| `REFINE_NCS_SIGMA` | The NCS restraints' sigma (Å) the weight trials kept; 0: they were not kept, and the copies were refined independently | +| `REFINE_NCS_RMSD_START` | The copies' C-alpha (P) rmsd after superposition (Å), averaged over the pairs restrained, before the restrained stage | +| `REFINE_NCS_RMSD` | The same, refined | +| `REFINE_NCS_RELEASED` | Residues of the refined model in which the copies differ: more than half of their restraints beyond the release distance (0.5 Å) | | `REFINE_ATOMS`, `REFINE_ATOMS_REFINED` | The model's atoms (hydrogens not counted) and those refined: the atoms of the residues the dictionary describes, with an occupancy above zero | | `REFINE_RESIDUES_FIXED`, `REFINE_RESIDUES_FIXED_NAMES` | The residues the embedded dictionary does not describe, held fixed but kept in the structure factors, the solvent mask and the clash terms, and their names | | `REFINE_BULK_SOLVENT_K_SOL`, `REFINE_BULK_SOLVENT_B_SOL` | *dev*. The flat bulk solvent of the refined model | diff --git a/image_analysis/structure_refinement/CMakeLists.txt b/image_analysis/structure_refinement/CMakeLists.txt index d7a5b82ef..e7bde871b 100644 --- a/image_analysis/structure_refinement/CMakeLists.txt +++ b/image_analysis/structure_refinement/CMakeLists.txt @@ -20,6 +20,7 @@ ADD_LIBRARY(JFJochStructureRefinement STATIC MonomerLibrary.cpp MonomerLibrary.h RefinementMath.h RefinementMinimiser.h RefinementRestraints.cpp RefinementRestraints.h + RefinementNcs.cpp RefinementNcs.h RefinementXray.cpp RefinementXray.h Refinement.cpp Refinement.h DensityBlobs.cpp DensityBlobs.h diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp index ed119cd3f..ac81328b5 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -34,6 +34,7 @@ #include "ModelScaling.h" #include "RefinementMath.h" #include "RefinementMinimiser.h" +#include "RefinementNcs.h" #include "RefinementXray.h" #ifdef JFJOCH_USE_CUDA #include "../../common/CUDAWrapper.h" @@ -58,6 +59,8 @@ constexpr double R_FREE_TIE = 0.001; // relative; the most a position (A) or a B (A^2) moves in the first trial of a line search. constexpr double MINIMISER_F_TOLERANCE = 1e-7; constexpr double MAX_STEP_XYZ = 0.3, MAX_STEP_B = 10.0; +// The loosest NCS restraint a trajectory step tries (A); looser than that, the copies are refined independently. +constexpr double NCS_SIGMA_MAX = 0.4; // The figure of merit of the likelihood, , at X = 2 Fo D Fc / var: I1(X)/I0(X) for an acentric // reflection, tanh(X / 2) for a centric one. @@ -986,12 +989,15 @@ void PlaceRiders(Context &ctx, RefinementTopology &topo) { } } -// The geometry and the repulsions of the model as it is. -RestraintSet RestraintsNow(Context &ctx, const RefinementTopology &topo) { +// The geometry and the repulsions of the model as it is, and where `ncs_sigma` is not 0 the NCS restraints at it. +RestraintSet RestraintsNow(Context &ctx, const RefinementTopology &topo, const NcsCopies *ncs = nullptr, + double ncs_sigma = 0) { const RefinementStage stage("repulsions"); RestraintSet set = topo.geometry; set.Append(topo.metal_bonds); set.Append(BuildRepulsions(topo, ctx.Model(), ctx.cell, *ctx.sg, REPULSION_MARGIN, ctx.nthreads)); + if (ncs != nullptr && ncs_sigma > 0) + set.Append(NcsRestraintsAt(*ncs, ncs_sigma)); set.Finalize(ctx.atoms.size()); return set; } @@ -1093,6 +1099,12 @@ std::vector OwnSlotDevices(const Context &main, int wanted, const Refinemen return std::vector(std::max(beside, 0), -1); } +// The weights of one trial: the positions' X-ray weight, the B factors' (0: held) and the NCS restraints' sigma (0: +// none). +struct TrialWeights { + double w = 0, adp = 0, ncs = 0; +}; + // One weight trial: `cycles` macro-cycles from pos0 and b0, of the positions at X-ray weight w and, where adp_w is // not 0, the B factors at X-ray weight adp_w. `fit`: what was fitted at pos0 and b0 already, null to fit it here. struct TrialOutcome { @@ -1106,7 +1118,7 @@ struct TrialOutcome { TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const std::vector &pos0, const std::vector &b0, const Context::Fit *fit, double w, double adp_w, int cycles, - const RefinementSettings &settings) { + const RefinementSettings &settings, const NcsCopies *ncs, double ncs_sigma) { Context &ctx = *slot.ctx; RefinementTopology &topo = *slot.topo; const RefinementStage stage(adp_w > 0 ? "rung.b" : "rung.xyz"); @@ -1129,7 +1141,7 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s if (adp_w > 0) ctx.MoveOverallB(); { - const RestraintSet set = RestraintsNow(ctx, topo); + const RestraintSet set = RestraintsNow(ctx, topo, ncs, ncs_sigma); const RefinementStage s1("set_restraints"); ctx.engine->SetRestraints(0, set); } @@ -1155,6 +1167,7 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s RefinementWeightTrial &t = out.trial; t.weight = w; t.adp_weight = adp_w; + t.ncs_sigma = ncs_sigma; const RFactors rf = ctx.Refit(); t.r_work = rf.r_work; t.r_free = rf.r_free; @@ -1165,9 +1178,14 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s &topo.metal_bonds); } t.within_bounds = t.geometry.rmsd_bonds <= MAX_RMSD_BONDS && t.geometry.rmsd_angles <= MAX_RMSD_ANGLES; - out.line = fmt::format("X-ray weight {:.3g}, B {}: R-work {:.4f}, R-free {:.4f}, rmsd bonds {:.4f} A, " + if (ncs != nullptr && !ncs->Empty()) + t.ncs_rmsd = NcsRmsd(*ncs, ctx.XyzB()); + out.line = fmt::format("X-ray weight {:.3g}, B {}{}: R-work {:.4f}, R-free {:.4f}, rmsd bonds {:.4f} A, " "angles {:.2f} deg{}{}", w, adp_w > 0 ? fmt::format("refined at weight {:.3g}", adp_w) : std::string("held"), + ncs == nullptr || ncs->Empty() ? std::string() + : ncs_sigma > 0 ? fmt::format(", NCS at {:.3g} A (copies {:.3f} A apart)", ncs_sigma, t.ncs_rmsd) + : fmt::format(", no NCS (copies {:.3f} A apart)", t.ncs_rmsd), rf.r_work, rf.r_free, t.geometry.rmsd_bonds, t.geometry.rmsd_angles, t.geometry.metal_bonds > 0 ? fmt::format(", metal bonds {:.3f} A", t.geometry.rmsd_metal) : std::string(), @@ -1183,10 +1201,10 @@ TrialOutcome RunTrial(TrialSlot slot, const std::vector &refinable, const s // its own (on its GPU), returned in the order given. Every trial starts from pos0 and b0 and takes nothing from the // one before it in its slot, so the outcome is the same for any number of slots. std::vector RunTrials(const std::vector &trial_slots, - const std::vector> &trials, + const std::vector &trials, const std::vector &refinable, const std::vector &pos0, const std::vector &b0, const Context::Fit *fit, int cycles, - const RefinementSettings &settings) { + const RefinementSettings &settings, const NcsCopies *ncs) { std::vector out(trials.size()); std::vector errors(trial_slots.size()); // the threads shared among the slots that run (the outcome does not depend on how many each has) @@ -1198,7 +1216,7 @@ std::vector RunTrials(const std::vector &trial_slots, std::vector order(trials.size()); std::iota(order.begin(), order.end(), 0); std::stable_sort(order.begin(), order.end(), - [&](size_t a, size_t b) { return (trials[a].second > 0) > (trials[b].second > 0); }); + [&](size_t a, size_t b) { return (trials[a].adp > 0) > (trials[b].adp > 0); }); auto run_slot = [&](size_t s) { try { #ifdef JFJOCH_USE_CUDA @@ -1206,8 +1224,8 @@ std::vector RunTrials(const std::vector &trial_slots, #endif for (size_t j = s; j < trials.size(); j += trial_slots.size()) if (const size_t i = order[j]; !trial_slots[s].ctx->Cancelled()) - out[i] = RunTrial(trial_slots[s], refinable, pos0, b0, fit, trials[i].first, trials[i].second, - cycles, settings); + out[i] = RunTrial(trial_slots[s], refinable, pos0, b0, fit, trials[i].w, trials[i].adp, cycles, + settings, ncs, trials[i].ncs); } catch (...) { errors[s] = std::current_exception(); } @@ -1285,6 +1303,12 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vectorOnGPU(); @@ -1348,7 +1372,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector(std::max({settings.weight_ladder.size() * settings.adp_ladder.size(), + int wanted = static_cast(std::max({settings.weight_ladder.size() * settings.adp_ladder.size() * + (ncs_p != nullptr ? 2 : 1), settings.final_ladder.size(), size_t{3}})); if (settings.max_trials_at_once > 0) wanted = std::min(wanted, static_cast(settings.max_trials_at_once)); @@ -1404,6 +1429,17 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector(ncs.distances.records.size()); + result.ncs_chains = ncs.description; + result.ncs_rmsd_start = NcsRmsd(ncs, ctx.XyzB()); + logger.Info("Refinement: NCS, {} molecule{} in {} copies ({}), {} local distance restraints; copies {:.3f} A " + "apart (C-alpha/P rmsd after superposition)", ncs.groups, ncs.groups == 1 ? "" : "s", + ncs.copies, ncs.description, result.ncs_restraints, result.ncs_rmsd_start); + } + ctx.engine->SetRestraints(1, topo.adp); { const RestraintSet set = RestraintsNow(ctx, topo); @@ -1459,13 +1495,17 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector> ladder; + // held and refined at their own balanced weight; where there are copies, each without NCS restraints and + // then with them. + std::vector ladder; for (double rung : settings.weight_ladder) - for (double adp : settings.adp_ladder) - ladder.emplace_back(w0 * rung, adp_w0 * adp); + for (double adp : settings.adp_ladder) { + ladder.push_back({w0 * rung, adp_w0 * adp, 0}); + if (ncs_p != nullptr) + ladder.push_back({w0 * rung, adp_w0 * adp, settings.ncs_sigma}); + } std::vector outcomes = RunTrials(trial_slots, ladder, refinable, pos0, b0, nullptr, - settings.macro_cycles, settings); + settings.macro_cycles, settings, ncs_p); int k = choose(outcomes, "ladder"); TrialOutcome state; if (k >= 0) @@ -1474,7 +1514,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector from_best; bool state_is_best = true; for (int step = 0; step < settings.max_steps && state.ran && !ctx.Cancelled(); step++) { - const double w = state.trial.weight, adp = state.trial.adp_weight; - std::vector> trials = {{w, adp}}; + const double w = state.trial.weight, adp = state.trial.adp_weight, ncs_s = state.trial.ncs_sigma; + std::vector trials = {{w, adp, ncs_s}}; if (step % 2 == 1 && adp > 0) { - trials.emplace_back(w, adp / 4); - trials.emplace_back(w, 4 * adp); + trials.push_back({w, adp / 4, ncs_s}); + trials.push_back({w, 4 * adp, ncs_s}); } else { - trials.emplace_back(w / 2, adp); + trials.push_back({w / 2, adp, ncs_s}); if (2 * w <= w_max) - trials.emplace_back(2 * w, adp); + trials.push_back({2 * w, adp, ncs_s}); + } + if (ncs_p != nullptr) { + const double changed = ncs_s == 0 ? settings.ncs_sigma : step % 2 == 0 ? ncs_s / 2 : 2 * ncs_s; + trials.push_back({w, adp, changed <= NCS_SIGMA_MAX ? changed : 0}); } outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, &state.fit, settings.macro_cycles, - settings); + settings, ncs_p); if (state_is_best) from_best = outcomes; k = choose(outcomes, fmt::format("step {}", step + 1).c_str()); @@ -1517,22 +1562,23 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector> trials; + const double w = state.trial.weight, adp = state.trial.adp_weight, ncs_s = state.trial.ncs_sigma; + std::vector trials; for (double f : settings.final_ladder) if (w * f <= w_max) - trials.emplace_back(w * f, adp); + trials.push_back({w * f, adp, ncs_s}); outcomes.assign(trials.size(), TrialOutcome()); - std::vector> missing; + std::vector missing; for (size_t i = 0; i < trials.size(); i++) { for (const TrialOutcome &o : from_best) - if (o.ran && o.trial.weight == trials[i].first && o.trial.adp_weight == trials[i].second) + if (o.ran && o.trial.weight == trials[i].w && o.trial.adp_weight == trials[i].adp && + o.trial.ncs_sigma == trials[i].ncs) outcomes[i] = o; if (!outcomes[i].ran) missing.push_back(trials[i]); } std::vector ran = RunTrials(trial_slots, missing, refinable, state.pos, state.b, &state.fit, - settings.macro_cycles, settings); + settings.macro_cycles, settings, ncs_p); for (size_t i = 0, m = 0; i < trials.size(); i++) if (!outcomes[i].ran) outcomes[i] = std::move(ran[m++]); @@ -1555,6 +1601,15 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector 0 ? fmt::format("restrained at {:.3g} A", result.ncs_sigma) + : std::string("not kept, the copies refined independently"), + result.ncs_rmsd, result.ncs_rmsd_start, result.ncs_released); + } result.r_work = state.trial.r_work; result.r_free = state.trial.r_free; result.geometry = state.trial.geometry; diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h index cefde1799..3615a85b9 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -59,6 +59,10 @@ struct RefinementSettings { std::vector final_ladder = {0.5, 1.0, 2.0}; double max_weight = 1.0 / 8; // the positions' X-ray weight at most, times the balanced one bool riding_hydrogens = true; // the dictionary's hydrogens added, riding on their parents + // Where the model has copies of a molecule (RefinementNcs.h), the trials try them restrained to each other, at + // NCS restraint sigmas from ncs_sigma, beside the trials without; R-free decides. + bool ncs = true; + double ncs_sigma = 0.1; // A // Waters added at the difference map's peaks and removed where the maps no longer hold them, between refinements // (SolventBuild.h), to SOLVENT_BUILD_D_MIN. bool solvent = true; @@ -110,6 +114,8 @@ struct RefinementScore { struct RefinementWeightTrial { double weight = 0; // the X-ray term's weight double adp_weight = 0; // the B refinement's X-ray weight; 0: the B factors were held + double ncs_sigma = 0; // the NCS restraints' sigma (A); 0: none + double ncs_rmsd = NAN; // the copies' C-alpha (P) rmsd after superposition, where there are copies double r_work = 1, r_free = 1; GeometryStatistics geometry; bool within_bounds = false; // the geometry bounds the choice is made inside @@ -138,6 +144,14 @@ struct RefinementResult { double weight = 0; double adp_weight = 0; // 0: the B factors kept as placed std::vector trials; + // Non-crystallographic symmetry: the molecules present more than once, their chains, the restraints between the + // copies, and the sigma chosen (0: the copies refined independently); the copies' C-alpha (P) rmsd after + // superposition, averaged over the pairs restrained, as given and refined; the residues whose restraints let go. + int ncs_groups = 0, ncs_copies = 0, ncs_restraints = 0; + std::string ncs_chains; + double ncs_sigma = 0; + double ncs_rmsd_start = NAN, ncs_rmsd = NAN; + int ncs_released = 0; GeometryStatistics geometry_start, geometry; int atoms = 0, refined_atoms = 0; diff --git a/image_analysis/structure_refinement/RefinementGPU.cu b/image_analysis/structure_refinement/RefinementGPU.cu index 90d9672ec..097e60a0f 100644 --- a/image_analysis/structure_refinement/RefinementGPU.cu +++ b/image_analysis/structure_refinement/RefinementGPU.cu @@ -267,6 +267,11 @@ __global__ void restraint_kernel(const RestraintRecord *records, int n, const in xyz(xyzb, r.atom[3]), r.ideal, r.sigma, g); na = 4; break; + case RestraintKind::NcsDistance: + e = refine_math::NcsDistanceTerm(xyz(xyzb, r.atom[0]), xyz(xyzb, r.atom[1]), xyz(xyzb, r.atom[2]), + xyz(xyzb, r.atom[3]), r.sigma, r.ideal, g); + na = 4; + break; case RestraintKind::Repulsion: { const V3f q = xyz(xyzb, r.atom[1]); const V3f image{r.rot[0] * q.x + r.rot[1] * q.y + r.rot[2] * q.z + r.tran[0], diff --git a/image_analysis/structure_refinement/RefinementMath.h b/image_analysis/structure_refinement/RefinementMath.h index af87c6bc6..3d0240de3 100644 --- a/image_analysis/structure_refinement/RefinementMath.h +++ b/image_analysis/structure_refinement/RefinementMath.h @@ -17,13 +17,14 @@ #endif // What a restraint is, for the evaluators on either side (RefinementRestraints.h builds them). -enum class RestraintKind : int { Bond = 0, Angle, Torsion, Chiral, Plane, Repulsion, BSimilarity }; +enum class RestraintKind : int { Bond = 0, Angle, Torsion, Chiral, Plane, Repulsion, BSimilarity, NcsDistance }; // One restraint. E = ((value - ideal) / sigma)^2, summed over a set (BSimilarity: sigma is kappa, BSimilarityTerm). struct RestraintRecord { int kind = 0; int atom[4] = {0, 0, 0, 0}; // Plane: atom[0] the first of its atoms in plane_atoms, atom[1] their number - float ideal = 0; // A, degrees, A^3 (chiral) or A (repulsion: the contact distance) + float ideal = 0; // A, degrees, A^3 (chiral) or A (repulsion: the contact distance; NcsDistance: the + // release distance) float sigma = 1; int period = 0; // Torsion int slot = 0; // the first of its gradient slots, one per atom @@ -199,6 +200,27 @@ REFINE_HD inline T RepulsionTerm(V3 a, V3 b, T d0, T sigma, V3 g[2]) { return r * r; } +// A distance a-b of one copy of a molecule against the same distance c-d of another copy, by the robust loss of +// Geman and McClure: E = u / (1 + (delta / release)^2), u = (delta / sigma)^2, delta = |a - b| - |c - d| - the +// squared deviation where the two distances agree, and no more than (release / sigma)^2 where they differ, so that a +// region in which the copies really differ lets itself go. +// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 (local NCS restraints on corresponding interatomic +// distances, with the Geman-McClure function) +template +REFINE_HD inline T NcsDistanceTerm(V3 a, V3 b, V3 c, V3 d, T sigma, T release, V3 g[4]) { + const V3 u = a - b, v = c - d; + const T lu = Norm(u), lv = Norm(v); + const T delta = lu - lv; + const T q = (delta / release) * (delta / release); + const T f = T(2.0) * delta / (sigma * sigma * (T(1.0) + q) * (T(1.0) + q)); + const T fu = lu > T(0) ? f / lu : T(0), fv = lv > T(0) ? f / lv : T(0); + g[0] = fu * u; + g[1] = (-fu) * u; + g[2] = (-fv) * v; + g[3] = fv * v; + return (delta / sigma) * (delta / sigma) / (T(1.0) + q); +} + // Two atoms' B alike in proportion to their size: the symmetric Kullback-Leibler divergence of two isotropic // Gaussians of variances B_i and B_j over kappa^2 (up to a factor), E = (B_i - B_j)^2 / (kappa^2 B_i B_j) - a // difference restraint whose sigma is kappa sqrt(B_i B_j). Sets dE/dB_i and dE/dB_j. diff --git a/image_analysis/structure_refinement/RefinementNcs.cpp b/image_analysis/structure_refinement/RefinementNcs.cpp new file mode 100644 index 000000000..afdb2e252 --- /dev/null +++ b/image_analysis/structure_refinement/RefinementNcs.cpp @@ -0,0 +1,249 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RefinementNcs.h" + +#include +#include +#include +#include + +#include "gemmi/polyheur.hpp" // check_polymer_type +#include "gemmi/qcp.hpp" // superpose_positions +#include "gemmi/seqalign.hpp" // align_string_sequences + +#include "RefinementMath.h" + +namespace { + +// The distances restrained: those shorter than this in either copy, the bonded and angle-related pairs left out. +// Following Nicholls, Long & Murshudov (2012) Acta Cryst. D68, 404-417 +constexpr double NCS_DISTANCE_MAX = 4.2; // A +// Where a pair's restraint lets go (NcsDistanceTerm): copies whose distance differs by this much are taken to differ. +constexpr float NCS_RELEASE = 0.5f; // A +// Copies: more than this many residues aligned, at this identity at least. +// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 +constexpr int NCS_MIN_ALIGNED = 15; +constexpr double NCS_MIN_IDENTITY = 80.0; // % + +// One polymer chain: its residues and, per residue, its atoms' indices by name. +struct NcsChain { + std::string name; + gemmi::PolymerType type = gemmi::PolymerType::Unknown; + std::vector sequence; + std::vector residue; // numbered over the model + std::vector> atoms; // name -> atom index + std::vector anchor; // CA or P, -1 for none +}; + +// Two copies' corresponding atoms (aligned residues of the same name, atoms of the same name), the residue of each in +// the first copy, and their anchors. +void Correspond(const NcsChain &a, const NcsChain &b, std::vector> &atoms, + std::vector &residue_of, std::vector> &anchors) { + const gemmi::AlignmentResult al = gemmi::align_string_sequences(a.sequence, b.sequence, {}, + gemmi::AlignmentScoring::simple()); + size_t i = 0, j = 0; + for (const gemmi::AlignmentResult::Item &item : al.cigar) + for (uint32_t k = 0; k < item.len(); k++) { + const char op = item.op(); + if (op == 'M' && a.sequence[i] == b.sequence[j]) { + for (const auto &[name, ia] : a.atoms[i]) { + auto it = b.atoms[j].find(name); + if (it != b.atoms[j].end()) { + atoms.push_back({ia, it->second}); + residue_of.push_back(a.residue[i]); + } + } + if (a.anchor[i] >= 0 && b.anchor[j] >= 0) + anchors.push_back({a.anchor[i], b.anchor[j]}); + } + if (op == 'M' || op == 'I') + i++; + if (op == 'M' || op == 'D') + j++; + } +} + +bool AreCopies(const NcsChain &a, const NcsChain &b) { + if (a.type != b.type || a.sequence.size() <= NCS_MIN_ALIGNED || b.sequence.size() <= NCS_MIN_ALIGNED) + return false; + const gemmi::AlignmentResult al = gemmi::align_string_sequences(a.sequence, b.sequence, {}, + gemmi::AlignmentScoring::simple()); + return al.match_count > NCS_MIN_ALIGNED && al.calculate_identity() >= NCS_MIN_IDENTITY; +} + +// The pairs (i, j), i < j, of `atoms` closer than `reach` at `pos`, by a grid of cells of that size. +std::set> PairsWithin(const std::vector &atoms, const std::vector &pos, + double reach) { + std::map, std::vector> cells; // cell -> positions in `atoms` + auto cell_of = [&](const gemmi::Position &p) { + return std::array{static_cast(std::floor(p.x / reach)), static_cast(std::floor(p.y / reach)), + static_cast(std::floor(p.z / reach))}; + }; + for (size_t k = 0; k < atoms.size(); k++) + cells[cell_of(pos[atoms[k]])].push_back(static_cast(k)); + std::set> out; + for (size_t k = 0; k < atoms.size(); k++) { + const std::array c = cell_of(pos[atoms[k]]); + for (int dx = -1; dx <= 1; dx++) + for (int dy = -1; dy <= 1; dy++) + for (int dz = -1; dz <= 1; dz++) { + auto it = cells.find({c[0] + dx, c[1] + dy, c[2] + dz}); + if (it == cells.end()) + continue; + for (int l : it->second) + if (atoms[l] > atoms[k] && pos[atoms[k]].dist(pos[atoms[l]]) < reach) + out.insert({atoms[k], atoms[l]}); + } + } + return out; +} + +} // namespace + +NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo) { + NcsCopies out; + std::vector pos; + std::vector chains; + int atom = 0, residue = 0; + for (const gemmi::Chain &ch : model.chains) { + NcsChain c; + c.name = ch.name; + const gemmi::ConstResidueSpan polymer = ch.get_polymer(); + c.type = gemmi::check_polymer_type(polymer); + const bool na = gemmi::is_polynucleotide(c.type); + for (const gemmi::Residue &r : ch.residues) { + const bool in_polymer = !polymer.empty() && &r >= &polymer.front() && &r <= &polymer.back(); + std::map names; + int anchor = -1; + for (const gemmi::Atom &a : r.atoms) { + pos.push_back(a.pos); + if (in_polymer && !a.is_hydrogen() && a.altloc == '\0' && a.occ > 0 && topo.refinable[atom]) { + names[a.name] = atom; + if (a.name == (na ? "P" : "CA")) + anchor = atom; + } + atom++; + } + if (in_polymer) { + c.sequence.push_back(r.name); + c.residue.push_back(residue); + c.atoms.push_back(std::move(names)); + c.anchor.push_back(anchor); + } + residue++; + } + if (gemmi::is_polypeptide(c.type) || gemmi::is_polynucleotide(c.type)) + chains.push_back(std::move(c)); + } + + // The molecules: each chain joins the first group whose first chain it is a copy of. + std::vector> groups; + for (size_t k = 0; k < chains.size(); k++) { + bool placed = false; + for (std::vector &g : groups) + if (AreCopies(chains[g[0]], chains[k])) { + g.push_back(static_cast(k)); + placed = true; + break; + } + if (!placed) + groups.push_back({static_cast(k)}); + } + + // The distances the dictionary fixes: bonded atoms and those bonded to a common atom. + std::set> fixed; + for (const RestraintRecord &r : topo.geometry.records) + if (r.kind == static_cast(RestraintKind::Bond)) + fixed.insert(std::minmax(r.atom[0], r.atom[1])); + else if (r.kind == static_cast(RestraintKind::Angle)) + fixed.insert(std::minmax(r.atom[0], r.atom[2])); + + for (const std::vector &g : groups) { + if (g.size() < 2) + continue; + out.groups++; + out.copies += static_cast(g.size()); + std::string names; + for (int k : g) + names += (names.empty() ? "" : "=") + chains[k].name; + out.description += (out.description.empty() ? "" : " ") + names; + // each copy with the next, and the last with the first where that is a pair not made already + const size_t n_pairs = g.size() == 2 ? 1 : g.size(); + for (size_t p = 0; p < n_pairs; p++) { + const NcsChain &a = chains[g[p]], &b = chains[g[(p + 1) % g.size()]]; + std::vector> corr, anchors; + std::vector residue_of; + Correspond(a, b, corr, residue_of, anchors); + out.anchors.push_back(std::move(anchors)); + std::unordered_map to_b, res_of; + std::vector atoms_a, atoms_b; + for (size_t k = 0; k < corr.size(); k++) { + to_b[corr[k][0]] = corr[k][1]; + res_of[corr[k][0]] = residue_of[k]; + atoms_a.push_back(corr[k][0]); + atoms_b.push_back(corr[k][1]); + } + std::set> pairs = PairsWithin(atoms_a, pos, NCS_DISTANCE_MAX); + std::unordered_map to_a; + for (const auto &[ia, ib] : to_b) + to_a[ib] = ia; + for (const auto &[i, j] : PairsWithin(atoms_b, pos, NCS_DISTANCE_MAX)) + pairs.insert(std::minmax(to_a[i], to_a[j])); + for (const auto &[i, j] : pairs) { + if (fixed.count({i, j}) > 0) + continue; + RestraintRecord r; + r.kind = static_cast(RestraintKind::NcsDistance); + r.atom[0] = i; + r.atom[1] = j; + r.atom[2] = to_b[i]; + r.atom[3] = to_b[j]; + r.ideal = NCS_RELEASE; + r.sigma = 1.0f; + out.distances.Add(r, 4); + out.record_residue.push_back(res_of[i]); + } + } + } + return out; +} + +RestraintSet NcsRestraintsAt(const NcsCopies &ncs, double sigma) { + RestraintSet set = ncs.distances; + for (RestraintRecord &r : set.records) + r.sigma = static_cast(sigma); + return set; +} + +double NcsRmsd(const NcsCopies &ncs, const std::vector> &xyzb) { + double sum = 0; + int n = 0; + for (const std::vector> &anchors : ncs.anchors) { + if (anchors.size() < 3) + continue; + std::vector pa, pb; + for (const std::array &p : anchors) { + pa.emplace_back(xyzb[p[0]][0], xyzb[p[0]][1], xyzb[p[0]][2]); + pb.emplace_back(xyzb[p[1]][0], xyzb[p[1]][1], xyzb[p[1]][2]); + } + sum += gemmi::superpose_positions(pa.data(), pb.data(), pa.size(), nullptr).rmsd; + n++; + } + return n > 0 ? sum / n : NAN; +} + +int NcsReleasedResidues(const NcsCopies &ncs, const std::vector> &xyzb) { + std::map> count; // residue -> restraints released, all + for (size_t k = 0; k < ncs.distances.records.size(); k++) { + const RestraintRecord &r = ncs.distances.records[k]; + const double delta = RestraintValue(ncs.distances, r, xyzb); + std::array &c = count[ncs.record_residue[k]]; + c[0] += std::fabs(delta) > r.ideal; + c[1]++; + } + int released = 0; + for (const auto &[res, c] : count) + released += 2 * c[0] > c[1]; + return released; +} diff --git a/image_analysis/structure_refinement/RefinementNcs.h b/image_analysis/structure_refinement/RefinementNcs.h new file mode 100644 index 000000000..239acb2f9 --- /dev/null +++ b/image_analysis/structure_refinement/RefinementNcs.h @@ -0,0 +1,48 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include + +#include "gemmi/model.hpp" + +#include "RefinementRestraints.h" + +// Non-crystallographic symmetry for restrained refinement (Refinement.h): the copies of one molecule in the +// asymmetric unit, found from the model alone, and local restraints that hold each copy's interatomic distances to +// the same distances in another copy. Whether they are used, and how tightly, the weight trials decide by R-free. + +// The copies of a model and the restraints between them. +struct NcsCopies { + // Per pair of copies: a distance in one, a-b, against the same distance in the other, c-d (NcsDistanceTerm), + // with sigma 1 A (NcsRestraintsAt sets it) and the release distance in `ideal`. + RestraintSet distances; + std::vector record_residue; // per record: its residue, numbered over the model in order + // Per pair of copies: their corresponding C-alpha (P for a nucleic acid) atoms, for the rmsd between them. + std::vector>> anchors; + int groups = 0; // molecules present more than once + int copies = 0; // the chains of those molecules + std::string description; // the groups as chain names, e.g. "A=B=C D=E" + bool Empty() const { return distances.records.empty(); } +}; + +// The copies of `model` (whose entities are set up, as BuildRefinementTopology leaves them): polymer chains whose +// sequences align with more than 15 residues and 80% identity (REFMAC's criterion) are copies of one molecule. Their +// residues of the same name aligned to each other correspond, and in them the atoms of the same name - neither a +// hydrogen nor an alternative conformation, all refinable (topo.refinable). Each copy is paired with the next of its +// molecule (the last with the first, where there are three or more), and each pair restrained on the distances +// within a chain shorter than NCS_DISTANCE_MAX in either copy, those the dictionary's bonds and angles fix excepted. +// Atoms are numbered in model order, as ModelPositions() lists them. +NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo); + +// The restraints at `sigma` (A). +RestraintSet NcsRestraintsAt(const NcsCopies &ncs, double sigma); + +// The rmsd of corresponding C-alpha (P) atoms of two copies after superposition, averaged over the pairs of copies. +double NcsRmsd(const NcsCopies &ncs, const std::vector> &xyzb); + +// The residues let go: those whose restraints are more than half beyond their release distance. +int NcsReleasedResidues(const NcsCopies &ncs, const std::vector> &xyzb); diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp index e4cfc3eff..00d025a9a 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.cpp +++ b/image_analysis/structure_refinement/RefinementRestraints.cpp @@ -293,6 +293,9 @@ double RestraintValue(const RestraintSet &set, const RestraintRecord &r, const s return refine_math::Norm(At(xyzb, r.atom[0]) - Image(r, At(xyzb, r.atom[1]))); case RestraintKind::BSimilarity: return xyzb[r.atom[0]][3] - xyzb[r.atom[1]][3]; + case RestraintKind::NcsDistance: + return refine_math::Norm(At(xyzb, r.atom[0]) - At(xyzb, r.atom[1])) - + refine_math::Norm(At(xyzb, r.atom[2]) - At(xyzb, r.atom[3])); case RestraintKind::Plane: break; } @@ -328,6 +331,11 @@ double EvaluateRestraints(const RestraintSet &set, const std::vector 0) { + Add(s, KeyInt("REFINE_NCS_GROUPS", rf.ncs_groups)); + Add(s, KeyInt("REFINE_NCS_COPIES", rf.ncs_copies)); + Add(s, KeyText("REFINE_NCS_CHAINS", rf.ncs_chains)); + Add(s, KeyInt("REFINE_NCS_RESTRAINTS", rf.ncs_restraints, true)); + Add(s, KeyReal("REFINE_NCS_SIGMA", rf.ncs_sigma, "{:.3g}")); + Add(s, KeyReal("REFINE_NCS_RMSD_START", rf.ncs_rmsd_start, "{:.3f}")); + Add(s, KeyReal("REFINE_NCS_RMSD", rf.ncs_rmsd, "{:.3f}")); + Add(s, KeyInt("REFINE_NCS_RELEASED", rf.ncs_released)); + } Add(s, KeyInt("REFINE_ATOMS", rf.atoms)); Add(s, KeyInt("REFINE_ATOMS_REFINED", rf.refined_atoms)); Add(s, KeyInt("REFINE_RESIDUES_FIXED", rf.fixed_residues)); diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index d9fd7a09f..1911bef09 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -11,6 +11,7 @@ #include #include +#include #include #include "../common/CUDAWrapper.h" @@ -19,6 +20,7 @@ #include "../image_analysis/structure_refinement/Refinement.h" #include "../image_analysis/structure_refinement/RefinementMath.h" #include "../image_analysis/structure_refinement/RefinementMinimiser.h" +#include "../image_analysis/structure_refinement/RefinementNcs.h" #include "../image_analysis/structure_refinement/RefinementRestraints.h" #include "../image_analysis/structure_refinement/RefinementXray.h" #ifdef JFJOCH_USE_CUDA @@ -836,3 +838,107 @@ TEST_CASE("Refinement_restrained_recovers_model", "[Refinement]") { CHECK_FALSE(r.cancelled); } } + +// Two copies of one RNA chain in the asymmetric unit, the second turned and moved: found as copies of one molecule, held +// to each other by local distance restraints whose gradient is right on the CPU and the GPU, and let go where the +// copies really differ. A single copy has none. +TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { + gemmi::Structure st = gemmi::read_structure_gz("../../tests/test_data/ligands/1ehz_psu.pdb.gz"); + st.models.resize(1); + gemmi::remove_hydrogens(st.models[0]); + st.cell = gemmi::UnitCell(120, 120, 120, 90, 90, 90); + st.spacegroup_hm = "P 1"; + st.setup_cell_images(); + { + gemmi::Structure single = st; + const RefinementTopology topo = BuildRefinementTopology(single); + CHECK(FindNcsCopies(single.models[0], topo).Empty()); + } + gemmi::Chain copy = st.models[0].chains.at(0); + copy.name = "B"; + const double a = 0.7; + const gemmi::Mat33 rot(std::cos(a), -std::sin(a), 0, std::sin(a), std::cos(a), 0, 0, 0, 1); + for (gemmi::Residue &r : copy.residues) + for (gemmi::Atom &at : r.atoms) + at.pos = gemmi::Position(rot.multiply(at.pos)) + gemmi::Position(40, 10, -5); + st.models[0].chains.push_back(copy); + st.connections.clear(); // the LINK records name chain A only + const RefinementTopology topo = BuildRefinementTopology(st); + const NcsCopies ncs = FindNcsCopies(st.models[0], topo); + REQUIRE_FALSE(ncs.Empty()); + CHECK(ncs.groups == 1); + CHECK(ncs.copies == 2); + CHECK(ncs.description == "A=B"); + CHECK(ncs.anchors.size() == 1); + CHECK(ncs.anchors[0].size() > 25); + CAPTURE(ncs.distances.records.size()); + CHECK(ncs.distances.records.size() > 1000); + + std::vector atoms = atoms_of(st.models[0]); + std::vector> xyzb; + for (gemmi::Atom *at : atoms) + xyzb.push_back({at->pos.x, at->pos.y, at->pos.z, at->b_iso}); + RestraintSet set = NcsRestraintsAt(ncs, 0.1); + set.Finalize(xyzb.size()); + // the same molecule twice: nothing to restrain + CHECK(EvaluateRestraints(set, xyzb, nullptr) == Catch::Approx(0.0).margin(1e-6)); + CHECK(NcsRmsd(ncs, xyzb) < 1e-3); + CHECK(NcsReleasedResidues(ncs, xyzb) == 0); + + // the second copy shaken: the gradient against central differences + std::mt19937 rng(7); + std::normal_distribution d(0, 0.1); + const size_t half = xyzb.size() / 2; + for (size_t i = half; i < xyzb.size(); i++) + for (int c = 0; c < 3; c++) + xyzb[i][c] += d(rng); + std::vector> grad(xyzb.size(), {0, 0, 0, 0}); + const double e = EvaluateRestraints(set, xyzb, &grad); + CHECK(e > 0); + const double h = 1e-6; + for (size_t i = 0; i < xyzb.size(); i += 17) + for (int c = 0; c < 3; c++) { + auto x = xyzb; + x[i][c] += h; + const double fp = EvaluateRestraints(set, x, nullptr); + x[i][c] -= 2 * h; + const double fm = EvaluateRestraints(set, x, nullptr); + CAPTURE(i, c); + CHECK(grad[i][c] == Catch::Approx((fp - fm) / (2 * h)).epsilon(1e-4).margin(1e-3)); + } + +#ifdef JFJOCH_USE_CUDA + if (get_gpu_count() > 0) { + const gemmi::SpaceGroup &p1 = *gemmi::find_spacegroup_by_name("P 1"); + const gemmi::UnitCell wide_cell(40, 40, 40, 90, 90, 90); + std::vector merged = asu_reflections(wide_cell, p1, 3.0); + RefinementReflections refl = PrepareRefinementReflections(merged, wide_cell, p1, 3.0); + auto gpu = MakeRefinementEngineGPU(0, wide_cell, p1, 3.0, refl, st.models[0], 100.0, SIZE_MAX); + REQUIRE(gpu); + gpu->SetRestraints(0, set); + std::vector> ggpu; + const double egpu = gpu->EvaluateRestraints(0, xyzb, &ggpu); + CHECK(egpu == Catch::Approx(e).epsilon(1e-3)); + for (size_t i = 0; i < xyzb.size(); i++) + for (int c = 0; c < 3; c++) + CHECK(ggpu[i][c] == Catch::Approx(grad[i][c]).epsilon(2e-3).margin(2e-2)); + } +#endif + + // one residue of the second copy (its 16th) distorted, each atom moved by 1.5 A in a direction of its own: its + // restraints let go, the others hold + size_t first = gemmi::count_atom_sites(st.models[0].chains[0]); + for (int k = 0; k < 15; k++) + first += st.models[0].chains[1].residues[k].atoms.size(); + std::vector> moved_xyzb = xyzb; + std::uniform_real_distribution u(-1.0, 1.0); + for (size_t i = first; i < first + st.models[0].chains[1].residues[15].atoms.size(); i++) { + const gemmi::Vec3 v = gemmi::Vec3(u(rng), u(rng), u(rng)).normalized() * 1.5; + moved_xyzb[i][0] += v.x; + moved_xyzb[i][1] += v.y; + moved_xyzb[i][2] += v.z; + } + CHECK(NcsReleasedResidues(ncs, xyzb) == 0); + CHECK(NcsReleasedResidues(ncs, moved_xyzb) >= 1); + CHECK(NcsReleasedResidues(ncs, moved_xyzb) <= 3); +} From 696a066442cb831b69f1b1ffd54664cc73436a3c Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 27 Oct 2026 16:33:31 +0100 Subject: [PATCH 2/2] RESULT-CHANGING: RefinementNcs: copies need 80% identity over the longer chain, not the shorter Over the shorter chain a short RNA aligned with gaps into a long one finds most of its four letters matched: on the 70S ribosome (4ybb) the 5S, 16S and 23S rRNAs fell into one group, and their restraints (1073 residues released) were noise. The test adds a fragment of the molecule that must not be taken for a copy. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/ACKNOWLEDGEMENT.md | 2 +- docs/CPU_DATA_ANALYSIS_REFINEMENT.md | 2 +- .../structure_refinement/RefinementNcs.cpp | 8 ++++++-- .../structure_refinement/RefinementNcs.h | 13 +++++++------ tests/RefinementTest.cpp | 15 ++++++++++++--- 5 files changed, 27 insertions(+), 13 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 10332fba5..5c1e3570d 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -597,7 +597,7 @@ nucleic acids, F. Leonarski, L. D'Ascenzo and P. Auffinger, "Mg2+ ions: do they nitrogens?" (2017), Nucleic Acids Res. 45, 987-1004 [doi:10.1093/nar/gkw1175](https://doi.org/10.1093/nar/gkw1175). **NCS restraints** — copies of a molecule in the model are found by sequence alignment with REFMAC5's criterion -(more than 15 residues aligned, 80% identity; Murshudov et al. 2011, above) and restrained locally, on pairs of +(more than 15 residues aligned, 80% identity - here of the longer chain; Murshudov et al. 2011, above) and restrained locally, on pairs of corresponding interatomic distances under the Geman-McClure robust loss, as REFMAC5's local NCS restraints are, with the distance range of R. A. Nicholls, F. Long and G. N. Murshudov, "Low-resolution refinement tools in REFMAC5" (2012), Acta Cryst. D68, 404-417 [doi:10.1107/S090744491105606X](https://doi.org/10.1107/S090744491105606X), and diff --git a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md index d7168120d..a92164881 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -199,7 +199,7 @@ the ceiling of §15.1. **NCS restraints.** Where the model holds several copies of one molecule, they are found from the model alone: two polymer chains of one kind (protein, or nucleic acid) are copies where their sequences align with more than 15 -residues at 80% identity or more, REFMAC5's criterion. In the residues aligned to a residue of the same name, the +residues identical and identical in 80% of the longer chain or more - REFMAC5's criterion, taken on the longer chain. In the residues aligned to a residue of the same name, the atoms of the same name correspond - not a hydrogen, not an alternative conformation, both refinable. Each copy is paired with the next of its molecule (and the last with the first, where there are three or more), and each pair is restrained **locally**, on its interatomic distances rather than on a superposition: every distance within a chain diff --git a/image_analysis/structure_refinement/RefinementNcs.cpp b/image_analysis/structure_refinement/RefinementNcs.cpp index afdb2e252..7590a7c08 100644 --- a/image_analysis/structure_refinement/RefinementNcs.cpp +++ b/image_analysis/structure_refinement/RefinementNcs.cpp @@ -3,6 +3,7 @@ #include "RefinementNcs.h" +#include #include #include #include @@ -21,7 +22,7 @@ namespace { constexpr double NCS_DISTANCE_MAX = 4.2; // A // Where a pair's restraint lets go (NcsDistanceTerm): copies whose distance differs by this much are taken to differ. constexpr float NCS_RELEASE = 0.5f; // A -// Copies: more than this many residues aligned, at this identity at least. +// Copies: more than this many residues aligned, at this identity at least (of the longer chain). // Following Murshudov et al. (2011) Acta Cryst. D67, 355-367 constexpr int NCS_MIN_ALIGNED = 15; constexpr double NCS_MIN_IDENTITY = 80.0; // % @@ -69,7 +70,10 @@ bool AreCopies(const NcsChain &a, const NcsChain &b) { return false; const gemmi::AlignmentResult al = gemmi::align_string_sequences(a.sequence, b.sequence, {}, gemmi::AlignmentScoring::simple()); - return al.match_count > NCS_MIN_ALIGNED && al.calculate_identity() >= NCS_MIN_IDENTITY; + // identity over the longer chain: over the shorter one, a short RNA aligned with gaps into a long one finds most of + // its four letters matched + const size_t longer = std::max(a.sequence.size(), b.sequence.size()); + return al.match_count > NCS_MIN_ALIGNED && 100.0 * al.match_count >= NCS_MIN_IDENTITY * double(longer); } // The pairs (i, j), i < j, of `atoms` closer than `reach` at `pos`, by a grid of cells of that size. diff --git a/image_analysis/structure_refinement/RefinementNcs.h b/image_analysis/structure_refinement/RefinementNcs.h index 239acb2f9..0d8592ca4 100644 --- a/image_analysis/structure_refinement/RefinementNcs.h +++ b/image_analysis/structure_refinement/RefinementNcs.h @@ -30,12 +30,13 @@ struct NcsCopies { }; // The copies of `model` (whose entities are set up, as BuildRefinementTopology leaves them): polymer chains whose -// sequences align with more than 15 residues and 80% identity (REFMAC's criterion) are copies of one molecule. Their -// residues of the same name aligned to each other correspond, and in them the atoms of the same name - neither a -// hydrogen nor an alternative conformation, all refinable (topo.refinable). Each copy is paired with the next of its -// molecule (the last with the first, where there are three or more), and each pair restrained on the distances -// within a chain shorter than NCS_DISTANCE_MAX in either copy, those the dictionary's bonds and angles fix excepted. -// Atoms are numbered in model order, as ModelPositions() lists them. +// sequences align with more than 15 residues identical, and identical in 80% of the longer chain's residues, are +// copies of one molecule (REFMAC's criterion, taken on the longer chain). Their residues of the same name aligned to +// each other correspond, and in them the atoms of the same name - neither a hydrogen nor an alternative conformation, +// all refinable (topo.refinable). Each copy is paired with the next of its molecule (the last with the first, where +// there are three or more), and each pair restrained on the distances within a chain shorter than NCS_DISTANCE_MAX in +// either copy, those the dictionary's bonds and angles fix excepted. Atoms are numbered in model order, as +// ModelPositions() lists them. NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo); // The restraints at `sigma` (A). diff --git a/tests/RefinementTest.cpp b/tests/RefinementTest.cpp index 1911bef09..ca5229d4f 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -839,7 +839,8 @@ TEST_CASE("Refinement_restrained_recovers_model", "[Refinement]") { } } -// Two copies of one RNA chain in the asymmetric unit, the second turned and moved: found as copies of one molecule, held +// Two copies of one RNA chain in the asymmetric unit, the second turned and moved, and a fragment of it: the two found +// as copies of one molecule, the fragment not, held // to each other by local distance restraints whose gradient is right on the CPU and the GPU, and let go where the // copies really differ. A single copy has none. TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { @@ -862,6 +863,14 @@ TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { for (gemmi::Atom &at : r.atoms) at.pos = gemmi::Position(rot.multiply(at.pos)) + gemmi::Position(40, 10, -5); st.models[0].chains.push_back(copy); + // and a fragment of the molecule, its first 18 residues: not a copy of the whole + gemmi::Chain fragment = st.models[0].chains.at(0); + fragment.name = "C"; + fragment.residues.resize(18); + for (gemmi::Residue &r : fragment.residues) + for (gemmi::Atom &at : r.atoms) + at.pos += gemmi::Position(-30, 40, 20); + st.models[0].chains.push_back(fragment); st.connections.clear(); // the LINK records name chain A only const RefinementTopology topo = BuildRefinementTopology(st); const NcsCopies ncs = FindNcsCopies(st.models[0], topo); @@ -888,8 +897,8 @@ TEST_CASE("RefinementNcs_two_copies", "[Refinement]") { // the second copy shaken: the gradient against central differences std::mt19937 rng(7); std::normal_distribution d(0, 0.1); - const size_t half = xyzb.size() / 2; - for (size_t i = half; i < xyzb.size(); i++) + const size_t half = gemmi::count_atom_sites(st.models[0].chains[0]); + for (size_t i = half; i < 2 * half; i++) for (int c = 0; c < 3; c++) xyzb[i][c] += d(rng); std::vector> grad(xyzb.size(), {0, 0, 0, 0});