From fe6dbb6cd9fe73011896d114a042351e438a7a66 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Wed, 13 Jun 2018 18:02:36 +0100 Subject: [PATCH 01/11] Added rotation transformation, tests, docs and update CHANGELOG --- package/CHANGELOG | 1 + .../MDAnalysis/transformations/__init__.py | 4 + package/MDAnalysis/transformations/rotate.py | 115 ++++++++++++++++++ .../transformations/test_rotate.py | 56 +++++++++ 4 files changed, 176 insertions(+) create mode 100644 package/MDAnalysis/transformations/rotate.py create mode 100644 testsuite/MDAnalysisTests/transformations/test_rotate.py diff --git a/package/CHANGELOG b/package/CHANGELOG index 7aaaeb1cadc..7b05f427682 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,6 +19,7 @@ The rules for this file: * 0.18.1 Enhancements + * Added a rotation coordinate transformation * Added a on-the-fly trajectory transformations API and a coordinate translation function (Issue #786) * Added various duecredit stubs diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index eaf3d0a74eb..c95103ec73b 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -53,6 +53,8 @@ def wrapped(ts): Currently implemented transformations are: - translate: translate the coordinates of a given trajectory frame by a given vector. + - rotate: rotates the coordinates by a given angle arround an axis formed by a direction + and a point @@ -83,5 +85,7 @@ def wrapped(ts): from __future__ import absolute_import from .translate import translate +from .rotate import rotate + diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py new file mode 100644 index 00000000000..9121e69f8bf --- /dev/null +++ b/package/MDAnalysis/transformations/rotate.py @@ -0,0 +1,115 @@ +# -*- 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 ..lib.transformations import rotation_matrix +from ..core.groups import AtomGroup + +def rotate(angle, direction, center="geometry", **kwargs): + ''' + Rotates the trajectory by a given angle on a given axis. The axis is defined by + the user, combining the direction vector and a position. This position 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 pi/2 on a x axis centered on a given atom group: + + .. code-block:: python + + ts = u.trajectory.ts + angle = math.pi/2 + ag = u.atoms() + rotated = MDAnalysis.transformations.rotate(angle, a)(ts) + + e.g. rotate the coordinates by a custom axis: + + .. code-block:: python + + ts = u.trajectory.ts + angle = math.pi/2 + p = [1,2,3] + d = [0,0,1] + rotated = MDAnalysis.transformations.rotate(angle, point=point, direction=d)(ts) + + Parameters + ---------- + angle: float + rotation angle in radians + direction: list + 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' + pbc: bool or None, optional + If True, all the atoms from the given AtomGroup will be moved to the unit cell + before calculating the center + point: list, optional + list of the coordinates of the point from where a custom axis of rotation will + be defined. + ts: Timestep + frame that will be transformed + + Returns + ------- + :class:`~MDAnalysis.coordinates.base.Timestep` object + + ''' + theta = np.float32(angle) + vector = np.float32(direction) + position = kwargs.pop("position", []) + ag = kwargs.pop("ag", None) + pbc_arg = kwargs.pop("pbc", None) + if len(position)>2: + position = np.float32(position) + elif isinstance(ag, AtomGroup): + if center == "geometry": + position = ag.center_of_geometry(pbc = pbc_arg) + elif center == "mass": + position = ag.center_of_mass(pbc = pbc_arg) + else: + raise ValueError('{} is not a valid argument for center'.format(center)) + else: + raise ValueError('A position or an AtomGroup must be specified') + + def wrapper(ts): + rotation = rotation_matrix(theta, vector, position)[:3, :3] + ts.positions= np.dot(ts.positions, rotation) + + return ts + + return wrapper + \ 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..a8dc3efc8f3 --- /dev/null +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -0,0 +1,56 @@ +# -*- 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 +import math +from numpy.testing import assert_array_almost_equal + +import MDAnalysis as mda +from MDAnalysis.transformations import rotate +from MDAnalysis.lib.transformations import rotation_matrix +from MDAnalysisTests.datafiles import GRO + +def test_translate_custom_position(): + u = mda.Universe(GRO) + vector = np.float32([1,0,0]) + pos = np.float32([0,0,0]) + matrix = rotation_matrix(math.pi/2, vector, pos)[:3, :3] + ref = u.trajectory.ts + ref.positions = np.dot(ref.positions, matrix) + transformed = rotate(math.pi/2, vector, position = pos)(ref) + assert_array_almost_equal(transformed, ref.positions, decimal=6) + +def test_translate_atomgroup(): + u = mda.Universe(GRO) + vector = np.float32([1,0,0]) + ref_pos = np.float32([52.953686, 44.179996, 29.397898]) + matrix = rotation_matrix(math.pi/2, vector, ref_pos)[:3, :3] + ref = u.trajectory.ts + ref.positions = np.dot(ref.positions, matrix) + selection = u.atoms.select_atoms('resid 1') + transformed = rotate(math.pi/2, vector, ag = selection, center='geometry')(ref) + assert_array_almost_equal(transformed, ref.positions, decimal=6) + + \ No newline at end of file From 81808d8776faf41e906ee62c72367f71997111e9 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Mon, 18 Jun 2018 12:05:08 +0100 Subject: [PATCH 02/11] changed name to rotateby; fixed typos ; dropped kwargs --- package/MDAnalysis/transformations/__init__.py | 2 +- package/MDAnalysis/transformations/rotate.py | 12 ++++-------- .../MDAnalysisTests/transformations/test_rotate.py | 10 +++++----- 3 files changed, 10 insertions(+), 14 deletions(-) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index c95103ec73b..3ca410b0e24 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -85,7 +85,7 @@ def wrapped(ts): from __future__ import absolute_import from .translate import translate -from .rotate import rotate +from .rotate import rotateby diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 9121e69f8bf..2a79ff6a51a 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -36,7 +36,7 @@ from ..lib.transformations import rotation_matrix from ..core.groups import AtomGroup -def rotate(angle, direction, center="geometry", **kwargs): +def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[]): ''' Rotates the trajectory by a given angle on a given axis. The axis is defined by the user, combining the direction vector and a position. This position can be the center @@ -49,7 +49,7 @@ def rotate(angle, direction, center="geometry", **kwargs): ts = u.trajectory.ts angle = math.pi/2 ag = u.atoms() - rotated = MDAnalysis.transformations.rotate(angle, a)(ts) + rotated = MDAnalysis.transformations.rotate(angle, ag)(ts) e.g. rotate the coordinates by a custom axis: @@ -77,11 +77,9 @@ def rotate(angle, direction, center="geometry", **kwargs): pbc: bool or None, optional If True, all the atoms from the given AtomGroup will be moved to the unit cell before calculating the center - point: list, optional + position: list, optional list of the coordinates of the point from where a custom axis of rotation will be defined. - ts: Timestep - frame that will be transformed Returns ------- @@ -90,9 +88,7 @@ def rotate(angle, direction, center="geometry", **kwargs): ''' theta = np.float32(angle) vector = np.float32(direction) - position = kwargs.pop("position", []) - ag = kwargs.pop("ag", None) - pbc_arg = kwargs.pop("pbc", None) + pbc_arg = pbc if len(position)>2: position = np.float32(position) elif isinstance(ag, AtomGroup): diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index a8dc3efc8f3..b25c41b4a8f 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -28,21 +28,21 @@ from numpy.testing import assert_array_almost_equal import MDAnalysis as mda -from MDAnalysis.transformations import rotate +from MDAnalysis.transformations import rotateby from MDAnalysis.lib.transformations import rotation_matrix from MDAnalysisTests.datafiles import GRO -def test_translate_custom_position(): +def test_rotate_custom_position(): u = mda.Universe(GRO) vector = np.float32([1,0,0]) pos = np.float32([0,0,0]) matrix = rotation_matrix(math.pi/2, vector, pos)[:3, :3] ref = u.trajectory.ts ref.positions = np.dot(ref.positions, matrix) - transformed = rotate(math.pi/2, vector, position = pos)(ref) + transformed = rotateby(math.pi/2, vector, position = pos)(ref) assert_array_almost_equal(transformed, ref.positions, decimal=6) -def test_translate_atomgroup(): +def test_rotate_atomgroup(): u = mda.Universe(GRO) vector = np.float32([1,0,0]) ref_pos = np.float32([52.953686, 44.179996, 29.397898]) @@ -50,7 +50,7 @@ def test_translate_atomgroup(): ref = u.trajectory.ts ref.positions = np.dot(ref.positions, matrix) selection = u.atoms.select_atoms('resid 1') - transformed = rotate(math.pi/2, vector, ag = selection, center='geometry')(ref) + transformed = rotateby(math.pi/2, vector, ag = selection, center='geometry')(ref) assert_array_almost_equal(transformed, ref.positions, decimal=6) \ No newline at end of file From fc03cad6ac466355fe346a1f5e2176a08e9d5798 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 19 Jun 2018 11:31:25 +0100 Subject: [PATCH 03/11] typecasting to np.float32 causes FLOP errors with rotation_matrix; tests now use make_Universe --- package/MDAnalysis/transformations/rotate.py | 10 ++-- .../transformations/test_rotate.py | 47 +++++++++++-------- 2 files changed, 33 insertions(+), 24 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 2a79ff6a51a..6c023fd3f1b 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -86,16 +86,16 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[] :class:`~MDAnalysis.coordinates.base.Timestep` object ''' - theta = np.float32(angle) - vector = np.float32(direction) + theta = angle + vector = direction pbc_arg = pbc if len(position)>2: - position = np.float32(position) + position = position elif isinstance(ag, AtomGroup): if center == "geometry": - position = ag.center_of_geometry(pbc = pbc_arg) + position = ag.center_of_geometry(pbc=pbc_arg) elif center == "mass": - position = ag.center_of_mass(pbc = pbc_arg) + position = ag.center_of_mass(pbc=pbc_arg) else: raise ValueError('{} is not a valid argument for center'.format(center)) else: diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index b25c41b4a8f..349ceeb4bab 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -24,33 +24,42 @@ import numpy as np import pytest -import math 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.datafiles import GRO +from MDAnalysisTests import make_Universe -def test_rotate_custom_position(): - u = mda.Universe(GRO) - vector = np.float32([1,0,0]) - pos = np.float32([0,0,0]) - matrix = rotation_matrix(math.pi/2, vector, pos)[:3, :3] - ref = u.trajectory.ts +def rotateby_universes(): + # create the Universe objects for the tests + reference = make_Universe(trajectory=True) + transformed = make_Universe(trajectory=True) + return reference, transformed + +def test_rotateby_custom_position(): + # what happens when we use a custom position for the axis of rotation? + ref_u, trans_u = rotateby_universes() + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1,0,0] + pos = [0,0,0] + matrix = rotation_matrix(np.pi / 2, vector, pos)[:3, :3] ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(math.pi/2, vector, position = pos)(ref) - assert_array_almost_equal(transformed, ref.positions, decimal=6) + transformed = rotateby(np.pi / 2, vector, position=pos)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) -def test_rotate_atomgroup(): - u = mda.Universe(GRO) - vector = np.float32([1,0,0]) - ref_pos = np.float32([52.953686, 44.179996, 29.397898]) - matrix = rotation_matrix(math.pi/2, vector, ref_pos)[:3, :3] - ref = u.trajectory.ts +def test_rotateby_atomgroup(): + # what happens when we rotate arround the center of geometry of a residue? + ref_u, trans_u = rotateby_universes() + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + center_pos = [6,7,8] + vector = [1,0,0] + matrix = rotation_matrix(np.pi, vector, center_pos)[:3, :3] ref.positions = np.dot(ref.positions, matrix) - selection = u.atoms.select_atoms('resid 1') - transformed = rotateby(math.pi/2, vector, ag = selection, center='geometry')(ref) - assert_array_almost_equal(transformed, ref.positions, decimal=6) + selection = trans_u.residues[0].atoms + transformed = rotateby(np.pi, vector, ag=selection, center='geometry')(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) \ No newline at end of file From fbe09eab442f879a0643661ee2735af65e755f3d Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 19 Jun 2018 11:34:25 +0100 Subject: [PATCH 04/11] some code cleanup --- package/MDAnalysis/transformations/rotate.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 6c023fd3f1b..b7b9828070b 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -86,8 +86,6 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[] :class:`~MDAnalysis.coordinates.base.Timestep` object ''' - theta = angle - vector = direction pbc_arg = pbc if len(position)>2: position = position @@ -102,7 +100,7 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[] raise ValueError('A position or an AtomGroup must be specified') def wrapper(ts): - rotation = rotation_matrix(theta, vector, position)[:3, :3] + rotation = rotation_matrix(angle, direction, position)[:3, :3] ts.positions= np.dot(ts.positions, rotation) return ts From cfe7b78e40bf4fb6c5c608b923852546966ed0e7 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 19 Jun 2018 15:13:08 +0100 Subject: [PATCH 05/11] more cleanup --- package/MDAnalysis/transformations/rotate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index b7b9828070b..edb6a0b277f 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -99,11 +99,11 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[] else: raise ValueError('A position or an AtomGroup must be specified') - def wrapper(ts): + def wrapped(ts): rotation = rotation_matrix(angle, direction, position)[:3, :3] ts.positions= np.dot(ts.positions, rotation) return ts - return wrapper + return wrapped \ No newline at end of file From 76e6b0312ca27edc870b4f93a263f2508391fe5c Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Tue, 19 Jun 2018 16:58:49 +0100 Subject: [PATCH 06/11] some fixes; position now defaults to None --- package/MDAnalysis/transformations/rotate.py | 4 ++-- testsuite/MDAnalysisTests/transformations/test_rotate.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index edb6a0b277f..05a0efc29ac 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -36,7 +36,7 @@ from ..lib.transformations import rotation_matrix from ..core.groups import AtomGroup -def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[]): +def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=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 position. This position can be the center @@ -87,7 +87,7 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[] ''' pbc_arg = pbc - if len(position)>2: + if position and len(position)>2: position = position elif isinstance(ag, AtomGroup): if center == "geometry": diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 349ceeb4bab..7eb79e60807 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -31,6 +31,7 @@ from MDAnalysis.lib.transformations import rotation_matrix from MDAnalysisTests import make_Universe +@pytest.fixture() def rotateby_universes(): # create the Universe objects for the tests reference = make_Universe(trajectory=True) From 2d9f989a63c75cc9904677855f3925850277ada3 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Sun, 24 Jun 2018 23:30:18 +0100 Subject: [PATCH 07/11] updated docs; center is now updated every frame; added more tests --- package/MDAnalysis/transformations/rotate.py | 48 +++--- .../transformations/test_rotate.py | 140 ++++++++++++++++-- 2 files changed, 159 insertions(+), 29 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 05a0efc29ac..0a89dd6c33e 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -32,11 +32,12 @@ import math import numpy as np +from functools import partial from ..lib.transformations import rotation_matrix from ..core.groups import AtomGroup -def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=None): +def rotateby(angle, direction, position=None, center="geometry", pbc=None, 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 position. This position can be the center @@ -64,8 +65,8 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=No Parameters ---------- angle: float - rotation angle in radians - direction: list + 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 @@ -76,31 +77,44 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=No or 'mass' pbc: bool or None, optional If True, all the atoms from the given AtomGroup will be moved to the unit cell - before calculating the center - position: list, optional + before calculating the center. Warning: Wrapping/unwrapping the trajectory or + performing PBC corrections may not be possible after rotating the trajectory. + position: 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 - + ''' pbc_arg = pbc - if position and len(position)>2: - position = position - elif isinstance(ag, AtomGroup): - if center == "geometry": - position = ag.center_of_geometry(pbc=pbc_arg) - elif center == "mass": - position = ag.center_of_mass(pbc=pbc_arg) - else: - raise ValueError('{} is not a valid argument for center'.format(center)) + angle = np.deg2rad(angle) + if position: + if len(position)!=3: + raise ValueError('{} is not a valid point'.format(position)) + 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 position or an AtomGroup must be specified') + raise ValueError('A position or an AtomGroup must be specified') def wrapped(ts): - rotation = rotation_matrix(angle, direction, position)[:3, :3] + if position: + point = position + else: + point = center_method() + rotation = rotation_matrix(angle, direction, point)[:3, :3] ts.positions= np.dot(ts.positions, rotation) return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 7eb79e60807..e3b50e50d18 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -32,35 +32,151 @@ from MDAnalysisTests import make_Universe @pytest.fixture() -def rotateby_universes(): +def rotate_universes(): # create the Universe objects for the tests reference = make_Universe(trajectory=True) - transformed = 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_rotateby_custom_position(): +def test_rotateby_custom_position(rotate_universes): # what happens when we use a custom position for the axis of rotation? - ref_u, trans_u = rotateby_universes() + 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] - matrix = rotation_matrix(np.pi / 2, vector, pos)[:3, :3] + angle = 90 + matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(np.pi / 2, vector, position=pos)(trans) + transformed = rotateby(angle, vector, position=pos)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) -def test_rotateby_atomgroup(): - # what happens when we rotate arround the center of geometry of a residue? - ref_u, trans_u = rotateby_universes() +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] - matrix = rotation_matrix(np.pi, vector, center_pos)[:3, :3] + 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(np.pi, vector, ag=selection, center='geometry')(trans) + 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', pbc=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', pbc=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, position=bad_position)(ts) - \ No newline at end of file +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, pbc=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 position 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 From 67142c051dfff5dd9c0198637035e5ef7fdad2f8 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Mon, 25 Jun 2018 11:18:13 +0100 Subject: [PATCH 08/11] no longer checking for lengths ; added test to rotation matrix function --- package/MDAnalysis/transformations/rotate.py | 40 +++++++++++-------- .../transformations/test_rotate.py | 30 ++++++++++++-- 2 files changed, 49 insertions(+), 21 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 0a89dd6c33e..d17a5111a97 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -37,10 +37,10 @@ from ..lib.transformations import rotation_matrix from ..core.groups import AtomGroup -def rotateby(angle, direction, position=None, center="geometry", pbc=None, ag=None): +def rotateby(angle, direction, point=None, center="geometry", pbc=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 position. This position can be the center + 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 pi/2 on a x axis centered on a given atom group: @@ -48,7 +48,7 @@ def rotateby(angle, direction, position=None, center="geometry", pbc=None, ag=No .. code-block:: python ts = u.trajectory.ts - angle = math.pi/2 + angle = 90 ag = u.atoms() rotated = MDAnalysis.transformations.rotate(angle, ag)(ts) @@ -57,7 +57,7 @@ def rotateby(angle, direction, position=None, center="geometry", pbc=None, ag=No .. code-block:: python ts = u.trajectory.ts - angle = math.pi/2 + angle = 90 p = [1,2,3] d = [0,0,1] rotated = MDAnalysis.transformations.rotate(angle, point=point, direction=d)(ts) @@ -75,24 +75,30 @@ def rotateby(angle, direction, position=None, center="geometry", pbc=None, ag=No center: str, optional used to choose the method of centering on the given atom group. Can be 'geometry' or 'mass' - pbc: bool or None, optional - If True, all the atoms from the given AtomGroup will be moved to the unit cell - before calculating the center. Warning: Wrapping/unwrapping the trajectory or - performing PBC corrections may not be possible after rotating the trajectory. - position: array-like, optional + pbc: 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 = pbc angle = np.deg2rad(angle) - if position: - if len(position)!=3: - raise ValueError('{} is not a valid point'.format(position)) + 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': @@ -107,14 +113,14 @@ def rotateby(angle, direction, position=None, center="geometry", pbc=None, ag=No else: raise ValueError('{} is not an AtomGroup object'.format(ag)) else: - raise ValueError('A position or an AtomGroup must be specified') + raise ValueError('A point or an AtomGroup must be specified') def wrapped(ts): - if position: - point = position + if point is None: + position = center_method() else: - point = center_method() - rotation = rotation_matrix(angle, direction, point)[:3, :3] + position = point + rotation = rotation_matrix(angle, direction, position)[:3, :3] ts.positions= np.dot(ts.positions, rotation) return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index e3b50e50d18..20d82a765fc 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -39,8 +39,30 @@ def rotate_universes(): 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 position for the axis of rotation? + # 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 @@ -50,7 +72,7 @@ def test_rotateby_custom_position(rotate_universes): angle = 90 matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(angle, vector, position=pos)(trans) + 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): @@ -136,7 +158,7 @@ def test_rotateby_bad_position(rotate_universes): vector = [0, 0, 1] bad_position = [1] with pytest.raises(ValueError): - rotateby(angle, vector, position=bad_position)(ts) + rotateby(angle, vector, point=bad_position)(ts) def test_rotateby_bad_pbc(rotate_universes): # this universe as a box size zero @@ -176,7 +198,7 @@ def test_rotateby_no_args(rotate_universes): ts = rotate_universes[0].trajectory.ts angle = 90 vector = [0, 0, 1] - # if no position or AtomGroup are passed to the function + # 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 From cf9c6c5f068e741a235fefc4b171cfe892232374 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Mon, 25 Jun 2018 11:20:00 +0100 Subject: [PATCH 09/11] typo in doc --- package/MDAnalysis/transformations/rotate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index d17a5111a97..a3644e74e14 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -43,7 +43,7 @@ def rotateby(angle, direction, point=None, center="geometry", pbc=False, ag=None 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 pi/2 on a x axis centered on a given atom group: + e.g. rotate the coordinates by 90 degrees on a x axis centered on a given atom group: .. code-block:: python From f3e74c948e5012222de4f1fa4a6cd58000bec05a Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 6 Jul 2018 15:35:56 +0100 Subject: [PATCH 10/11] changed pbc to wrap --- package/MDAnalysis/transformations/rotate.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index a3644e74e14..a9712d2a44a 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -35,9 +35,8 @@ from functools import partial from ..lib.transformations import rotation_matrix -from ..core.groups import AtomGroup -def rotateby(angle, direction, point=None, center="geometry", pbc=False, ag=None): +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 @@ -75,7 +74,7 @@ def rotateby(angle, direction, point=None, center="geometry", pbc=False, ag=None center: str, optional used to choose the method of centering on the given atom group. Can be 'geometry' or 'mass' - pbc: bool, optional + 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. @@ -93,7 +92,7 @@ def rotateby(angle, direction, point=None, center="geometry", pbc=False, ag=None after rotating the trajectory. ''' - pbc_arg = pbc + pbc_arg = wrap angle = np.deg2rad(angle) if point: point = np.asarray(point, np.float32) From 5ed599f1ffb47dff76257a1e9831a68e512d7286 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 6 Jul 2018 22:33:30 +0100 Subject: [PATCH 11/11] updated tests --- testsuite/MDAnalysisTests/transformations/test_rotate.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 20d82a765fc..f22e7771770 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -120,7 +120,7 @@ def test_rotateby_atomgroup_cog_pbc(rotate_universes): 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', pbc=True)(trans) + 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): @@ -136,7 +136,7 @@ def test_rotateby_atomgroup_com_pbc(rotate_universes): 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', pbc=True)(trans) + 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): @@ -169,7 +169,7 @@ def test_rotateby_bad_pbc(rotate_universes): angle = 90 vector = [0, 0, 1] with pytest.raises(ValueError): - rotateby(angle, vector, ag = ag, pbc=True)(ts) + rotateby(angle, vector, ag = ag, wrap=True)(ts) def test_rotateby_bad_center(rotate_universes): # this universe as a box size zero