From 68c2509df869f6b60def5039467c326b09c918be Mon Sep 17 00:00:00 2001 From: Manuel Nuno Melo Date: Thu, 10 Jun 2021 12:58:53 +0100 Subject: [PATCH] Fixed wrong atom sorting in _split_by_compound_indices Closes #3353 --- package/CHANGELOG | 2 ++ package/MDAnalysis/core/groups.py | 12 +++++---- testsuite/MDAnalysisTests/core/test_unwrap.py | 26 ++++++++++++++++++- 3 files changed, 34 insertions(+), 6 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 6b2eb40f379..3ddbb121bcf 100644 --- a/package/CHANGELOG +++ b/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 diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 796c2a7c151..55324188bea 100644 --- a/package/MDAnalysis/core/groups.py +++ b/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 = [] diff --git a/testsuite/MDAnalysisTests/core/test_unwrap.py b/testsuite/MDAnalysisTests/core/test_unwrap.py index b72f435395d..8749ecae5bc 100644 --- a/testsuite/MDAnalysisTests/core/test_unwrap.py +++ b/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)