From 4819464618de18bc0460eb11bb6dfba132193f48 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Mon, 17 Jun 2019 13:15:32 +0530
Subject: [PATCH 1/7] Adds unwrap to center

---
 package/MDAnalysis/core/groups.py | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index e9dcfecc458..4ee7af04608 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -649,7 +649,7 @@ def isunique(self):
         return not np.count_nonzero(mask)
 
     @warn_if_not_unique
-    def center(self, weights, pbc=None, compound='group'):
+    def center(self, weights, pbc=None, compound='group', unwrap=False):
         """Weighted center of (compounds of) the group
 
         Computes the weighted center of :class:`Atoms<Atom>` in the group.
@@ -679,6 +679,11 @@ def center(self, weights, pbc=None, compound='group'):
             will be returned as an array of position vectors, i.e. a 2d array.
             Note that, in any case, *only* the positions of :class:`Atoms<Atom>`
             *belonging to the group* will be taken into account.
+        unwrap : bool or None, optional
+            If ``True`` and `compound` is ``'group'``, the atoms will be unwrapped
+            before calculations. If ``True`` and `compound` is
+            ``'segments'`` or ``'residues'``, all molecules will be unwrapped before
+             calculation to keep the compounds intact.
 
         Returns
         -------
@@ -734,6 +739,8 @@ def center(self, weights, pbc=None, compound='group'):
 
         comp = compound.lower()
         if comp == 'group':
+            if unwrap:
+                coords = atoms.unwrap(inplace=False)
             if pbc:
                 coords = atoms.pack_into_box(inplace=False)
             else:
@@ -769,11 +776,19 @@ def center(self, weights, pbc=None, compound='group'):
                              " one of 'group', 'residues', 'segments', "
                              "'molecules', or 'fragments'.".format(compound))
 
+
+
         # Sort positions and weights by compound index and promote to dtype if
         # required:
         sort_indices = np.argsort(compound_indices)
         compound_indices = compound_indices[sort_indices]
-        coords = atoms.positions[sort_indices]
+
+        # Unwrap Atoms
+        if unwrap:
+            coords = atoms.unwrap(inplace=False)
+        else:
+            coords = atoms.positions[sort_indices]
+
         if weights is None:
             coords = coords.astype(dtype, copy=False)
         else:

From 300c955a98b331186f6ccb815e7b33f410a78290 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Tue, 18 Jun 2019 23:26:28 +0530
Subject: [PATCH 2/7] Adds tests

---
 package/MDAnalysis/core/groups.py             |  5 +-
 .../MDAnalysisTests/core/test_atomgroup.py    | 27 +++++++-
 testsuite/MDAnalysisTests/core/util.py        | 65 +++++++++++++++++++
 3 files changed, 93 insertions(+), 4 deletions(-)

diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index 4ee7af04608..45fc061a2c6 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -740,7 +740,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
         comp = compound.lower()
         if comp == 'group':
             if unwrap:
-                coords = atoms.unwrap(inplace=False)
+                coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
             if pbc:
                 coords = atoms.pack_into_box(inplace=False)
             else:
@@ -785,10 +785,9 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
 
         # Unwrap Atoms
         if unwrap:
-            coords = atoms.unwrap(inplace=False)
+            coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
         else:
             coords = atoms.positions[sort_indices]
-
         if weights is None:
             coords = coords.astype(dtype, copy=False)
         else:
diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py
index 617cdb7b959..6da44f5b6d5 100644
--- a/testsuite/MDAnalysisTests/core/test_atomgroup.py
+++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py
@@ -52,7 +52,7 @@
     GRO
 )
 from MDAnalysisTests import make_Universe, no_deprecated_call
-
+from MDAnalysisTests.core.util import UnWrapUniverse
 import pytest
 
 
@@ -518,6 +518,31 @@ def test_center_wrong_shape(self, ag):
         with pytest.raises(ValueError):
             ag.center(weights)
 
+    @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments'))
+    @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues',
+                                          'group', 'segments'))
+    @pytest.mark.parametrize('is_triclinic', (False, True))
+    def test_center_unwrap(self, level, compound, is_triclinic):
+        u = UnWrapUniverse(is_triclinic=is_triclinic)
+        # select group appropriate for compound:
+        if compound == 'group':
+            group = u.atoms[39:47] # molecule 12
+        elif compound == 'segments':
+            group = u.atoms[23:47] # molecules 10, 11, 12
+        else:
+            group = u.atoms
+        # select topology level:
+        if level == 'residues':
+            group = group.residues
+        elif level == 'segments':
+            group = group.segments
+
+        # get the expected results
+        center = group.center(weights=None, compound=compound, unwrap=True)
+
+        ref_center = u.center(compound=compound)
+        assert_almost_equal(ref_center, center, decimal=4)
+
 
 class TestSplit(object):
 
diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py
index 664dbabbef9..ad2b911f38b 100644
--- a/testsuite/MDAnalysisTests/core/util.py
+++ b/testsuite/MDAnalysisTests/core/util.py
@@ -333,6 +333,7 @@ def __new__(cls, have_bonds=True, have_masses=True, have_molnums=True,
         # bind custom methods to universe:
         u.unwrapped_coords = cls.unwrapped_coords.__get__(u)
         u.wrapped_coords = cls.wrapped_coords.__get__(u)
+        u.center = cls.center.__get__(u)
 
         return u
 
@@ -545,3 +546,67 @@ def wrapped_coords(self, compound, center):
         positions = relpos * np.array([a, a, a])
 
         return positions.astype(np.float32)
+
+    def center(self, compound):
+        """Returns centers which correspond to the unwrapped system.
+
+        Parameters
+        ----------
+        compound : {'atoms', 'group', 'segments', 'residues', 'molecules', \
+                    'fragments'}
+            Which type of component is unwrapped. Note that for ``'group'``,
+            the result will only be correct *if the group is the entire system*.
+
+        Note
+        ----
+        This function assumes that all atom masses are equal. Therefore, the
+        returned coordinates for ``center='com'`` and ``center='cog'`` are
+        identical.
+        """
+
+        relpos = self.unwrapped_coords(compound, reference=None)
+
+        comp = compound.lower()
+        if comp not in ['group', 'segments', 'residues', 'molecules',
+                        'fragments']:
+            raise ValueError("Unknown unwrap compound: {}".format(compound))
+
+        pos = 0
+
+        if compound=="residues":
+            center_pos = np.zeros((15, 3), dtype=np.float32)
+        else:
+            center_pos = np.zeros((12, 3), dtype=np.float32)
+
+        for base in range(3):
+            loc_centre = relpos[base, :]
+            center_pos[pos,:] = loc_centre
+            pos+=1
+
+        for base in range(3, 15, 3):
+            loc_centre = np.mean(relpos[base:base + 3, :], axis=0)
+            center_pos[pos,:] = loc_centre
+            pos+=1
+
+        if compound=="residues":
+            for base in range(15, 47, 4):
+                loc_centre = np.mean(relpos[base:base + 4, :], axis=0)
+                center_pos[pos,:] = loc_centre
+                pos+=1
+        else:
+            for base in range(15, 23, 4):
+                loc_centre = np.mean(relpos[base:base + 4, :], axis=0)
+                center_pos[pos,:] = loc_centre
+                pos+=1
+            for base in range(23, 47, 8):
+                loc_centre = np.mean(relpos[base:base + 8, :], axis=0)
+                center_pos[pos,:] = loc_centre
+                pos+=1
+
+        if compound == "group":
+            center_pos = center_pos[11]
+        elif compound == "segments":
+            center_pos = center_pos[9:]
+
+        return center_pos
+

From 6f21b3143c60199e95fc50bb7d95de878aeabf33 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Wed, 19 Jun 2019 10:41:37 +0530
Subject: [PATCH 3/7] Review changes

---
 package/MDAnalysis/core/groups.py             | 12 +++++------
 .../MDAnalysisTests/core/test_atomgroup.py    |  7 +++++++
 testsuite/MDAnalysisTests/core/util.py        | 20 +++++++++----------
 3 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index 45fc061a2c6..1e8ec7df630 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -679,11 +679,8 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
             will be returned as an array of position vectors, i.e. a 2d array.
             Note that, in any case, *only* the positions of :class:`Atoms<Atom>`
             *belonging to the group* will be taken into account.
-        unwrap : bool or None, optional
-            If ``True`` and `compound` is ``'group'``, the atoms will be unwrapped
-            before calculations. If ``True`` and `compound` is
-            ``'segments'`` or ``'residues'``, all molecules will be unwrapped before
-             calculation to keep the compounds intact.
+        unwrap : bool, optional
+            If ``True``, compounds will be unwrapped before computing their centers.
 
         Returns
         -------
@@ -739,6 +736,9 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
 
         comp = compound.lower()
         if comp == 'group':
+            if unwrap and pbc:
+                raise ValueError("'unwrap' and 'pbc' cannot be true at the same time "
+                                 "for compound='group")
             if unwrap:
                 coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
             if pbc:
@@ -776,8 +776,6 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
                              " one of 'group', 'residues', 'segments', "
                              "'molecules', or 'fragments'.".format(compound))
 
-
-
         # Sort positions and weights by compound index and promote to dtype if
         # required:
         sort_indices = np.argsort(compound_indices)
diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py
index 6da44f5b6d5..c0e9a1b58cc 100644
--- a/testsuite/MDAnalysisTests/core/test_atomgroup.py
+++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py
@@ -543,6 +543,13 @@ def test_center_unwrap(self, level, compound, is_triclinic):
         ref_center = u.center(compound=compound)
         assert_almost_equal(ref_center, center, decimal=4)
 
+    def test_center_unwrap_pbc_true(self):
+        u = UnWrapUniverse(is_triclinic=False)
+        # select group appropriate for compound:
+        group = u.atoms[39:47]  # molecule 12
+        with pytest.raises(ValueError):
+            group.center(weights=None, compound="group", unwrap=True, pbc=True)
+
 
 class TestSplit(object):
 
diff --git a/testsuite/MDAnalysisTests/core/util.py b/testsuite/MDAnalysisTests/core/util.py
index ad2b911f38b..653c743d19c 100644
--- a/testsuite/MDAnalysisTests/core/util.py
+++ b/testsuite/MDAnalysisTests/core/util.py
@@ -579,28 +579,28 @@ def center(self, compound):
             center_pos = np.zeros((12, 3), dtype=np.float32)
 
         for base in range(3):
-            loc_centre = relpos[base, :]
-            center_pos[pos,:] = loc_centre
+            loc_center = relpos[base, :]
+            center_pos[pos,:] = loc_center
             pos+=1
 
         for base in range(3, 15, 3):
-            loc_centre = np.mean(relpos[base:base + 3, :], axis=0)
-            center_pos[pos,:] = loc_centre
+            loc_center = np.mean(relpos[base:base + 3, :], axis=0)
+            center_pos[pos,:] = loc_center
             pos+=1
 
         if compound=="residues":
             for base in range(15, 47, 4):
-                loc_centre = np.mean(relpos[base:base + 4, :], axis=0)
-                center_pos[pos,:] = loc_centre
+                loc_center = np.mean(relpos[base:base + 4, :], axis=0)
+                center_pos[pos,:] = loc_center
                 pos+=1
         else:
             for base in range(15, 23, 4):
-                loc_centre = np.mean(relpos[base:base + 4, :], axis=0)
-                center_pos[pos,:] = loc_centre
+                loc_center = np.mean(relpos[base:base + 4, :], axis=0)
+                center_pos[pos,:] = loc_center
                 pos+=1
             for base in range(23, 47, 8):
-                loc_centre = np.mean(relpos[base:base + 8, :], axis=0)
-                center_pos[pos,:] = loc_centre
+                loc_center = np.mean(relpos[base:base + 8, :], axis=0)
+                center_pos[pos,:] = loc_center
                 pos+=1
 
         if compound == "group":

From 55d0a84392e2e31e129e640e6a01d63aba00f11c Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Wed, 19 Jun 2019 18:26:59 +0530
Subject: [PATCH 4/7] add unwrap to 'cog'

---
 package/MDAnalysis/core/groups.py             | 11 +++++---
 .../MDAnalysisTests/core/test_atomgroup.py    | 26 +++++++++++++++++++
 2 files changed, 33 insertions(+), 4 deletions(-)

diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index 1e8ec7df630..1b9c06a5955 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -739,12 +739,13 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
             if unwrap and pbc:
                 raise ValueError("'unwrap' and 'pbc' cannot be true at the same time "
                                  "for compound='group")
-            if unwrap:
-                coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
             if pbc:
                 coords = atoms.pack_into_box(inplace=False)
             else:
                 coords = atoms.positions
+            if unwrap:
+                coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
+
             # If there's no atom, return its (empty) coordinates unchanged.
             if len(atoms) == 0:
                 return coords
@@ -816,7 +817,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
         return centers
 
     @warn_if_not_unique
-    def center_of_geometry(self, pbc=None, compound='group'):
+    def center_of_geometry(self, pbc=None, compound='group', unwrap=False):
         """Center of geometry of (compounds of) the group.
 
         Computes the center of geometry (a.k.a. centroid) of
