Adds unwrap keyword to shape_parameter

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)

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)
 

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)
 
 

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":