Fixed wrong atom sorting in _split_by_compound_indices

package/CHANGELOG

@@ -23,6 +23,8 @@ The rules for this file:
   * 2.0.0
 
 Fixes
+  * Fixed sometimes wrong sorting of atoms into fragments when unwrapping
+    (Issue #3352, PR #3353)
   * Fixed missing array flatten preventing PCA from running when mean positions  
     provided and bugs causing tests to incorrectly pass (Issue #2728, PR #3296)
   * DL_POLY HISTORY file may contain unit cell dimensions even if there are

package/MDAnalysis/core/groups.py

@@ -847,6 +847,11 @@ def _split_by_compound_indices(self, compound, stable_sort=False):
         # essentially a non-split?
 
         compound_indices = self._get_compound_indices(compound)
+        compound_sizes = np.bincount(compound_indices)
+        size_per_atom = compound_sizes[compound_indices]
+        compound_sizes = compound_sizes[compound_sizes != 0]
+        unique_compound_sizes = unique_int_1d(compound_sizes)
+
         # Are we already sorted? argsorting and fancy-indexing can be expensive
         # so we do a quick pre-check.
         needs_sorting = np.any(np.diff(compound_indices) < 0)
@@ -858,11 +863,8 @@ def _split_by_compound_indices(self, compound, stable_sort=False):
             else:
                 # Quicksort
                 sort_indices = np.argsort(compound_indices)
-
-        compound_sizes = np.bincount(compound_indices)
-        size_per_atom = compound_sizes[compound_indices]
-        compound_sizes = compound_sizes[compound_sizes != 0]
-        unique_compound_sizes = unique_int_1d(compound_sizes)
+            # We must sort size_per_atom accordingly (Issue #3352).
+            size_per_atom = size_per_atom[sort_indices]
 
         compound_masks = []
         atom_masks = []

testsuite/MDAnalysisTests/core/test_unwrap.py

@@ -21,12 +21,14 @@
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 #
 import numpy as np
-from numpy.testing import assert_almost_equal, assert_array_equal
+from numpy.testing import (assert_raises, assert_almost_equal,
+                           assert_array_equal)
 import pytest
 
 import MDAnalysis as mda
 from MDAnalysis import NoDataError
 from MDAnalysisTests.core.util import UnWrapUniverse
+from MDAnalysis.tests.datafiles import CONECT
 
 
 class TestUnwrap(object):
@@ -356,3 +358,25 @@ def test_unwrap_no_molnums_exception_safety(self, level, reference):
             group.unwrap(compound='molecules', reference=reference,
                          inplace=True)
         assert_array_equal(group.atoms.positions, orig_pos)
+
+
+def test_uncontiguous():
+    """Real-life case of fragment sparsity that triggers Issue 3352
+    """
+    precision = 5
+    displacement_vec = [14.7, 0., 0.]
+    u = mda.Universe(CONECT)
+    # This is one of the few residues that has bonds
+    ag = u.residues[66].atoms
+    ref_pos = ag.positions
+    # Let's break it by placing it over the box boundary and re-packing
+    u.atoms.positions -= displacement_vec
+    u.atoms.pack_into_box()
+    # Let's make sure we really broke the fragment
+    assert_raises(AssertionError, assert_almost_equal,
+                  ref_pos, ag.positions+displacement_vec,
+                  decimal=precision)
+    # Ok, let's make it whole again and check that we're good
+    u.atoms.unwrap()
+    assert_almost_equal(ref_pos, ag.positions+displacement_vec,
+                        decimal=precision)