@@ -840,6 +841,8 @@ def center_of_geometry(self, pbc=None, compound='group'):
             will be returned as an array of position vectors, i.e. a 2d array.
             Note that, in any case, *only* the positions of :class:`Atoms<Atom>`
             *belonging to the group* will be taken into account.
+        unwrap : bool, optional
+            If ``True``, compounds will be unwrapped before computing their centers.
 
         Returns
         -------
@@ -862,7 +865,7 @@ def center_of_geometry(self, pbc=None, compound='group'):
         .. versionchanged:: 0.20.0 Added ``'molecules'`` and ``'fragments'``
             compounds
         """
-        return self.center(None, pbc=pbc, compound=compound)
+        return self.center(None, pbc=pbc, compound=compound, unwrap=unwrap)
 
     centroid = center_of_geometry
 
diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py
index c0e9a1b58cc..cb4e211f35b 100644
--- a/testsuite/MDAnalysisTests/core/test_atomgroup.py
+++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py
@@ -834,6 +834,32 @@ def test_chi1_nodep(self, PSFDCD):
             sel = PSFDCD.segments[0].residues[12].chi1_selection()  # LYS
 
 
+class TestUnwrapFlag(object):
+
+    prec = 3
+
+    @pytest.fixture()
+    def ref_noUnwrap(self):
+        return {
+            'COG': np.array([5.1, 7.5, 7. ], dtype=np.float32),
+        }
+
+    @pytest.fixture()
+    def ref_Unwrap(self):
+        return {
+            'COG': np.array([10.1,  7.5,  7. ], dtype=np.float32),
+        }
+
+    def test_default(self, ref_noUnwrap):
+        u = UnWrapUniverse(is_triclinic=False)
+        group = u.atoms[31:39]  # molecules  11
+        assert_almost_equal(group.center_of_geometry(), ref_noUnwrap['COG'], self.prec)
+
+    def test_UnWrapFlag(self, ref_Unwrap):
+        u = UnWrapUniverse(is_triclinic=False)
+        group = u.atoms[31:39]  # molecules  11
+        assert_almost_equal(group.center_of_geometry(unwrap=True), ref_Unwrap['COG'], self.prec)
+
 class TestPBCFlag(object):
 
     prec = 3

From 6ae780bea1b64b3542104c0521f1c19809047fc2 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Wed, 3 Jul 2019 22:13:26 +0530
Subject: [PATCH 5/7] check for pbc and unwrap

---
 package/CHANGELOG                 | 1 +
 package/MDAnalysis/core/groups.py | 2 ++
 2 files changed, 3 insertions(+)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index 997e420616f..65684073954 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -57,6 +57,7 @@ Enhancements
   * survival probability additions: residues, intermittency, step with performance,
     (PR #2226)   
   * added unwrap keyword to center (PR #2275)
+  * added unwrap keyword to center_of_geometry (PR #2279)
 
 Changes
   * added official support for Python 3.7 (PR #1963)
diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index 129f2cb5de6..0b338132ab7 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -833,6 +833,7 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
         return centers
 
     @warn_if_not_unique
+    @check_pbc_and_unwrap
     def center_of_geometry(self, pbc=None, compound='group', unwrap=False):
         """Center of geometry of (compounds of) the group.
 
