Rotation transformation by davidercruz · Pull Request #1937 · MDAnalysis/mdanalysis

davidercruz · 2018-06-13T17:57:43Z

Adds a rotation transformation to the transformations module

Changes made in this Pull Request:

trajectory coordinates may now be rotated by a given angle and an axis. The axis
is formed by a given position or the center of mass/(geometry of an AtomGroup
added tests for both ways of defining a position

PR Checklist

Tests?
Docs?
CHANGELOG updated?
Issue raised/referenced?

coveralls · 2018-06-13T18:25:18Z

Coverage decreased (-0.004%) to 89.963% when pulling 5ed599f on davidercruz:rotation-transform into 22dfcdb on MDAnalysis:develop.

davidercruz · 2018-06-15T17:42:37Z

I just realised that this transformation may cause issues with molecule wrapping/unwrapping. An easy example would be a membrane system rotated by 90 degrees on an axis which is parallel to the plane of the membrane, in a non cubic box. For the system to "fill" the box would require taking atoms from its periodic images, which would change its composition.
Should I add a warning to the docs of the function?
I was looking at the wrapping functions wrap, pack_into_box and apply_PBC, and I don't think they will fail graciously in this case 🤔

jbarnoud · 2018-06-15T18:09:26Z

Rotating a system is always tricky with periodic conditions. @Tsjerk told me once that you can solve the problem by rotating the box definition as well. It may be worth investigating, else a warning in the doc could do. Le 15 juin 2018 7:42 PM, Davide Cruz <notifications@github.com> a écrit :I just realised that this transformation may cause issues with molecule wrapping/unwrapping. An easy example would be a membrane system rotated by 90 degrees on an axis which is parallel to the plane of the membrane, in a non cubic box. For the system to "fill" the box would require taking atoms from its periodic images, which would change it's composition. Should I add a warning to the docs of the function? I was looking at the wrapping functions wrap, pack_into_box and apply_PBC, and I don't think they will fail graciously in this case 🤔 —You are receiving this because your review was requested.Reply to this email directly, view it on GitHub, or mute the thread.

Tsjerk · 2018-06-16T07:23:45Z

The problem is that any rotation of the coordinates causes mismatch with the PBC. As mentioned, the solution is to also rotate the PBC. However, (un)wrapping routines typically rely on certain conditions of the PBC, such as having a first vector aligned with X, the second in the XY plane and the third with positive Z (Gromacs rules). Then (un)wrapping routines may still fail after rotation. The ultimate solution is to store the cumulative rotation of the system (due to fitting etc) alongside PBC and for every (un)wrapping operation, including distance calculations and clustering, rotate the coordinates to match the PBC and rotate back afterwards if so required. Hope it helps, Tsjerk

…

On Fri, Jun 15, 2018, 20:09 Jonathan Barnoud ***@***.***> wrote: Rotating a system is always tricky with periodic conditions. @Tsjerk told me once that you can solve the problem by rotating the box definition as well. It may be worth investigating, else a warning in the doc could do. Le 15 juin 2018 7:42 PM, Davide Cruz ***@***.***> a écrit :I just realised that this transformation may cause issues with molecule wrapping/unwrapping. An easy example would be a membrane system rotated by 90 degrees on an axis which is parallel to the plane of the membrane, in a non cubic box. For the system to "fill" the box would require taking atoms from its periodic images, which would change it's composition. Should I add a warning to the docs of the function? I was looking at the wrapping functions wrap, pack_into_box and apply_PBC, and I don't think they will fail graciously in this case 🤔 —You are receiving this because your review was requested.Reply to this email directly, view it on GitHub, or mute the thread. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1937 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABIX971ezUfryFGQqxTUcW38wlNgKf7Lks5t8_hXgaJpZM4UmqZ0> .

davidercruz · 2018-06-16T13:07:57Z

@Tsjerk in the context of the API, the cumulative rotation could be stored in the TImestep object. Then every one of those operations would read that property of the timestep and perform the rotation.
Another solution is just adding the rotation to the timestep without actually changing the coordinates. Then at the end of the workflow a hidden transformation would rotate the coordinates if a rotation is present in the timestep.

For now I think it's better to warn the user in the docs. What do you think @jbarnoud ?

jbarnoud · 2018-06-16T13:45:05Z

I think the best is to warn the user in the doc. Latter on we can think about a smart unwrap method that knows about the box. Le 16 juin 2018 3:07 PM, Davide Cruz <notifications@github.com> a écrit :@Tsjerk in the context of the API, the cumulative rotation could be stored in the TImestep object. Then every one of those operations would read that property of the timestep and perform the rotation. Another solution is just adding the rotation to the timestep without actually changing the coordinates. Then at the end of the workflow a hidden transformation would rotate the coordinates if a rotation is present in the timestep. For now I think it's better to warn the user in the docs. What do you think @jbarnoud ? —You are receiving this because you were mentioned.Reply to this email directly, view it on GitHub, or mute the thread.

jbarnoud · 2018-06-16T14:36:43Z

package/MDAnalysis/transformations/rotate.py

+from ..lib.transformations import rotation_matrix
+from ..core.groups import AtomGroup
+
+def rotate(angle, direction, center="geometry", **kwargs):


There is no reason to use **kwargs here. It is better to write explicitly all the arguments in the signature. The main issue with **kwargs is that it makes a bunch of argument validation moot: there will not be an error if I make a typo in an argument name for instance, or if I give too many arguments. Also, having a full signature is better on a documentation standpoint.

jbarnoud · 2018-06-16T14:38:12Z

package/MDAnalysis/transformations/rotate.py

+        ts = u.trajectory.ts
+        angle = math.pi/2
+        ag = u.atoms()
+        rotated = MDAnalysis.transformations.rotate(angle, a)(ts)


a is not defined, and ag is not used. Typo I guess.

jbarnoud · 2018-06-16T14:39:54Z

package/MDAnalysis/transformations/rotate.py

+    point: list, optional
+        list of the coordinates of the point from where a custom axis of rotation will
+        be defined. 
+    ts: Timestep


The function does not take ts as an argument.

jbarnoud · 2018-06-16T14:44:29Z

package/MDAnalysis/transformations/rotate.py

+from ..lib.transformations import rotation_matrix
+from ..core.groups import AtomGroup
+
+def rotate(angle, direction, center="geometry", **kwargs):


It may be better to call the function rotateby to be consistent with AtomGroup.rotateby. Also, it may be easier to use AtomGroup.rotateby in your function rather than reimplementing half of it.

richardjgowers · 2018-06-18T13:46:11Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

+from MDAnalysisTests.datafiles import GRO
+
+def test_rotate_custom_position():
+    u = mda.Universe(GRO)


so reading the GRO system requires file i/o etc. better would be to use make_Universe() (in core/test_groups etc). I can't remember if this gives nonzero positions, if it doesn't you can make your own fixture that does.

richardjgowers · 2018-06-18T13:46:25Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

+    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]


