Add calculations for shape metrics (mdtraj#1471)

* Craeted shape.py and did some refactoring

* Added shape metrics, added tests, mentioned in docs

* changed documentation

docs/analysis.rst

@@ -48,6 +48,11 @@ Shape Metrics
     :toctree: api/generated/
 
     compute_gyration_tensor
+    principal_moments
+    asphericity
+    acylindricity
+    relative_shape_antisotropy
+
 
 
 Surface Area, Radius of Gyration and Inertia

mdtraj/geometry/__init__.py

@@ -31,7 +31,8 @@
            'wernet_nilsson', 'compute_dssp', 'compute_neighbors',
            'compute_neighborlist', 'compute_rdf', 'compute_nematic_order',
            'compute_inertia_tensor', 'compute_gyration_tensor', 'find_closest_contact',
-           'compute_directors',
+           'compute_directors', 'principal_moments', 'asphericity', 'acylindricity', 
+           'relative_shape_antisotropy',
 
            # from thermodynamic_properties
            'dipole_moments', 'static_dielectric', 'isothermal_compressability_kappa_T',
@@ -52,3 +53,4 @@
 from .thermodynamic_properties import *
 from .rdf import *
 from .order import *
+from .shape import *

mdtraj/geometry/order.py

@@ -32,7 +32,7 @@
 from mdtraj.utils.six import string_types
 
 
-__all__ = ['compute_nematic_order', 'compute_inertia_tensor', 'compute_directors', 'compute_gyration_tensor']
+__all__ = ['compute_nematic_order', 'compute_inertia_tensor', 'compute_directors']
 
 
 def compute_nematic_order(traj, indices='chains'):
@@ -187,36 +187,6 @@ def compute_inertia_tensor(traj):
     A = np.einsum("i, kij->k", masses, xyz ** 2).reshape(traj.n_frames, 1, 1)
     B = np.einsum("ij..., ...jk->...ki", masses[:, np.newaxis] * xyz.T, xyz)
     return A * eyes - B
-
-
-def compute_gyration_tensor(traj):
-    """Compute the gyration tensor of a trajectory.
-
-    For each frame,
-
-    .. math::
-
-        S_{xy} = sum_{i_atoms} r^{i}_x r^{i}_y
-
-    Parameters
-    ----------
-    traj : Trajectory
-        Trajectory to compute gyration tensor of.
-
-    Returns
-    -------
-    S_xy:  np.ndarray, shape=(traj.n_frames, 3, 3), dtype=float64
-        Gyration tensors for each frame.
-    
-    References
-    ----------
-    .. [1] https://isg.nist.gov/deepzoomweb/measurement3Ddata_help#shape-metrics-formulas
-
-    """
-    center_of_geom = np.expand_dims(compute_center_of_geometry(traj), axis=1)
-    xyz = traj.xyz - center_of_geom
-    return np.einsum('...ji,...jk->...ik', xyz, xyz) / traj.n_atoms
-
 
 
 def _get_indices(traj, indices):
@@ -365,26 +335,4 @@ def _compute_inertia_tensor_slow(traj):
         I_ab[n, 1, 0] = I_ab[n, 0, 1]
         I_ab[n, 2, 0] = I_ab[n, 0, 2]
         I_ab[n, 2, 1] = I_ab[n, 1, 2]
-    return I_ab
-
-def _compute_gyration_tensor_slow(traj):
-    """Compute the gyration tensor of a trajectory. """
-    xyz = traj.xyz
-    center_of_geom = np.expand_dims(compute_center_of_geometry(traj), axis=1)
-    centered_xyz = xyz - center_of_geom
-
-    S_nm = np.zeros(shape=(traj.n_frames, 3, 3), dtype=np.float64)
-    for n, xyz in enumerate(centered_xyz):
-        N = xyz.shape[0]
-        for r in xyz:
-            S_nm[n, 0, 0] += r[0] * r[0]
-            S_nm[n, 1, 1] += r[1] * r[1]
-            S_nm[n, 2, 2] += r[2] * r[2]
-            S_nm[n, 0, 1] += r[0] * r[1]
-            S_nm[n, 0, 2] += r[0] * r[2]
-            S_nm[n, 1, 2] += r[1] * r[2]
-        S_nm[n, 1, 0] = S_nm[n, 0, 1]
-        S_nm[n, 2, 0] = S_nm[n, 0, 2]
-        S_nm[n, 2, 1] = S_nm[n, 1, 2]
-        S_nm[n, :, :] /= N
-    return S_nm
+    return I_ab

mdtraj/geometry/shape.py

@@ -0,0 +1,174 @@
+##############################################################################
+# MDTraj: A Python Library for Loading, Saving, and Manipulating
+#         Molecular Dynamics Trajectories.
+# Copyright 2012-2015 Stanford University and the Authors
+#
+# Authors: Christoph Klein
+# Contributors: Tim Moore
+#
+# MDTraj is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as
+# published by the Free Software Foundation, either version 2.1
+# of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with MDTraj. If not, see <http://www.gnu.org/licenses/>.
+##############################################################################
+
+from __future__ import print_function, division
+
+import numpy as np
+
+from mdtraj.geometry.distance import compute_center_of_geometry
+
+__all__ = ['compute_gyration_tensor', 'principal_moments', 'asphericity', 'acylindricity', 'relative_shape_antisotropy']
+
+
+def compute_gyration_tensor(traj):
+    """Compute the gyration tensor of a trajectory.
+
+    For every frame,
+
+    .. math::
+
+        S_{xy} = sum_{i_atoms} r^{i}_x r^{i}_y
+
+    Parameters
+    ----------
+    traj : Trajectory
+        Trajectory to compute gyration tensor of.
+
+    Returns
+    -------
+    S_xy:  np.ndarray, shape=(traj.n_frames, 3, 3), dtype=float64
+        Gyration tensors for each frame.
+    
+    References
+    ----------
+    .. [1] https://isg.nist.gov/deepzoomweb/measurement3Ddata_help#shape-metrics-formulas
+
+    """
+    center_of_geom = np.expand_dims(compute_center_of_geometry(traj), axis=1)
+    xyz = traj.xyz - center_of_geom
+    return np.einsum('...ji,...jk->...ik', xyz, xyz) / traj.n_atoms
+
+def principal_moments(traj):
+    """Compute the principal moments of a trajectory.
+
+    For each frame calculate the eigenvalues of the gyration tensor.
+
+    Parameters
+    ----------
+    traj : Trajectory
+        Trajectory to compute gyration tensor of.
+
+    Returns
+    -------
+    evecs:  np.ndarray, shape=(traj.n_frames, 3), dtype=float64
+        Principal moments for each frame in ascending order.
+
+    """
+    gyration_tensor = compute_gyration_tensor(traj)
+    return np.linalg.eigvalsh(gyration_tensor)
+
+def asphericity(traj):
+    """Compute the asphericity of a trajectory.
+
+    For each frame compute the principal moments then,
+
+    .. math::
+
+        b = \frac{1}{2}(\lambda_1^2 + \lambda_2^2)
+
+
+    Parameters
+    ----------
+    traj : Trajectory
+        Trajectory to compute gyration tensor of.
+
+    Returns
+    -------
+    b:  np.ndarray, shape=(traj.n_frames, 1), dtype=float64
+        Asphericity of each frame of the trajectory.
+
+    """
+    pm = principal_moments(traj)
+    b = pm[:,2] - (pm[:,0] + pm[:,1]) / 2.0
+    return b
+
+def acylindricity(traj):
+    """Compute the acylindricity of a trajectory.
+
+    For each frame compute the principal moments then,
+
+    .. math::
+
+        c = \lambda_2^2 - \lambda_1^2
+
+
+    Parameters
+    ----------
+    traj : Trajectory
+        Trajectory to compute gyration tensor of.
+
+    Returns
+    -------
+    c:  np.ndarray, shape=(traj.n_frames, 1), dtype=float64
+        Acylindricity of each frame of the trajectory.
+
+    """
+    pm = principal_moments(traj)
+    c = pm[:,1] - pm[:,0]
+    return c
+
+def relative_shape_antisotropy(traj):
+    """Compute the acylindricity of a trajectory.
+
+    For each frame compute the principal moments then,
+
+    .. math::
+
+        \kappa^2 = \frac{3}{2}\frac{\lambda_1^4 + \lambda_2^4 + \lambda_3^4}{(\lambda_1^2 + \lambda_2^2 + \lambda_3^2)^2} - \frac{1}{2}
+
+
+    Parameters
+    ----------
+    traj : Trajectory
+        Trajectory to compute gyration tensor of.
+
+    Returns
+    -------
+    c:  np.ndarray, shape=(traj.n_frames, 1), dtype=float64
+        Acylindricity of each frame of the trajectory.
+
+    """
+    pm = principal_moments(traj)
+    kappa2 = 1.5 * np.square(pm).sum(axis=1) / np.square(pm.sum(axis=1)) - 0.5
+    return kappa2
+
+def _compute_gyration_tensor_slow(traj):
+    """Compute the gyration tensor of a trajectory. """
+    xyz = traj.xyz
+    center_of_geom = np.expand_dims(compute_center_of_geometry(traj), axis=1)
+    centered_xyz = xyz - center_of_geom
+
+    S_nm = np.zeros(shape=(traj.n_frames, 3, 3), dtype=np.float64)
+    for n, xyz in enumerate(centered_xyz):
+        N = xyz.shape[0]
+        for r in xyz:
+            S_nm[n, 0, 0] += r[0] * r[0]
+            S_nm[n, 1, 1] += r[1] * r[1]
+            S_nm[n, 2, 2] += r[2] * r[2]
+            S_nm[n, 0, 1] += r[0] * r[1]
+            S_nm[n, 0, 2] += r[0] * r[2]
+            S_nm[n, 1, 2] += r[1] * r[2]
+        S_nm[n, 1, 0] = S_nm[n, 0, 1]
+        S_nm[n, 2, 0] = S_nm[n, 0, 2]
+        S_nm[n, 2, 1] = S_nm[n, 1, 2]
+        S_nm[n, :, :] /= N
+    return S_nm

tests/test_order.py

@@ -103,14 +103,7 @@ def test_inertia(traj1, traj2, traj3):
     assert eq(order.compute_inertia_tensor(traj3),
               order._compute_inertia_tensor_slow(traj3))
 
-def test_gyration(traj1, traj2, traj3):
-    assert eq(order.compute_gyration_tensor(traj1),
-              order._compute_gyration_tensor_slow(traj1))
-    assert eq(order.compute_gyration_tensor(traj2),
-              order._compute_gyration_tensor_slow(traj2))
-    assert eq(order.compute_gyration_tensor(traj3),
-              order._compute_gyration_tensor_slow(traj3))
-
+
 def test_nematic_order_args(traj2):
     with pytest.raises(ValueError):
         order.compute_nematic_order(traj2, indices='dog')

tests/test_shape.py

@@ -0,0 +1,101 @@
+##############################################################################
+# MDTraj: A Python Library for Loading, Saving, and Manipulating
+#         Molecular Dynamics Trajectories.
+# Copyright 2012-2017 Stanford University and the Authors
+#
+# Authors: Christoph Klein
+# Contributors: Tim Moore
+#
+# MDTraj is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as
+# published by the Free Software Foundation, either version 2.1
+# of the License, or (at your option) any later version.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with MDTraj. If not, see <http://www.gnu.org/licenses/>.
+##############################################################################
+
+import mdtraj as md
+import numpy as np
+import pytest
+from mdtraj.geometry import shape
+from mdtraj.testing import eq
+
+"""The trajectories `2lines.pdb` and `3lines.pdb` contain several frames of,
+respectively, 2 and 3 "residues" each consisting of two atoms in different
+orientations relative to one another.
+
+`2lines.pdb`
+Frame 0: || in x         - S2 should be 1.0
+Frame 1: || in y         - S2 should be 1.0
+Frame 2: || in z         - S2 should be 1.0
+Frame 3: |- in x/y       - S2 should be 0.25
+Frame 4: |- in y/z       - S2 should be 0.25
+
+`3lines.pdb`
+Frame 0: ||| in x        - S2 should be 1.0
+Frame 1: ||| in y        - S2 should be 1.0
+Frame 2: ||| in z        - S2 should be 1.0
+Frame 3: at right angles - S2 should be 0.0
+"""
+
+@pytest.fixture()
+def traj1(get_fn):
+    return md.load(get_fn('1line.pdb'))
+
+
+@pytest.fixture()
+def traj2(get_fn):
+    return md.load(get_fn('2lines.pdb'))
+
+
+@pytest.fixture()
+def traj3(get_fn):
+    return md.load(get_fn('3lines.pdb'))
+
+@pytest.fixture()
+def traj4(get_fn):
+    return md.load(get_fn('2koc.pdb'))
+
+def test_gyration(traj1, traj2, traj3):
+    assert eq(shape.compute_gyration_tensor(traj1),
+              shape._compute_gyration_tensor_slow(traj1))
+    assert eq(shape.compute_gyration_tensor(traj2),
+              shape._compute_gyration_tensor_slow(traj2))
+    assert eq(shape.compute_gyration_tensor(traj3),
+              shape._compute_gyration_tensor_slow(traj3))
+
+def test_principal_moments(traj4):
+    rg_actual = md.compute_rg(traj4)
+
+    principal_moments = shape.principal_moments(traj4)
+
+    rg_computed = np.sqrt(principal_moments.sum(axis=1))
+
+    assert eq(rg_actual,rg_computed)
+
+def test_asphericity(traj4):
+    b_computed = shape.asphericity(traj4)
+
+    pm = shape.principal_moments(traj4)
+    rg = md.compute_rg(traj4)
+    b_actual = 1.5 * pm[:,2] - rg**2 / 2.0
+
+    assert eq(b_actual,b_computed)
+
+def test_shape_metrics(traj4):
+    b = shape.asphericity(traj4)
+    c = shape.acylindricity(traj4)
+    rg = md.compute_rg(traj4)
+
+    kappa_actual = (b**2 + 0.75*c**2)/(rg**4)
+    kappa_computed = shape.relative_shape_antisotropy(traj4)
+
+    assert eq(kappa_actual,kappa_computed)
+
+