@@ -880,6 +881,7 @@ def center_of_geometry(self, pbc=None, compound='group', unwrap=False):
         .. versionchanged:: 0.19.0 Added `compound` parameter
         .. versionchanged:: 0.20.0 Added ``'molecules'`` and ``'fragments'``
             compounds
+        .. versionchanged:: 0.20.0 Added `unwrap` parameter
         """
         return self.center(None, pbc=pbc, compound=compound, unwrap=unwrap)
 

From 6d33fe940b4069f1f930064893fd53fd221cd792 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Thu, 4 Jul 2019 03:04:34 +0530
Subject: [PATCH 6/7] Removes unnecessary lines

---
 package/MDAnalysis/core/groups.py | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py
index 0b338132ab7..ef80087d794 100644
--- a/package/MDAnalysis/core/groups.py
+++ b/package/MDAnalysis/core/groups.py
@@ -750,18 +750,12 @@ def center(self, weights, pbc=None, compound='group', unwrap=False):
 
         comp = compound.lower()
         if comp == 'group':
-            if unwrap and pbc:
-                raise ValueError("'unwrap' and 'pbc' cannot be true at the same time "
-                                 "for compound='group")
             if pbc:
                 coords = atoms.pack_into_box(inplace=False)
             elif unwrap:
                 coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
             else:
                 coords = atoms.positions
-            if unwrap:
-                coords = atoms.unwrap(compound=comp, reference=None, inplace=False)
-
             # If there's no atom, return its (empty) coordinates unchanged.
             if len(atoms) == 0:
                 return coords

From a6bdacdbedf304f442f7da6772248fd87ad44d46 Mon Sep 17 00:00:00 2001
From: Ninad Bhat <bhat.ninadmb@gmail.com>
Date: Tue, 9 Jul 2019 18:01:09 +0530
Subject: [PATCH 7/7] adds unwrap to principal_axes

---
 package/MDAnalysis/core/topologyattrs.py      | 14 +++++++++---
 .../MDAnalysisTests/core/test_atomgroup.py    | 22 ++++++++++++++++++-
 2 files changed, 32 insertions(+), 4 deletions(-)

diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py
index 566d70788c6..2d6a1dcf295 100644
--- a/package/MDAnalysis/core/topologyattrs.py
+++ b/package/MDAnalysis/core/topologyattrs.py
@@ -880,6 +880,10 @@ def moment_of_inertia(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
         ----
@@ -1089,7 +1093,8 @@ def asphericity(group, pbc=None):
         ('asphericity', asphericity))
 
     @warn_if_not_unique
-    def principal_axes(group, pbc=None):
+    @check_pbc_and_unwrap
+    def principal_axes(group, pbc=None, unwrap=False, compound='group'):
         """Calculate the principal axes from the moment of inertia.
 
         e1,e2,e3 = AtomGroup.principal_axes()
