diff --git a/package/CHANGELOG b/package/CHANGELOG index 5f99db203b0..d17be062b56 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,6 +19,8 @@ The rules for this file: * 0.18.1 Enhancements + + * Added a rotation coordinate transformation (PR #1937) * Added a box centering trajectory transformation (PR #1946) * Added a on-the-fly trajectory transformations API and a coordinate translation function (Issue #786) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index 94444791fbf..8c7b388c20a 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -56,9 +56,9 @@ def wrapped(ts): - translate: translate the coordinates of a given trajectory frame by a given vector. - center_in_box: translate the coordinates of a given trajectory frame so that a given - AtomGroup is centered in the unit cell - - + AtomGroup is centered in the unit cell + - rotateby: rotates the coordinates by a given angle arround an axis formed by a direction + and a point Examples -------- @@ -93,6 +93,8 @@ def wrapped(ts): from __future__ import absolute_import from .translate import translate, center_in_box +from .rotate import rotateby + diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py new file mode 100644 index 00000000000..a9712d2a44a --- /dev/null +++ b/package/MDAnalysis/transformations/rotate.py @@ -0,0 +1,128 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +"""\ +Rotate trajectory --- :mod:`MDAnalysis.transformations.translate` +================================================================= + +Rotates the coordinates by a given angle arround an axis formed by a direction and a +point + +""" +from __future__ import absolute_import + +import math +import numpy as np +from functools import partial + +from ..lib.transformations import rotation_matrix + +def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=None): + ''' + Rotates the trajectory by a given angle on a given axis. The axis is defined by + the user, combining the direction vector and a point. This point can be the center + of geometry or the center of mass of a user defined AtomGroup, or a list defining custom + coordinates. + e.g. rotate the coordinates by 90 degrees on a x axis centered on a given atom group: + + .. code-block:: python + + ts = u.trajectory.ts + angle = 90 + ag = u.atoms() + rotated = MDAnalysis.transformations.rotate(angle, ag)(ts) + + e.g. rotate the coordinates by a custom axis: + + .. code-block:: python + + ts = u.trajectory.ts + angle = 90 + p = [1,2,3] + d = [0,0,1] + rotated = MDAnalysis.transformations.rotate(angle, point=point, direction=d)(ts) + + Parameters + ---------- + angle: float + rotation angle in degrees + direction: array-like + vector that will define the direction of a custom axis of rotation from the + provided point. + ag: AtomGroup, optional + use this to define the center of mass or geometry as the point from where the + rotation axis will be defined + center: str, optional + used to choose the method of centering on the given atom group. Can be 'geometry' + or 'mass' + wrap: bool, optional + If `True`, all the atoms from the given AtomGroup will be moved to the unit cell + before calculating the center of mass or geometry. Default is `False`, no changes + to the atom coordinates are done before calculating the center of the AtomGroup. + point: array-like, optional + list of the coordinates of the point from where a custom axis of rotation will + be defined. + + Returns + ------- + :class:`~MDAnalysis.coordinates.base.Timestep` object + + Warning + ------- + Wrapping/unwrapping the trajectory or performing PBC corrections may not be possible + after rotating the trajectory. + + ''' + pbc_arg = wrap + angle = np.deg2rad(angle) + if point: + point = np.asarray(point, np.float32) + if point.shape != (3, ) and point.shape != (1, 3): + raise ValueError('{} is not a valid point'.format(point)) + elif ag: + try: + if center == 'geometry': + center_method = partial(ag.center_of_geometry, pbc=pbc_arg) + elif center == 'mass': + center_method = partial(ag.center_of_mass, pbc=pbc_arg) + else: + raise ValueError('{} is not a valid argument for center'.format(center)) + except AttributeError: + if center == 'mass': + raise AttributeError('{} is not an AtomGroup object with masses'.format(ag)) + else: + raise ValueError('{} is not an AtomGroup object'.format(ag)) + else: + raise ValueError('A point or an AtomGroup must be specified') + + def wrapped(ts): + if point is None: + position = center_method() + else: + position = point + rotation = rotation_matrix(angle, direction, position)[:3, :3] + ts.positions= np.dot(ts.positions, rotation) + + return ts + + return wrapped + \ No newline at end of file diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py new file mode 100644 index 00000000000..f22e7771770 --- /dev/null +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -0,0 +1,204 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +from __future__ import absolute_import + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +import MDAnalysis as mda +from MDAnalysis.transformations import rotateby +from MDAnalysis.lib.transformations import rotation_matrix +from MDAnalysisTests import make_Universe + +@pytest.fixture() +def rotate_universes(): + # create the Universe objects for the tests + reference = make_Universe(trajectory=True) + transformed = make_Universe(['masses'], trajectory=True) + transformed.trajectory.ts.dimensions = np.array([372., 373., 374., 90, 90, 90]) + return reference, transformed + +def test_rotation_matrix(): + # test if the rotation_matrix function is working properly + # simple angle and unit vector on origin + angle = 180 + vector = [0, 0, 1] + pos = [0, 0, 0] + ref_matrix = np.asarray([[-1, 0, 0], + [0, -1, 0], + [0, 0, 1]], np.float64) + matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] + assert_array_almost_equal(matrix, ref_matrix, decimal=6) + # another angle in a custom axis + angle = 60 + vector = [1, 2, 3] + pos = [1, 2, 3] + ref_matrix = np.asarray([[ 0.53571429, -0.6229365 , 0.57005291], + [ 0.76579365, 0.64285714, -0.01716931], + [-0.35576719, 0.44574074, 0.82142857]], np.float64) + matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] + assert_array_almost_equal(matrix, ref_matrix, decimal=6) + + +def test_rotateby_custom_position(rotate_universes): + # what happens when we use a custom point for the axis of rotation? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1,0,0] + pos = [0,0,0] + angle = 90 + matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] + ref.positions = np.dot(ref.positions, matrix) + transformed = rotateby(angle, vector, point=pos)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + +def test_rotateby_atomgroup_cog_nopbc(rotate_universes): + # what happens when we rotate arround the center of geometry of a residue + # without pbc? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + center_pos = [6,7,8] + vector = [1,0,0] + angle = 90 + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] + ref.positions = np.dot(ref.positions, matrix) + selection = trans_u.residues[0].atoms + transformed = rotateby(angle, vector, ag=selection, center='geometry')(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + +def test_rotateby_atomgroup_com_nopbc(rotate_universes): + # what happens when we rotate arround the center of mass of a residue + # without pbc? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1,0,0] + angle = 90 + selection = trans_u.residues[0].atoms + center_pos = selection.center_of_mass() + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] + ref.positions = np.dot(ref.positions, matrix) + transformed = rotateby(angle, vector, ag=selection, center='mass')(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + +def test_rotateby_atomgroup_cog_pbc(rotate_universes): + # what happens when we rotate arround the center of geometry of a residue + # with pbc? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1,0,0] + angle = 90 + selection = trans_u.residues[0].atoms + center_pos = selection.center_of_geometry(pbc=True) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] + ref.positions = np.dot(ref.positions, matrix) + transformed = rotateby(angle, vector, ag=selection, center='geometry', wrap=True)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + +def test_rotateby_atomgroup_com_pbc(rotate_universes): + # what happens when we rotate arround the center of mass of a residue + # with pbc? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1,0,0] + angle = 90 + selection = trans_u.residues[0].atoms + center_pos = selection.center_of_mass(pbc=True) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] + ref.positions = np.dot(ref.positions, matrix) + transformed = rotateby(angle, vector, ag=selection, center='mass', wrap=True)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + +def test_rotateby_bad_ag(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + ag = rotate_universes[0].residues[0].atoms + # what happens if something other than an AtomGroup is given? + angle = 90 + vector = [0, 0, 1] + bad_ag = 1 + with pytest.raises(ValueError): + rotateby(angle, vector, ag = bad_ag)(ts) + +def test_rotateby_bad_position(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + # what if the box is in the wrong format? + angle = 90 + vector = [0, 0, 1] + bad_position = [1] + with pytest.raises(ValueError): + rotateby(angle, vector, point=bad_position)(ts) + +def test_rotateby_bad_pbc(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + ag = rotate_universes[0].residues[0].atoms + # is pbc passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size + angle = 90 + vector = [0, 0, 1] + with pytest.raises(ValueError): + rotateby(angle, vector, ag = ag, wrap=True)(ts) + +def test_rotateby_bad_center(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + ag = rotate_universes[0].residues[0].atoms + # what if a wrong center type name is passed? + angle = 90 + vector = [0, 0, 1] + bad_center = " " + with pytest.raises(ValueError): + rotateby(angle, vector, ag = ag, center=bad_center)(ts) + +def test_rotateby_no_masses(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + ag = rotate_universes[0].residues[0].atoms + # if the universe has no masses and `mass` is passed as the center arg + angle = 90 + vector = [0, 0, 1] + bad_center = "mass" + with pytest.raises(AttributeError): + rotateby(angle, vector, ag = ag, center=bad_center)(ts) + +def test_rotateby_no_args(rotate_universes): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + angle = 90 + vector = [0, 0, 1] + # if no point or AtomGroup are passed to the function + # it should raise a ValueError + with pytest.raises(ValueError): + rotateby(angle, vector)(ts) \ No newline at end of file