Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 11 additions & 3 deletions package/MDAnalysis/core/topologyattrs.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
----
Expand Down Expand Up @@ -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()
Expand All @@ -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
-------
Expand All @@ -1119,10 +1128,9 @@ def principal_axes(group, pbc=None):
.. versionchanged:: 0.8 Added *pbc* keyword

"""
atomgroup = group.atoms
Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@richardjgowers @jbarnoud Is there a reason for using atomgroup = group.atoms
in place of directly using group.moment_of_inertia?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here it should be fine since MOI is doing the conversion. In general these methods are attached to all groups, so ResidueGroup too, but they want to operate on atoms. So for a RG, atomgroup = group.atoms is converting to an AtomGroup within the method.

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]
Expand Down
22 changes: 21 additions & 1 deletion testsuite/MDAnalysisTests/core/test_atomgroup.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]])

}

Expand All @@ -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()
Expand All @@ -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()
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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):

Expand Down