Adds unwrap keyword to moment of inertia

package/CHANGELOG

@@ -59,6 +59,7 @@ Enhancements
   * added unwrap keyword to center (PR #2275)
   * 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)
 
 Changes
   * added official support for Python 3.7 (PR #1963)

package/MDAnalysis/core/__init__.py

@@ -447,7 +447,6 @@ def __doc__(self):
            such as :meth:`MDAnalysis.core.groups.AtomGroup.center_of_mass`
            and :meth:`MDAnalysis.core.groups.AtomGroup.center_of_geometry`!
         """),
-
 ]
 
 #: Global flag registry for :mod:`MDAnalysis.core`.

package/MDAnalysis/core/topologyattrs.py

@@ -870,6 +870,7 @@ def total_mass(group, compound='group'):
         ('total_mass', total_mass))
 
     @warn_if_not_unique
+    @check_pbc_and_unwrap
     def moment_of_inertia(group, **kwargs):
         """Tensor moment of inertia relative to center of mass as 3x3 numpy
         array.
@@ -887,15 +888,22 @@ def moment_of_inertia(group, **kwargs):
 
 
         .. versionchanged:: 0.8 Added *pbc* keyword
+        .. versionchanged:: 0.20.0 Added `unwrap` parameter
 
         """
         atomgroup = group.atoms
         pbc = kwargs.pop('pbc', flags['use_pbc'])
+        unwrap = kwargs.pop('unwrap', False)
+        compound = kwargs.pop('compound', 'group')
+
+        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()
 
-        # Convert to local coordinates
-        com = atomgroup.center_of_mass(pbc=pbc)
         if pbc:
             pos = atomgroup.pack_into_box(inplace=False) - com
+        elif unwrap:
+            pos = atomgroup.unwrap(compound=compound, inplace=False) - com
         else:
             pos = atomgroup.positions - com
 

testsuite/MDAnalysisTests/core/test_atomgroup.py

@@ -838,20 +838,75 @@ class TestUnwrapFlag(object):
 
     prec = 3
 
+    @pytest.fixture()
+    def ag(self):
+        universe = mda.Universe(TRZ_psf, TRZ)
+        group = universe.residues[0:3]
+        group.wrap(inplace=True)
+        return group
+
+    @pytest.fixture()
+    def ref_noUnwrap_residues(self):
+        return {
+            '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 ]]),
+            'MOI': np.array([
+                [7333.79167791, -211.8997285, -721.50785456],
+                [-211.8997285, 7059.07470427, -91.32156884],
+                [-721.50785456, -91.32156884, 6509.31735029]]),
+
+        }
+
+    @pytest.fixture()
+    def ref_Unwrap_residues(self):
+        return {
+            '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]])
+        }
+
     @pytest.fixture()
     def ref_noUnwrap(self):
         return {
             'COG': np.array([5.1, 7.5, 7. ], dtype=np.float32),
             'COM': np.array([6.48785, 7.5, 7.0], dtype=np.float32),
+            'MOI': np.array([
+                [0.0, 0.0, 0.0],
+                [0.0, 98.6542, 0.0],
+                [0.0, 0.0, 98.65421327]]),
         }
 
     @pytest.fixture()
     def ref_Unwrap(self):
         return {
             'COG': np.array([10.1,  7.5,  7. ], dtype=np.float32),
             'COM': np.array([6.8616, 7.5, 7.], dtype=np.float32),
+            'MOI': np.array([
+                [0.0, 0.0, 0.0],
+                [0.0, 132.673, 0.0],
+                [0.0, 0.0, 132.673]]),
         }
 
+    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)
+
+    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)
+
     def test_default(self, ref_noUnwrap):
         u = UnWrapUniverse(is_triclinic=False)
         group = u.atoms[31:39]  # molecules  11
@@ -860,6 +915,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)
 
     def test_UnWrapFlag(self, ref_Unwrap):
         u = UnWrapUniverse(is_triclinic=False)
@@ -868,6 +924,7 @@ def test_UnWrapFlag(self, ref_Unwrap):
         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)
 
 class TestPBCFlag(object):