diff --git a/package/CHANGELOG b/package/CHANGELOG index 83a6dbe1228..d054a1b40b3 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -62,6 +62,7 @@ Enhancements * 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) * stopped official support of Python 3.4 (#2066, PR #2174) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 22a71302b30..7b5654ebbf2 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -945,6 +945,10 @@ def radius_of_gyration(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ---- @@ -974,6 +978,7 @@ def radius_of_gyration(group, **kwargs): ('radius_of_gyration', radius_of_gyration)) @warn_if_not_unique + @check_pbc_and_unwrap def shape_parameter(group, **kwargs): """Shape parameter. @@ -984,6 +989,10 @@ def shape_parameter(group, **kwargs): pbc : bool, optional If ``True``, move all atoms within the primary unit cell before calculation. [``False``] + 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 ---- @@ -1004,21 +1013,28 @@ def shape_parameter(group, **kwargs): .. versionchanged:: 0.8 Added *pbc* keyword """ - atomgroup = group.atoms pbc = kwargs.pop('pbc', flags['use_pbc']) - masses = atomgroup.masses + masses = group.atoms.masses + unwrap = kwargs.pop('unwrap', False) + compound = kwargs.pop('compound', 'group') + + com = group.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() - com = atomgroup.center_of_mass(pbc=pbc) if pbc: - recenteredpos = atomgroup.pack_into_box(inplace=False) - com + recenteredpos = group.pack_into_box(inplace=False) - com + elif unwrap: + recenteredpos = group.unwrap(compound=compound) - com else: - recenteredpos = atomgroup.positions - com + recenteredpos = group.atoms.positions - com + tensor = np.zeros((3, 3)) for x in range(recenteredpos.shape[0]): tensor += masses[x] * np.outer(recenteredpos[x, :], recenteredpos[x, :]) - tensor /= atomgroup.total_mass() + tensor /= group.total_mass() eig_vals = np.linalg.eigvalsh(tensor) shape = 27.0 * np.prod(eig_vals - np.mean(eig_vals)) / np.power(np.sum(eig_vals), 3) diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py index c17c3013233..45fe2c8c806 100644 --- a/testsuite/MDAnalysisTests/core/test_atomgroup.py +++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py @@ -858,7 +858,9 @@ def ref_noUnwrap_residues(self): [7333.79167791, -211.8997285, -721.50785456], [-211.8997285, 7059.07470427, -91.32156884], [-721.50785456, -91.32156884, 6509.31735029]]), + 'Shape': 0.00173905, 'Asph': 0.02060121, + 'Shape': 0.00173905, } @pytest.fixture() @@ -873,7 +875,9 @@ def ref_Unwrap_residues(self): 'MOI': np.array([[16747.486, -1330.489, 2938.243], [-1330.489, 19315.253, 3306.212], [ 2938.243, 3306.212, 8990.481]]), + 'Shape': 0.27397619, 'Asph': 0.2969491080, + 'Shape': 0.27397619, } @pytest.fixture() @@ -885,7 +889,9 @@ def ref_noUnwrap(self): [0.0, 0.0, 0.0], [0.0, 98.6542, 0.0], [0.0, 0.0, 98.65421327]]), + 'Shape': 2.0, 'Asph': 1.0, + 'Shape': 2.0, } @pytest.fixture() @@ -897,20 +903,23 @@ def ref_Unwrap(self): [0.0, 0.0, 0.0], [0.0, 132.673, 0.0], [0.0, 0.0, 132.673]]), + 'Shape': 2.0, 'Asph': 1.0, + 'Shape': 2.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.shape_parameter(compound='residues'), ref_noUnwrap_residues['Shape'], 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) + assert_almost_equal(ag.shape_parameter(unwrap=True, compound='residues'), ref_Unwrap_residues['Shape'], self.prec) def test_default(self, ref_noUnwrap): u = UnWrapUniverse(is_triclinic=False) @@ -921,16 +930,19 @@ 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.shape_parameter(), ref_noUnwrap['Shape'], self.prec) assert_almost_equal(group.asphericity(), ref_noUnwrap['Asph'], self.prec) def test_UnWrapFlag(self, ref_Unwrap): u = UnWrapUniverse(is_triclinic=False) group = u.atoms[31:39] # molecules 11 + # Changing masses for center_of_mass group.masses = [100.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0] 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.shape_parameter(unwrap=True), ref_Unwrap['Shape'], self.prec) assert_almost_equal(group.asphericity(unwrap=True), ref_Unwrap['Asph'], self.prec) diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py index 653c743d19c..7aea81ad6de 100644 --- a/testsuite/MDAnalysisTests/core/util.py +++ b/testsuite/MDAnalysisTests/core/util.py @@ -585,7 +585,7 @@ def center(self, compound): for base in range(3, 15, 3): loc_center = np.mean(relpos[base:base + 3, :], axis=0) - center_pos[pos,:] = loc_center + center_pos[pos,:] = loc_center pos+=1 if compound=="residues":