diff --git a/package/CHANGELOG b/package/CHANGELOG index 75dcd5f534d..d0049fc206e 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -40,6 +40,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 capability to calculate MSD from frames with irregular (non-linear) @@ -82,6 +84,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 diff --git a/package/MDAnalysis/analysis/polymer.py b/package/MDAnalysis/analysis/polymer.py index 2bee3839710..7a0510eab03 100644 --- a/package/MDAnalysis/analysis/polymer.py +++ b/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] diff --git a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py b/testsuite/MDAnalysisTests/analysis/test_persistencelength.py index e17482dd3d4..54005c8a667 100644 --- a/testsuite/MDAnalysisTests/analysis/test_persistencelength.py +++ b/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)