np.pi not math.pi

jbarnoud

A few comments on top of Richard's ones. The one on the tests is the most serious: please try to break your code on purpose and see if your tests fail.

jbarnoud · 2018-06-18T13:45:58Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

    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) 


The conventions in PEP8 say that there should not be spaces around the = here, but there should be spaces around the /.

jbarnoud · 2018-06-18T13:49:12Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

    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) 


Here you are overwriting your reference. transformed and ref are the same object so the test cannot fail.

jbarnoud · 2018-06-18T13:49:55Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

 from MDAnalysisTests.datafiles import GRO

-def test_translate_custom_position():
+def test_rotate_custom_position():


Should prabably be test_rotateby_....

…sts now use make_Universe

davidercruz · 2018-06-19T10:46:26Z

As the commit says, the typecasting to numpy.float32 in the rotateby function was causing floating point operation errors, and the tests were failing.
When rotation_matrix is called from within rotateby, the precision of the rotation matrix is reduced in half compared to when it is called outside, e.g.

# matrix from within rotateby 
array([[ 1.000000e+00,  0.000000e+00,  0.000000e+00],
       [ 0.000000e+00, -4.371139e-08, -1.000000e+00],
       [ 0.000000e+00,  1.000000e+00, -4.371139e-08]])

# matrix when calling rotation_matrix outside the transformation
array([[ 1.000000e+00,  0.000000e+00,  0.000000e+00],
       [ 0.000000e+00,  6.123234e-17, -1.000000e+00],
       [ 0.000000e+00,  1.000000e+00,  6.123234e-17]])

And this was causing the tests to fail. I removed all the typecasting in my code, since the rotation_matrix already performs calculations in double precision (np.float64)

richardjgowers · 2018-06-19T15:15:31Z

testsuite/MDAnalysisTests/transformations/test_rotate.py

+from MDAnalysis.lib.transformations import rotation_matrix
+from MDAnalysisTests import make_Universe
+
+def rotateby_universes():


@pytest.fixture() def rotateby_universes():

to make this work. Pytest then inspects arguments to test functions and sees if any names for missing arguments match its known fixtures

richardjgowers · 2018-06-19T15:16:35Z

package/MDAnalysis/transformations/rotate.py

+from ..lib.transformations import rotation_matrix
+from ..core.groups import AtomGroup
+
+def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=[]):


iirc there's some weird stuff that can happen with position=[] as a kwarg. Like I think you can append to that list and the default drifts over time... use =None

richardjgowers · 2018-06-19T15:18:26Z

@davidercruz is there a list somewhere of the transforms you're doing next? ie when is wrap/unwrap happening? :)

