Merge ncs-fix into rc181: NCS copies also within 2.5 A C-alpha rms after superposition

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

docs/ACKNOWLEDGEMENT.md

@@ -631,7 +631,8 @@ what the data support rather than a rule by resolution: E. A. Merritt, "To B or
 resolution?" (2012), Acta Cryst. D68, 468-477 [doi:10.1107/S0907444911028320](https://doi.org/10.1107/S0907444911028320).
 
 **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
+(more than 15 residues aligned, 80% identity - here of the longer chain - and 2.5 Å rms after superposition;
+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

docs/CPU_DATA_ANALYSIS_REFINEMENT.md

@@ -236,11 +236,13 @@ 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
+residues identical and identical in 80% of the longer chain or more, and their C-alpha (P) atoms lie within 2.5 Å
+rms of each other after superposition - REFMAC5's criterion, the identity taken on the longer chain. The same
+sequence in a conformation further away than that, such as a quasi-equivalent subunit of a virus capsid, is another
+molecule. 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,
 
 $$

image_analysis/structure_refinement/RefinementNcs.cpp

@@ -5,6 +5,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <limits>
 #include <map>
 #include <set>
 #include <unordered_map>
@@ -26,6 +27,10 @@ constexpr float NCS_RELEASE = 0.5f;        // A
 // Following Murshudov et al. (2011) Acta Cryst. D67, 355-367
 constexpr int NCS_MIN_ALIGNED = 15;
 constexpr double NCS_MIN_IDENTITY = 80.0;  // %
+// and their C-alpha (P) atoms no further apart than this after superposition: the same sequence further away (a
+// quasi-equivalent subunit of a virus capsid) is another molecule.
+// Following Murshudov et al. (2011) Acta Cryst. D67, 355-367
+constexpr double NCS_MAX_RMSD = 2.5;       // A
 
 // One polymer chain: its residues and, per residue, its atoms' indices by name.
 struct NcsChain {
@@ -65,15 +70,30 @@ void Correspond(const NcsChain &a, const NcsChain &b, std::vector<std::array<int
         }
 }
 
-bool AreCopies(const NcsChain &a, const NcsChain &b) {
+// The C-alpha (P) rmsd of two chains after superposition where their sequences make them copies of one molecule,
+// infinity where they do not.
+double CopyRmsd(const NcsChain &a, const NcsChain &b, const std::vector<gemmi::Position> &pos) {
+    const double none = std::numeric_limits<double>::infinity();
     if (a.type != b.type || a.sequence.size() <= NCS_MIN_ALIGNED || b.sequence.size() <= NCS_MIN_ALIGNED)
-        return false;
+        return none;
     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);
+    if (al.match_count <= NCS_MIN_ALIGNED || 100.0 * al.match_count < NCS_MIN_IDENTITY * double(longer))
+        return none;
+    std::vector<std::array<int, 2>> atoms, anchors;
+    std::vector<int> residue_of;
+    Correspond(a, b, atoms, residue_of, anchors);
+    if (anchors.size() < 3)
+        return none;
+    std::vector<gemmi::Position> pa, pb;
+    for (const std::array<int, 2> &p : anchors) {
+        pa.push_back(pos[p[0]]);
+        pb.push_back(pos[p[1]]);
+    }
+    return gemmi::superpose_positions(pa.data(), pb.data(), pa.size(), nullptr).rmsd;
 }
 
 // The pairs (i, j), i < j, of `atoms` closer than `reach` at `pos`, by a grid of cells of that size.
@@ -146,7 +166,7 @@ NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &top
     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])) {
+            if (CopyRmsd(chains[g[0]], chains[k], pos) <= NCS_MAX_RMSD) {
                 g.push_back(static_cast<int>(k));
                 placed = true;
                 break;
@@ -154,7 +174,6 @@ NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &top
         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)

image_analysis/structure_refinement/RefinementNcs.h

@@ -30,13 +30,15 @@ struct NcsCopies {
 };
 
 // The copies of `model` (whose entities are set up, as BuildRefinementTopology leaves them): polymer chains whose
-// sequences align with more than 15 residues 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.
+// sequences align with more than 15 residues identical, identical in 80% of the longer chain's residues, and whose
+// C-alpha (P) atoms are within 2.5 A rms of each other after superposition are copies of one molecule (REFMAC's
+// criterion, the identity taken on the longer chain); the same sequence in a conformation further away - a
+// quasi-equivalent subunit of a virus capsid - is another molecule. Their residues of the same name aligned to each
+// other correspond, and in them the atoms of the same name - neither a hydrogen nor an alternative conformation, all
+// refinable (topo.refinable). Each copy is paired with the next of its molecule (the last with the first, where there
+// are three or more), and each pair restrained on the distances within a chain shorter than NCS_DISTANCE_MAX in either
+// copy, those the dictionary's bonds and angles fix excepted. Atoms are numbered in model order, as ModelPositions()
+// lists them.
 NcsCopies FindNcsCopies(const gemmi::Model &model, const RefinementTopology &topo);
 
 // The restraints at `sigma` (A).

tests/RefinementTest.cpp

@@ -997,8 +997,8 @@ TEST_CASE("Refinement_ion_form_factors", "[Refinement]") {
     }
 }
 
-// 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
+// Two copies of one RNA chain in the asymmetric unit, the second turned and moved, a fragment of it and the chain in
+// another conformation: the two found as copies of one molecule, the fragment and the other conformation 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]") {
@@ -1029,6 +1029,17 @@ TEST_CASE("RefinementNcs_two_copies", "[Refinement]") {
         for (gemmi::Atom &at : r.atoms)
             at.pos += gemmi::Position(-30, 40, 20);
     st.models[0].chains.push_back(fragment);
+    // and the whole molecule in another conformation, each residue moved 3 A its own way: the same sequence, not a copy
+    gemmi::Chain other = st.models[0].chains.at(0);
+    other.name = "D";
+    std::mt19937 shift_rng(11);
+    std::uniform_real_distribution<double> shift(-1.0, 1.0);
+    for (gemmi::Residue &r : other.residues) {
+        const gemmi::Vec3 v = gemmi::Vec3(shift(shift_rng), shift(shift_rng), shift(shift_rng)).normalized() * 3.0;
+        for (gemmi::Atom &at : r.atoms)
+            at.pos += gemmi::Position(v) + gemmi::Position(30, -40, 25);
+    }
+    st.models[0].chains.push_back(other);
     st.connections.clear();   // the LINK records name chain A only
     const RefinementTopology topo = BuildRefinementTopology(st);
     const NcsCopies ncs = FindNcsCopies(st.models[0], topo);