diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 7d32e9afb..8ecc19a82 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 - 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 +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 bf5c8dbb1..678375da7 100644 --- a/docs/CPU_DATA_ANALYSIS_REFINEMENT.md +++ b/docs/CPU_DATA_ANALYSIS_REFINEMENT.md @@ -211,6 +211,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 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 +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`), @@ -262,6 +282,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. @@ -443,7 +469,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 beyond the six species**: waters and ions are relabelled among water, Mg²⁺, Na⁺, K⁺, Ca²⁺ and Cl⁻ (§15.10); NH₄⁺ is only reported, transition metals (Zn, Mn, ...) are measured but never relabelled, and no @@ -479,6 +507,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 12fa07cf6..45a15b664 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 ce3d4e757..9cf473782 100644 --- a/image_analysis/structure_refinement/Refinement.cpp +++ b/image_analysis/structure_refinement/Refinement.cpp @@ -36,6 +36,7 @@ #include "MonomerLibrary.h" #include "RefinementMath.h" #include "RefinementMinimiser.h" +#include "RefinementNcs.h" #include "RefinementXray.h" #ifdef JFJOCH_USE_CUDA #include "../../common/CUDAWrapper.h" @@ -62,6 +63,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. @@ -1023,12 +1026,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; } @@ -1130,6 +1136,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 { @@ -1143,7 +1155,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"); @@ -1167,7 +1179,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); } @@ -1193,6 +1205,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; @@ -1203,9 +1216,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(), @@ -1221,10 +1239,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) @@ -1236,7 +1254,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 @@ -1244,8 +1262,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(); } @@ -1439,6 +1457,12 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector(ctx.atoms.size()) - pool.n; // the placeholders are the last atoms (no hydrogens) + // The copies of a molecule, restrained to each other in trials beside those that refine them independently. + NcsCopies ncs; + if (settings.restrained && settings.ncs) + ncs = FindNcsCopies(ctx.Model(), topo); + const NcsCopies *ncs_p = ncs.Empty() ? nullptr : &ncs; + MakeEngine(ctx, settings.gpu, settings, logger); result.on_gpu = ctx.engine->OnGPU(); @@ -1528,7 +1552,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)); @@ -1584,6 +1609,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); @@ -1639,13 +1675,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) @@ -1654,7 +1694,8 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const 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}); } std::vector outcomes = RunTrials(trial_slots, trials, refinable, state.pos, state.b, fresh && step == 0 ? nullptr : &state.fit, - settings.macro_cycles, settings); + settings.macro_cycles, settings, ncs_p); if (state_is_best) from_best = outcomes; const int k = choose(outcomes, fmt::format("step {}", step + 1).c_str()); @@ -1869,8 +1914,9 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector local = - RunTrials(trial_slots, {{state.trial.weight, state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0}}, - solvent, changed_pos, changed_b, nullptr, settings.macro_cycles, solvent_settings); + RunTrials(trial_slots, {{state.trial.weight, state.trial.adp_weight > 0 ? state.trial.adp_weight : adp_w0, + state.trial.ncs_sigma}}, + solvent, changed_pos, changed_b, nullptr, settings.macro_cycles, solvent_settings, ncs_p); logger.Info("Solvent build, round {}: the solvent refined alone: {}", round, local[0].line); TrialOutcome after = std::move(local[0]); if (!after.ran || !(after.trial.r_free < state.trial.r_free)) { @@ -1924,22 +1970,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++]); @@ -1962,6 +2009,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 d5c245401..4fa8f5821 100644 --- a/image_analysis/structure_refinement/Refinement.h +++ b/image_analysis/structure_refinement/Refinement.h @@ -60,6 +60,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, and waters and ions // relabelled, within the weight trajectory (SolventBuild.h), to SOLVENT_BUILD_D_MIN; off (rugnux // --no-refine-solvent) the solvent stays as given and is only audited. @@ -112,6 +116,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 @@ -140,6 +146,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..7590a7c08 --- /dev/null +++ b/image_analysis/structure_refinement/RefinementNcs.cpp @@ -0,0 +1,253 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "RefinementNcs.h" + +#include +#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 (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; // % + +// 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()); + // 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. +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..0d8592ca4 --- /dev/null +++ b/image_analysis/structure_refinement/RefinementNcs.h @@ -0,0 +1,49 @@ +// 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 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). +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 5397c4ca4..dd322e632 100644 --- a/image_analysis/structure_refinement/RefinementRestraints.cpp +++ b/image_analysis/structure_refinement/RefinementRestraints.cpp @@ -352,6 +352,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; } @@ -387,6 +390,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 a2573d925..fcad9f4f4 100644 --- a/tests/RefinementTest.cpp +++ b/tests/RefinementTest.cpp @@ -13,6 +13,7 @@ #include #include +#include #include #include "../common/CUDAWrapper.h" @@ -22,6 +23,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 @@ -967,3 +969,116 @@ TEST_CASE("Refinement_ion_form_factors", "[Refinement]") { CHECK(gpu.r_work_start == Catch::Approx(cpu.r_work_start).margin(1e-4)); } } + +// 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]") { + 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); + // 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); + 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 = 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}); + 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); +}