diff --git a/package/CHANGELOG b/package/CHANGELOG index f243351b7d7..83f2aab4358 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -15,11 +15,13 @@ The rules for this file: ------------------------------------------------------------------------------ ??/??/?? tylerjereddy, richardjgowers, IAlibay, hmacdope, orbeckst, cbouy, lilyminium, daveminh, jbarnoud, yuxuanzhuang, VOD555, ianmkenney, - calcraven，xiki-tempula, mieczyslaw + calcraven，xiki-tempula, mieczyslaw, manuel.nuno.melo * 2.0.0 Fixes + * AtomGroup.center now works correctly for compounds + unwrapping + (Issue #2984) * Avoid using deprecated array indexing in topology attributes (Issue #2990, PR #2991) * ParmedParser no longer guesses elements if they are not recognised, instead diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 0264617c469..274ace30947 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -830,12 +830,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 520519046e3..5d6421d1856 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -991,6 +991,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 { @@ -1013,12 +1019,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, } @@ -1058,6 +1064,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 @@ -1281,22 +1293,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'), @@ -1317,24 +1334,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):