@@ -1103,6 +1108,10 @@ def principal_axes(group, pbc=None):
         pbc : bool, optional
             If ``True``, move all atoms within the primary unit cell before
             calculation. If ``None`` use value defined in setup flags.
+        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.
 
         Returns
         -------
@@ -1119,10 +1128,9 @@ def principal_axes(group, pbc=None):
         .. versionchanged:: 0.8 Added *pbc* keyword
 
         """
-        atomgroup = group.atoms
         if pbc is None:
             pbc = flags['use_pbc']
-        e_val, e_vec = np.linalg.eig(atomgroup.moment_of_inertia(pbc=pbc))
+        e_val, e_vec = np.linalg.eig(group.moment_of_inertia(pbc=pbc, unwrap=unwrap, compound=compound))
 
         # Sort
         indices = np.argsort(e_val)[::-1]
diff --git a/testsuite/MDAnalysisTests/core/test_atomgroup.py b/testsuite/MDAnalysisTests/core/test_atomgroup.py
index c5070a01022..d6aa43ef315 100644
--- a/testsuite/MDAnalysisTests/core/test_atomgroup.py
+++ b/testsuite/MDAnalysisTests/core/test_atomgroup.py
@@ -858,6 +858,10 @@ def ref_noUnwrap_residues(self):
                 [7333.79167791, -211.8997285, -721.50785456],
                 [-211.8997285, 7059.07470427, -91.32156884],
                 [-721.50785456, -91.32156884, 6509.31735029]]),
+            'PAxes': np.array([
+                [-0.85911708, 0.19258726, 0.4741603],
+                [-0.07520116, -0.96394227, 0.25526473],
+                [-0.50622389, -0.18364489, -0.84262206]])
 
         }
 
@@ -872,7 +876,11 @@ 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]]),
+            'PAxes': np.array(
+                [[-0.15999399, 0.9581568, 0.23735515],
+                 [-0.93680726, -0.0715994, -0.34244077],
+                 [-0.31111746, -0.27714449, 0.90906372]]),
         }
 
     @pytest.fixture()
@@ -884,6 +892,10 @@ def ref_noUnwrap(self):
                 [0.0, 0.0, 0.0],
                 [0.0, 98.6542, 0.0],
                 [0.0, 0.0, 98.65421327]]),
+            'PAxes': np.array(
+                [[0., 0., 1.],
+                 [0., 1., 0.],
+                 [1., 0., 0.]]),
         }
 
     @pytest.fixture()
@@ -895,17 +907,23 @@ def ref_Unwrap(self):
                 [0.0, 0.0, 0.0],
                 [0.0, 132.673, 0.0],
                 [0.0, 0.0, 132.673]]),
+            'PAxes': np.array(
+                [[0., 0., 1.],
+                 [0., 1., 0.],
+                 [1., 0., 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.principal_axes(compound='residues'), ref_noUnwrap_residues['PAxes'], 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.principal_axes(unwrap=True, compound='residues'), ref_Unwrap_residues['PAxes'], self.prec)
 
     def test_default(self, ref_noUnwrap):
         u = UnWrapUniverse(is_triclinic=False)
@@ -916,6 +934,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.principal_axes(), ref_noUnwrap['PAxes'], self.prec)
 
     def test_UnWrapFlag(self, ref_Unwrap):
         u = UnWrapUniverse(is_triclinic=False)
@@ -925,6 +944,7 @@ 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.principal_axes(unwrap=True), ref_Unwrap['PAxes'], self.prec)
 
 class TestPBCFlag(object):