diff --git a/package/CHANGELOG b/package/CHANGELOG index 4589c6a5fec..812baead045 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -15,7 +15,7 @@ The rules for this file: ------------------------------------------------------------------------------ ??/??/20 richardjgowers, IAlibay, orbeckst, tylerjereddy, jbarnoud, yuxuanzhuang, lilyminium, VOD555, p-j-smith, bieniekmateusz, - calcraven, ianmkenney, rcrehuet + calcraven, ianmkenney, rcrehuet, manuel.nuno.melo * 1.0.1 @@ -25,6 +25,8 @@ Fixes * Testsuite does not use any more matplotlib.use('agg') (#2191) * The methods provided by topology attributes now appear in the documentation (Issue #1845) + * AtomGroup.center now works correctly for compounds + unwrapping + (Issue #2984) * In ChainReader, read_frame does not trigger change of iterating position. (Issue #2723, PR #2815) * empty_atomgroup.select_atoms('name *') now returns an empty diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 74edf05ab74..d47cdb35dd6 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -851,12 +851,18 @@ def center(self, weights, pbc=False, compound='group', unwrap=False): # Sort positions and weights by compound index and promote to dtype if # required: - sort_indices = np.argsort(compound_indices) + + # are we already sorted? argsorting and fancy-indexing can be expensive + if np.any(np.diff(compound_indices) < 0): + sort_indices = np.argsort(compound_indices) + else: + sort_indices = slice(None) compound_indices = compound_indices[sort_indices] # Unwrap Atoms if unwrap: - coords = atoms.unwrap(compound=comp, reference=None, inplace=False) + coords = atoms.unwrap(compound=comp, reference=None, + inplace=False)[sort_indices] else: coords = atoms.positions[sort_indices] if weights is None: diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index 24f8a0bd524..1069d4b681a 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -993,6 +993,12 @@ def ag(self): group.wrap(inplace=True) return group + @pytest.fixture() + def unordered_ag(self, ag): + ndx = np.arange(len(ag)) + np.random.shuffle(ndx) + return ag[ndx] + @pytest.fixture() def ref_noUnwrap_residues(self): return { @@ -1015,12 +1021,12 @@ def ref_Unwrap_residues(self): 'COG': np.array([[21.356, 41.685, 40.501], [44.577, 43.312, 79.039], [ 2.204, 27.722, 54.023]], dtype=np.float32), - 'COM': np.array([[20.815, 42.013, 39.802], - [44.918, 43.282, 79.325], - [2.045, 28.243, 54.127]], dtype=np.float32), - 'MOI': np.array([[16747.486, -1330.489, 2938.243], - [-1330.489, 19315.253, 3306.212], - [ 2938.243, 3306.212, 8990.481]]), + 'COM': np.array([[21.286, 41.664, 40.465], + [44.528, 43.426, 78.671], + [ 2.111, 27.871, 53.767]], dtype=np.float32), + 'MOI': np.array([[16687.941, -1330.617, 2925.883], + [-1330.617, 19256.178, 3354.832], + [ 2925.883, 3354.832, 8989.946]]), 'Asph': 0.2969491080, } @@ -1060,6 +1066,12 @@ def test_UnWrapFlag_residues(self, ag, ref_Unwrap_residues): assert_almost_equal(ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec) assert_almost_equal(ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec) + def test_UnWrapFlag_residues_unordered(self, unordered_ag, ref_Unwrap_residues): + assert_almost_equal(unordered_ag.center_of_geometry(unwrap=True, compound='residues'), ref_Unwrap_residues['COG'], self.prec) + assert_almost_equal(unordered_ag.center_of_mass(unwrap=True, compound='residues'), ref_Unwrap_residues['COM'], self.prec) + assert_almost_equal(unordered_ag.moment_of_inertia(unwrap=True, compound='residues'), ref_Unwrap_residues['MOI'], self.prec) + assert_almost_equal(unordered_ag.asphericity(unwrap=True, compound='residues'), ref_Unwrap_residues['Asph'], self.prec) + def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) group = u.atoms[31:39] # molecules 11 @@ -1283,22 +1295,27 @@ def test_center_of_mass_compounds(self, ag, name, compound): @pytest.mark.parametrize('name, compound', (('resids', 'residues'), ('segids', 'segments'))) - def test_center_of_geometry_compounds_pbc(self, ag, name, compound): + @pytest.mark.parametrize('unwrap', (True, False)) + def test_center_of_geometry_compounds_pbc(self, ag, name, compound, + unwrap): ag.dimensions = [50, 50, 50, 90, 90, 90] - ref = [a.center_of_geometry() for a in ag.groupby(name).values()] + ref = [a.center_of_geometry(unwrap=unwrap) + for a in ag.groupby(name).values()] ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32), - ag.dimensions) - cog = ag.center_of_geometry(pbc=True, compound=compound) + ag.dimensions) + cog = ag.center_of_geometry(pbc=True, compound=compound, unwrap=unwrap) assert_almost_equal(cog, ref, decimal=5) @pytest.mark.parametrize('name, compound', (('resids', 'residues'), ('segids', 'segments'))) - def test_center_of_mass_compounds_pbc(self, ag, name, compound): + @pytest.mark.parametrize('unwrap', (True, False)) + def test_center_of_mass_compounds_pbc(self, ag, name, compound, unwrap): ag.dimensions = [50, 50, 50, 90, 90, 90] - ref = [a.center_of_mass() for a in ag.groupby(name).values()] + ref = [a.center_of_mass(unwrap=unwrap) + for a in ag.groupby(name).values()] ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32), - ag.dimensions) - com = ag.center_of_mass(pbc=True, compound=compound) + ag.dimensions) + com = ag.center_of_mass(pbc=True, compound=compound, unwrap=unwrap) assert_almost_equal(com, ref, decimal=5) @pytest.mark.parametrize('name, compound', (('molnums', 'molecules'), @@ -1319,24 +1336,30 @@ def test_center_of_mass_compounds_special(self, ag_molfrg, @pytest.mark.parametrize('name, compound', (('molnums', 'molecules'), ('fragindices', 'fragments'))) + @pytest.mark.parametrize('unwrap', (True, False)) def test_center_of_geometry_compounds_special_pbc(self, ag_molfrg, - name, compound): + name, compound, unwrap): ag_molfrg.dimensions = [50, 50, 50, 90, 90, 90] - ref = [a.center_of_geometry() for a in ag_molfrg.groupby(name).values()] + ref = [a.center_of_geometry(unwrap=unwrap) + for a in ag_molfrg.groupby(name).values()] ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32), - ag_molfrg.dimensions) - cog = ag_molfrg.center_of_geometry(pbc=True, compound=compound) + ag_molfrg.dimensions) + cog = ag_molfrg.center_of_geometry(pbc=True, compound=compound, + unwrap=unwrap) assert_almost_equal(cog, ref, decimal=5) @pytest.mark.parametrize('name, compound', (('molnums', 'molecules'), ('fragindices', 'fragments'))) + @pytest.mark.parametrize('unwrap', (True, False)) def test_center_of_mass_compounds_special_pbc(self, ag_molfrg, - name, compound): + name, compound, unwrap): ag_molfrg.dimensions = [50, 50, 50, 90, 90, 90] - ref = [a.center_of_mass() for a in ag_molfrg.groupby(name).values()] + ref = [a.center_of_mass(unwrap=unwrap) + for a in ag_molfrg.groupby(name).values()] ref = distances.apply_PBC(np.asarray(ref, dtype=np.float32), - ag_molfrg.dimensions) - com = ag_molfrg.center_of_mass(pbc=True, compound=compound) + ag_molfrg.dimensions) + com = ag_molfrg.center_of_mass(pbc=True, compound=compound, + unwrap=unwrap) assert_almost_equal(com, ref, decimal=5) def test_center_wrong_compound(self, ag):