diff --git a/package/CHANGELOG b/package/CHANGELOG index 1471e70456e..83a6dbe1228 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -60,6 +60,7 @@ Enhancements * added unwrap keyword to center_of_geometry (PR #2279) * added unwrap keyword to center_of_mass (PR #2286) * added unwrap keyword to moment_of_inertia (PR #2287) + * added unwrap keyword to asphericity (PR #2290) Changes * added official support for Python 3.7 (PR #1963) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 566d70788c6..22a71302b30 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1028,7 +1028,8 @@ def shape_parameter(group, **kwargs): ('shape_parameter', shape_parameter)) @warn_if_not_unique - def asphericity(group, pbc=None): + @check_pbc_and_unwrap + def asphericity(group, pbc=None, unwrap=None, compound='group'): """Asphericity. See [Dima2004b]_ for background information. @@ -1039,6 +1040,10 @@ def asphericity(group, pbc=None): If ``True``, move all atoms within the primary unit cell before calculation. If ``None`` use value defined in MDAnalysis.core.flags['use_pbc'] + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their centers. + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + Which type of component to keep together during unwrapping. Note ---- @@ -1059,6 +1064,7 @@ def asphericity(group, pbc=None): .. versionadded:: 0.7.7 .. versionchanged:: 0.8 Added *pbc* keyword + .. versionchanged:: 0.20.0 Added *unwrap* and *compound* parameter """ atomgroup = group.atoms @@ -1066,12 +1072,16 @@ def asphericity(group, pbc=None): pbc = flags['use_pbc'] masses = atomgroup.masses + com = atomgroup.center_of_mass(pbc=pbc, unwrap=unwrap, compound=compound) + if compound is not 'group': + com = (com * group.masses[:, None]).sum(axis=0) / group.masses.sum() + if pbc: - recenteredpos = (atomgroup.pack_into_box(inplace=False) - - atomgroup.center_of_mass(pbc=True)) + recenteredpos = (atomgroup.pack_into_box(inplace=False) - com) + elif unwrap: + recenteredpos = (atomgroup.unwrap(inplace=False) - com) else: - recenteredpos = (atomgroup.positions - - atomgroup.center_of_mass(pbc=False)) + recenteredpos = (atomgroup.positions - com) tensor = np.zeros((3, 3)) for x in range(recenteredpos.shape[0]): diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c5070a01022..c17c3013233 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -848,17 +848,17 @@ def ag(self): @pytest.fixture() def ref_noUnwrap_residues(self): return { - 'COG': np.array([[21.356, 28.52 , 36.762], - [32.062, 36.16 , 27.679], + 'COG': np.array([[21.356, 28.52, 36.762], + [32.062, 36.16, 27.679], [27.071, 29.997, 28.506]], dtype=np.float32), 'COM': np.array([[21.286, 28.407, 36.629], [31.931, 35.814, 27.916], - [26.817, 29.41 , 29.05 ]]), + [26.817, 29.41, 29.05]]), 'MOI': np.array([ [7333.79167791, -211.8997285, -721.50785456], [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), - + 'Asph': 0.02060121, } @pytest.fixture() @@ -872,7 +872,8 @@ def ref_Unwrap_residues(self): [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]]) + [ 2938.243, 3306.212, 8990.481]]), + 'Asph': 0.2969491080, } @pytest.fixture() @@ -884,6 +885,7 @@ def ref_noUnwrap(self): [0.0, 0.0, 0.0], [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), + 'Asph': 1.0, } @pytest.fixture() @@ -895,17 +897,20 @@ def ref_Unwrap(self): [0.0, 0.0, 0.0], [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), + 'Asph': 1.0, } def test_default_residues(self, ag, ref_noUnwrap_residues): assert_almost_equal(ag.center_of_geometry(compound='residues'), ref_noUnwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(compound='residues'), ref_noUnwrap_residues['COM'], self.prec) assert_almost_equal(ag.moment_of_inertia(compound='residues'), ref_noUnwrap_residues['MOI'], self.prec) + assert_almost_equal(ag.asphericity(compound='residues'), ref_noUnwrap_residues['Asph'], self.prec) def test_UnWrapFlag_residues(self, ag, ref_Unwrap_residues): assert_almost_equal(ag.center_of_geometry(unwrap=True, compound='residues'), ref_Unwrap_residues['COG'], self.prec) assert_almost_equal(ag.center_of_mass(unwrap=True, compound='residues'), ref_Unwrap_residues['COM'], self.prec) 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_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -916,6 +921,7 @@ def test_default(self, ref_noUnwrap): assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(), ref_noUnwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(), ref_noUnwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) @@ -925,6 +931,8 @@ def test_UnWrapFlag(self, ref_Unwrap): assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec) assert_almost_equal(group.center_of_mass(unwrap=True), ref_Unwrap['COM'], self.prec) assert_almost_equal(group.moment_of_inertia(unwrap=True), ref_Unwrap['MOI'], self.prec) + assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) + class TestPBCFlag(object):