davidercruz · 2018-06-19T15:36:53Z

@richardjgowers there is a list in my GSOC project. Today I'll create a PR for a box centering transformation. PBC corrections come next. There is a wrap method in core/groups.py. I'll have to look at it, and see what I can use for the transformation.
Just to confirm, the expected unwrap behaviour is similar to -ur compact in gmx trjconv right?

richardjgowers · 2018-06-19T15:49:41Z

@davidercruz there's also lib.mdamath.make_whole but it's shamefully slow (but I could make it faster if we had a cool transformation that used it...). I'm not great at gromacs, but it'd be cool if it unchopped molecules that had crossed a box boundary and got written to the trajectory in a "packed" way

Tsjerk · 2018-06-19T18:43:11Z

FWIW, what -ur compact of Gromacs intends to do and a possible (un)wrap/boxwrap should do is to put everything in the Voronoi region of the lattice with respect to some reference point, like, e.g., the center of mass of a protein. This may not be the most optimal visualization still, and I think it's a good opportunity to advertise the 'molecular shaped box', which was also introduced in my thesis, and which is the Voronoi region of the lattice with respect to a group of reference points, like, e.g., a protein. To unwrap/repack molecules, a possible strategy is to calculate the circular center of mass and place all atoms of the molecule at the smallest circular distance from that center. The molecule can then also be repositioned to the desired position based on that center. T.

…

On Tue, Jun 19, 2018 at 5:49 PM Richard Gowers ***@***.***> wrote: @davidercruz <https://github.com/davidercruz> there's also lib.mdamath.make_whole but it's shamefully slow (but I could make it faster if we had a cool transformation that used it...). I'm not great at gromacs, but it'd be cool if it unchopped molecules that had crossed a box boundary and got written to the trajectory in a "packed" way — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#1937 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABIX9yRQ_a8MQ9C7VjrQQgh5G_khIw76ks5t-R2ZgaJpZM4UmqZ0> .

-- Tsjerk A. Wassenaar, Ph.D.

jbarnoud · 2018-06-24T14:11:47Z

package/MDAnalysis/transformations/rotate.py

+    pbc_arg = pbc
+    if position and len(position)>2:
+        position = position
+    elif isinstance(ag, AtomGroup):


This has the same issues I pointed in #1946:

it is better not to check for type, but instead try to do the thing (duck typing)

the center is not updated after the first frame

jbarnoud · 2018-06-24T14:12:40Z

package/MDAnalysis/transformations/rotate.py

