Merge ncs-refine into refine-merge: NCS local distance restraints with the solvent build in the trajectory

The solvent build now runs inside RefineOnce's trajectory (solvent-p1), so the NCS trials were carried into the
trajectory lambda: every step adds its "NCS changed" trial, and the solvent-only judge of a build round runs at the
NCS sigma the trajectory ended at. The copies are found once per refinement on the polymer atoms; the solvent
placeholders (occupancy 0, own chain) and relabelled waters/ions are never part of them, and atom indices do not
move when the build edits the slots in place, so the restraint records stay valid through the rounds.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

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

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).

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 |

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

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, <cos(dphi)>, 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<int> OwnSlotDevices(const Context &main, int wanted, const Refinemen
     return std::vector<int>(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<int> &refinable, const std::vector<gemmi::Position> &pos0,
                       const std::vector<double> &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<int> &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<int> &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<int> &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<int> &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<TrialOutcome> RunTrials(const std::vector<TrialSlot> &trial_slots,
-                                    const std::vector<std::pair<double, double>> &trials,
+                                    const std::vector<TrialWeights> &trials,
                                     const std::vector<int> &refinable, const std::vector<gemmi::Position> &pos0,
                                     const std::vector<double> &b0, const Context::Fit *fit, int cycles,
-                                    const RefinementSettings &settings) {
+                                    const RefinementSettings &settings, const NcsCopies *ncs) {
     std::vector<TrialOutcome> out(trials.size());
     std::vector<std::exception_ptr> 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<TrialOutcome> RunTrials(const std::vector<TrialSlot> &trial_slots,
     std::vector<size_t> 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<TrialOutcome> RunTrials(const std::vector<TrialSlot> &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<Me
     }
     pool.first = static_cast<int>(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<Me
     // on the main model and engine, and on copies with engines of their own; each with a share of the threads. The
     // copies are made while the placement runs: they start from the model as it was before the topology was made
     // and take nothing from the main context that the placement changes.
-    int wanted = static_cast<int>(std::max({settings.weight_ladder.size() * settings.adp_ladder.size(),
+    int wanted = static_cast<int>(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<int>(settings.max_trials_at_once));
@@ -1584,6 +1609,17 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
                     topo.fixed_names.empty() ? std::string() : " (" + Joined(topo.fixed_names) + ")",
                     topo.bonds, topo.metal_bonds.records.size(), topo.metal_links, topo.metal_links_refused, topo.angles, topo.torsions, topo.chirals, topo.planes);
 
+        if (ncs_p != nullptr) {
+            result.ncs_groups = ncs.groups;
+            result.ncs_copies = ncs.copies;
+            result.ncs_restraints = static_cast<int>(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<Me
         const double w_max = w0 * settings.max_weight;
 
         // 1. The weights from the placed model: the ladder's rungs, tighter geometry first, each with the B factors
-        //    held and refined at their own balanced weight.
-        std::vector<std::pair<double, double>> ladder;
+        //    held and refined at their own balanced weight; where there are copies, each without NCS restraints and
+        //    then with them.
+        std::vector<TrialWeights> 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<TrialOutcome> 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<Me
         // 2. The trajectory: steps of a few macro-cycles from where the last one ended, each trying the weights as they
         //    are and, in turn, the positions' X-ray weight halved and doubled (halved: the geometry loosens as the
         //    model improves at one weight, until it leaves the bounds) and the B factors' divided and multiplied by
-        //    four.
+        //    four; where there are copies, also the NCS restraints changed: switched on where they are off, and in
+        //    turn tightened and loosened (to off beyond NCS_SIGMA_MAX) where they are on.
         //    A step is taken while it lowers R-free or R-work by more than the tie and its R-free stays within the tie
         //    of the best so far; the trajectory ends at its best step (on a deposited model, whose free set it was
         //    fitted to, R-work falls while R-free rises, and it ends where it began).
@@ -1667,19 +1708,23 @@ RefinementResult RefineOnce(const gemmi::Structure &placed, const std::vector<Me
             from_best.clear();
             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<std::pair<double, double>> trials = {{w, adp}};
+                const double w = state.trial.weight, adp = state.trial.adp_weight, ncs_s = state.trial.ncs_sigma;
+                std::vector<TrialWeights> 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<TrialOutcome> 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<Me
                 solvent_settings.xyz_iterations = SOLVENT_ITERATIONS;
                 solvent_settings.b_iterations = SOLVENT_ITERATIONS;
                 std::vector<TrialOutcome> 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<Me
         // 3. The final weights: a short ladder from where the trajectory ended, tighter geometry first, its choice
         //    kept where it displaces that. A rung the step from there tried already is taken from that step.
         if (state.ran && !ctx.Cancelled()) {
-            const double w = state.trial.weight, adp = state.trial.adp_weight;
-            std::vector<std::pair<double, double>> trials;
+            const double w = state.trial.weight, adp = state.trial.adp_weight, ncs_s = state.trial.ncs_sigma;
+            std::vector<TrialWeights> 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<std::pair<double, double>> missing;
+            std::vector<TrialWeights> 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<TrialOutcome> 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<Me
             result.restrained = true;
             result.weight = state.trial.weight;
             result.adp_weight = state.trial.adp_weight;
+            result.ncs_sigma = state.trial.ncs_sigma;
+            if (ncs_p != nullptr) {
+                result.ncs_rmsd = NcsRmsd(ncs, ctx.XyzB());
+                result.ncs_released = NcsReleasedResidues(ncs, ctx.XyzB());
+                logger.Info("Refinement: NCS {}; copies {:.3f} A apart (as placed {:.3f} A), {} residues released",
+                            result.ncs_sigma > 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;

image_analysis/structure_refinement/Refinement.h

@@ -60,6 +60,10 @@ struct RefinementSettings {
     std::vector<double> 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<RefinementWeightTrial> 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;

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],

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<T> a, V3<T> b, T d0, T sigma, V3<T> 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 <typename T>
+REFINE_HD inline T NcsDistanceTerm(V3<T> a, V3<T> b, V3<T> c, V3<T> d, T sigma, T release, V3<T> g[4]) {
+    const V3<T> 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.

image_analysis/structure_refinement/RefinementNcs.cpp

@@ -0,0 +1,253 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "RefinementNcs.h"
+
+#include <algorithm>
+#include <cmath>
+#include <map>
+#include <set>
+#include <unordered_map>
+
+#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<std::string> sequence;
+    std::vector<int> residue;                           // numbered over the model
+    std::vector<std::map<std::string, int>> atoms;      // name -> atom index
+    std::vector<int> 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<std::array<int, 2>> &atoms,
+                std::vector<int> &residue_of, std::vector<std::array<int, 2>> &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<std::pair<int, int>> PairsWithin(const std::vector<int> &atoms, const std::vector<gemmi::Position> &pos,
+                                          double reach) {
+    std::map<std::array<int, 3>, std::vector<int>> cells;   // cell -> positions in `atoms`
+    auto cell_of = [&](const gemmi::Position &p) {
+        return std::array<int, 3>{static_cast<int>(std::floor(p.x / reach)), static_cast<int>(std::floor(p.y / reach)),
+                                  static_cast<int>(std::floor(p.z / reach))};
+    };
+    for (size_t k = 0; k < atoms.size(); k++)
+        cells[cell_of(pos[atoms[k]])].push_back(static_cast<int>(k));
+    std::set<std::pair<int, int>> out;
+    for (size_t k = 0; k < atoms.size(); k++) {
+        const std::array<int, 3> 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<gemmi::Position> pos;
+    std::vector<NcsChain> 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<std::string, int> 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<std::vector<int>> groups;
+    for (size_t k = 0; k < chains.size(); k++) {
+        bool placed = false;
+        for (std::vector<int> &g : groups)
+            if (AreCopies(chains[g[0]], chains[k])) {
+                g.push_back(static_cast<int>(k));
+                placed = true;
+                break;
+            }
+        if (!placed)
+            groups.push_back({static_cast<int>(k)});
+    }
+
+    // The distances the dictionary fixes: bonded atoms and those bonded to a common atom.
+    std::set<std::pair<int, int>> fixed;
+    for (const RestraintRecord &r : topo.geometry.records)
+        if (r.kind == static_cast<int>(RestraintKind::Bond))
+            fixed.insert(std::minmax(r.atom[0], r.atom[1]));
+        else if (r.kind == static_cast<int>(RestraintKind::Angle))
+            fixed.insert(std::minmax(r.atom[0], r.atom[2]));
+
+    for (const std::vector<int> &g : groups) {
+        if (g.size() < 2)
+            continue;
+        out.groups++;
+        out.copies += static_cast<int>(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<std::array<int, 2>> corr, anchors;
+            std::vector<int> residue_of;
+            Correspond(a, b, corr, residue_of, anchors);
+            out.anchors.push_back(std::move(anchors));
+            std::unordered_map<int, int> to_b, res_of;
+            std::vector<int> 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<std::pair<int, int>> pairs = PairsWithin(atoms_a, pos, NCS_DISTANCE_MAX);
+            std::unordered_map<int, int> 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<int>(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<float>(sigma);
+    return set;
+}
+
+double NcsRmsd(const NcsCopies &ncs, const std::vector<std::array<double, 4>> &xyzb) {
+    double sum = 0;
+    int n = 0;
+    for (const std::vector<std::array<int, 2>> &anchors : ncs.anchors) {
+        if (anchors.size() < 3)
+            continue;
+        std::vector<gemmi::Position> pa, pb;
+        for (const std::array<int, 2> &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<std::array<double, 4>> &xyzb) {
+    std::map<int, std::array<int, 2>> 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<int, 2> &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;
+}

image_analysis/structure_refinement/RefinementNcs.h

@@ -0,0 +1,49 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <array>
+#include <string>
+#include <vector>
+
+#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<int> 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<std::vector<std::array<int, 2>>> 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<std::array<double, 4>> &xyzb);
+
+// The residues let go: those whose restraints are more than half beyond their release distance.
+int NcsReleasedResidues(const NcsCopies &ncs, const std::vector<std::array<double, 4>> &xyzb);

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<std::array<
                                             At(xyzb, r.atom[3]), ideal, sigma, g);
                 n = 4;
                 break;
+            case RestraintKind::NcsDistance:
+                e = refine_math::NcsDistanceTerm(At(xyzb, r.atom[0]), At(xyzb, r.atom[1]), At(xyzb, r.atom[2]),
+                                                 At(xyzb, r.atom[3]), sigma, ideal, g);
+                n = 4;
+                break;
             case RestraintKind::Repulsion:
                 e = refine_math::RepulsionTerm(At(xyzb, r.atom[0]), Image(r, At(xyzb, r.atom[1])), ideal, sigma, g);
                 g[1] = FromImage(r, g[1]);

rugnux/ResultReport.cpp

@@ -554,6 +554,16 @@ namespace {
                 Add(s, KeyReal("REFINE_RMSD_ANGLES_START", rf.geometry_start.rmsd_angles, "{:.2f}", true));
                 Add(s, KeyInt("REFINE_CLASHES", rf.geometry.clashes));
                 Add(s, KeyReal("REFINE_XRAY_WEIGHT", rf.weight, "{:.4g}", true));
+                if (rf.ncs_groups > 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));

tests/RefinementTest.cpp

@@ -13,6 +13,7 @@
 
 #include <gemmi/mmread_gz.hpp>
 #include <gemmi/neighbor.hpp>
+#include <gemmi/modify.hpp>
 #include <gemmi/symmetry.hpp>
 
 #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<gemmi::Atom *> atoms = atoms_of(st.models[0]);
+    std::vector<std::array<double, 4>> 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<double> 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<std::array<double, 4>> 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<MergedReflection> 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<std::array<double, 4>> 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<std::array<double, 4>> moved_xyzb = xyzb;
+    std::uniform_real_distribution<double> 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);
+}