Robust sort_backbone: detect disconnected selections and avoid infinite loops (#5113)

* fix #5112
* analysis.polymer: sort_backbone can hang indefinitely when the input AtomGroup is disconnected 
  (e.g., some backbone atoms were excluded by selection)
* changed logic for sort_backbone
* added test for internal missing
* update CHANGELOG

Author: DrDomenicoMarson

package/CHANGELOG

@@ -43,6 +43,8 @@ Fixes
  * Fixes the benchmark `SimpleRmsBench` in `benchmarks/analysis/rms.py`
    by changing the way weights for RMSD are calculated, instead of
    directly passing them. (Issue #3520, PR #5006)
+ * `analysis.polymer.sort_backbone` is now working for discontinuous polymers
+    and does not result in an infinite loop (Issue #5112, PR #5113)
 
 Enhancements
  * Added support for reading and processing streamed data in `coordinates.base`
@@ -90,6 +92,7 @@ Changes
    `analysis.lineardensity.LinearDensity.totalmass` (PR #5007)
  * Remove `default` channel from RTD conda env. (Issue # 5036, PR # 5037)   
 
+
 Deprecations
  * The RDKit converter parameter `NoImplicit` has been deprecated in favour of
    `implicit_hydrogens` and `inferrer` parameters. `max_iter` has been moved

package/MDAnalysis/analysis/polymer.py

@@ -67,37 +67,32 @@ def sort_backbone(backbone):
 
     .. versionadded:: 0.20.0
     """
-    if not backbone.n_fragments == 1:
-        raise ValueError(
-            "{} fragments found in backbone.  "
-            "backbone must be a single contiguous AtomGroup"
-            "".format(backbone.n_fragments)
-        )
+    degrees = [len(atom.bonded_atoms & backbone) for atom in backbone]
+    deg1_atoms = [atom for atom, d in zip(backbone, degrees) if d == 1]
+    wrong_atoms = [
+        atom for atom, d in zip(backbone, degrees) if d not in (1, 2)
+    ]
 
-    branches = [at for at in backbone if len(at.bonded_atoms & backbone) > 2]
-    if branches:
-        # find which atom has too many bonds for easier debug
+    if len(wrong_atoms) > 0:
         raise ValueError(
-            "Backbone is not linear.  "
-            "The following atoms have more than two bonds in backbone: {}."
-            "".format(",".join(str(a) for a in branches))
+            "Backbone contains atoms with connectivity degree not equal to 1 or 2. "
+            "This suggests branches or isolated atoms. Problematic atoms: {}."
+            "".format(", ".join(str(a) for a in wrong_atoms))
         )
 
-    caps = [
-        atom for atom in backbone if len(atom.bonded_atoms & backbone) == 1
-    ]
-    if not caps:
-        # cyclical structure
+    if len(deg1_atoms) != 2:
         raise ValueError(
-            "Could not find starting point of backbone, "
-            "is the backbone cyclical?"
+            "Backbone connectivity invalid: "
+            "expected exactly 2 atoms with connectivity degree 1 (caps). "
+            "Cyclical structures are not supported. "
+            "Atoms with connectivity degree 1 found: {}."
+            "".format(", ".join(str(a) for a in deg1_atoms))
         )
 
     # arbitrarily choose one of the capping atoms to be the startpoint
-    sorted_backbone = AtomGroup([caps[0]])
+    sorted_backbone = AtomGroup([deg1_atoms[0]])
 
-    # iterate until the sorted chain length matches the backbone size
-    while len(sorted_backbone) < len(backbone):
+    for _ in range(len(backbone) - 1):
         # current end of the chain
         end_atom = sorted_backbone[-1]
 

testsuite/MDAnalysisTests/analysis/test_persistencelength.py

@@ -149,25 +149,31 @@ def test_sortbb(self, u):
 
         assert_equal(s_ag.ids, [0, 1, 4, 6, 8])
 
-    def test_not_fragment(self, u):
-        # two fragments don't work
-        bad_ag = u.residues[0].atoms[:2] + u.residues[1].atoms[:2]
-        with pytest.raises(ValueError):
-            polymer.sort_backbone(bad_ag)
-
     def test_branches(self, u):
         # includes side branches, can't sort
         bad_ag = u.atoms[:10]  # include -H etc
 
-        with pytest.raises(ValueError):
+        with pytest.raises(ValueError, match="branches or isolated atoms"):
             polymer.sort_backbone(bad_ag)
 
+    def test_isolated(self, u):
+        u = mda.Universe.empty(4, trajectory=True)
+        bondlist = [(0, 1), (1, 2)]
+        u.add_TopologyAttr(Bonds(bondlist))
+        with pytest.raises(ValueError, match="branches or isolated atoms"):
+            polymer.sort_backbone(u.atoms)
+
+    def test_missing_internal(self, u):
+        u = mda.Universe.empty(4, trajectory=True)
+        bondlist = [(0, 1), (2, 3)]
+        u.add_TopologyAttr(Bonds(bondlist))
+        with pytest.raises(ValueError, match="Backbone connectivity invalid"):
+            polymer.sort_backbone(u.atoms)
+
     def test_circular(self):
         u = mda.Universe.empty(6, trajectory=True)
         # circular structure
         bondlist = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 0)]
         u.add_TopologyAttr(Bonds(bondlist))
-
-        with pytest.raises(ValueError) as ex:
+        with pytest.raises(ValueError, match="Cyclical"):
             polymer.sort_backbone(u.atoms)
-        assert "cyclical" in str(ex.value)