+from ..core.groups import AtomGroup
+
+def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=None):
+    '''


You mentioned a warning about periodic boundary conditions to add in the documentation. It should be added.

jbarnoud · 2018-06-24T14:14:06Z

package/MDAnalysis/transformations/rotate.py

+    Parameters
+    ----------
+    angle: float
+        rotation angle in radians


All other functions of MDAnalysis that work with angle expect degrees from the user. For consistency, this function should too.

jbarnoud · 2018-06-24T14:14:26Z

package/MDAnalysis/transformations/rotate.py

+    ----------
+    angle: float
+        rotation angle in radians
+    direction: list


jbarnoud · 2018-06-25T07:27:48Z

package/MDAnalysis/transformations/rotate.py

        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 


Thins is indeed the warning. Though, it has nothing to do with the PBC option. It is a more general problem that should be discussed in the body of the docstring. You can add a "Warnings" section.

The pbc option is about wrapping prior to calculating the center. Setting pcs to False is not going to change the fact that the coordinates are modified in a different way than the reference frame.

You're right. I don't know what I was thinking, writing that there 😕

jbarnoud · 2018-06-25T07:28:15Z

package/MDAnalysis/transformations/rotate.py

@@ -76,31 +77,44 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=No
        or 'mass'
    pbc: bool or None, optional


This option has the same issue as the one from the centering transformation.

jbarnoud · 2018-06-25T07:29:47Z

package/MDAnalysis/transformations/rotate.py

-    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


To be consistent with the other transformation, the option should be named point. The name of the user facing variable is more important than the name of the variable we manipulate inside; so the variable in the wrapper should be renamed so the argument is consistent to the user.

jbarnoud · 2018-06-25T07:36:21Z

package/MDAnalysis/transformations/rotate.py

-            raise ValueError('{} is not a valid argument for center'.format(center))
+    angle = np.deg2rad(angle)
+    if position:
+        if len(position)!=3:


I mentioned it already; len on an array is a terrible idea. An array with a (1, 3) shape will have length 3, the same goes for an array of shape (3, 3). First convert to an array with np.asarray, then compare the shape with position.shape != (3, ). You may be permissive and accept (1, 3) as a shape as well.

jbarnoud · 2018-06-25T07:39:15Z

I made a mistake and commented on a previous state. Still have a look at my last batch of comments, they are still all valid.

jbarnoud · 2018-06-25T07:43:40Z

A couple more things:

pay attention to keep you documentation in sync with the code: your examples are still in radian;
your code and your tests rely on the same code: if the rotation_matrix function is broken your tests will not break; make sure that there are tests for that function;
why are you slicing the rotattion matrix?

davidercruz · 2018-06-25T10:05:28Z

@jbarnoud I'm slicing the rotation matrix because the output matrix is 4 by 4 by default.
I don't know why it does so. This is what happens in the code:

def rotation_matrix(angle, direction, point=None):
    # things are done before
    M = np.identity(4) # M is 4 by 4
    M[:3, :3] = R
    # things are done after
    return M # M is still 4 by 4

orbeckst · 2018-06-26T00:03:31Z

@davidercruz note on "Why is M 4x4?": They are used for affine transformations, see #1939 (comment) for some notes. (That PR has nothing to do with this one but I happened to think about it a few days ago.)

jbarnoud · 2018-07-06T07:10:37Z

package/MDAnalysis/transformations/rotate.py

+from ..lib.transformations import rotation_matrix
+from ..core.groups import AtomGroup
+
+def rotateby(angle, direction, point=None, center="geometry", pbc=False, ag=None):


Shouldn't the pbc argument be wrap to be consistent with the remove center transformation?

jbarnoud · 2018-07-06T21:24:47Z

The tests are failing. I strongly encourage you to run the tests locally when you do changes: https://github.com/MDAnalysis/mdanalysis/wiki/UnitTests#developers.

Added rotation transformation, tests, docs and update CHANGELOG

fe6dbb6

jbarnoud self-requested a review June 14, 2018 11:08

jbarnoud requested changes Jun 16, 2018

View reviewed changes

changed name to rotateby; fixed typos ; dropped kwargs

81808d8

richardjgowers reviewed Jun 18, 2018

View reviewed changes

jbarnoud requested changes Jun 18, 2018

View reviewed changes

davidercruz added 2 commits June 19, 2018 11:31

typecasting to np.float32 causes FLOP errors with rotation_matrix; te…

fc03cad

…sts now use make_Universe

some code cleanup

fbe09ea

more cleanup

cfe7b78

richardjgowers reviewed Jun 19, 2018

View reviewed changes

some fixes; position now defaults to None

76e6b03

jbarnoud reviewed Jun 24, 2018

View reviewed changes

davidercruz and others added 2 commits June 24, 2018 23:30

updated docs; center is now updated every frame; added more tests

2d9f989

Merge branch 'develop' into rotation-transform

9e67c9d

jbarnoud reviewed Jun 25, 2018

View reviewed changes

davidercruz added 2 commits June 25, 2018 11:18

no longer checking for lengths ; added test to rotation matrix function

67142c0

typo in doc

cf9c6c5

Merge branch 'develop' into rotation-transform

39ef719

jbarnoud reviewed Jul 6, 2018

View reviewed changes

changed pbc to wrap

f3e74c9

updated tests

5ed599f

jbarnoud approved these changes Jul 7, 2018

View reviewed changes

jbarnoud merged commit 1a1b2f1 into MDAnalysis:develop Jul 7, 2018

jbarnoud added this to the 0.19.0 milestone Jul 7, 2018

jbarnoud self-assigned this Jul 7, 2018

davidercruz mentioned this pull request Jul 18, 2018

Added fitting transformations #1991

Merged

4 tasks

		@@ -76,31 +77,44 @@ def rotateby(angle, direction, center="geometry", pbc=None, ag=None, position=No
		or 'mass'
		pbc: bool or None, optional

Conversation

davidercruz commented Jun 13, 2018 • edited by orbeckst Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

PR Checklist

Uh oh!

coveralls commented Jun 13, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

davidercruz commented Jun 15, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jbarnoud commented Jun 15, 2018 via email

Uh oh!

Tsjerk commented Jun 16, 2018 via email

Uh oh!

davidercruz commented Jun 16, 2018

Uh oh!

jbarnoud commented Jun 16, 2018 via email

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jbarnoud left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

davidercruz commented Jun 19, 2018

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

richardjgowers commented Jun 19, 2018

Uh oh!

davidercruz commented Jun 19, 2018

Uh oh!

richardjgowers commented Jun 19, 2018

Uh oh!

Tsjerk commented Jun 19, 2018 via email

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jbarnoud commented Jun 25, 2018

Uh oh!

jbarnoud commented Jun 25, 2018

Uh oh!

davidercruz commented Jun 25, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

orbeckst commented Jun 26, 2018

davidercruz commented Jun 13, 2018 •

edited by orbeckst

Loading

coveralls commented Jun 13, 2018 •

edited

Loading

davidercruz commented Jun 15, 2018 •

edited

Loading

davidercruz commented Jun 25, 2018 •

edited

Loading