Box centering transformation by davidercruz · Pull Request #1946 · MDAnalysis/mdanalysis

davidercruz · 2018-06-19T18:20:31Z

Adds a centering function

This function takes an AtomGroup as argument and will translate it's center of mass/geometry to the center of the unit cell. Optional arguments include a user defined box, AtomGroup center type (mass or geometry) and whether to use PBC or not in AtomGroup center calculation.

PR Checklist

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

coveralls · 2018-06-19T18:44:33Z

Coverage increased (+0.02%) to 89.91% when pulling 94f0196 on davidercruz:center-transform into 691e785 on MDAnalysis:develop.

jbarnoud

Just some comments from a very quick read through. I'll review the PR in more details latter.

Your tests should cover erroneous uses to make sure it fails as expected. Please add tests for each exception you raise. There are many reasons to do that but here are the two that comes to my mind right away:

exceptions are part of the API, if the way the exceptions are raises changes we break some user's code and we want to know about it as soon as possible;
raising an exception is code, but it is hopefully rarely executed, this means that we are less likely to catch en error in that code. It became painfully obvious that raising an exception is code when we ported MDAnalysis to python 3 and replaced the old %-based string formating by the more recent format-based one; suddenly we had a lot of unusable errors while the tests were happily passing.

Also you need to test more cases:

is the pbc argument properly honored?
is the correct box used?
what if I give none of the optional argument?
what if I give only some?

The tests are not only there to make sure that your code works now. They are mostly there to make sure that if I do something stupid because I try to translate half of the code base to perl at 3 AM after a mixture of beer and coffee, everything I break becomes red!

jbarnoud · 2018-06-20T07:12:03Z

package/MDAnalysis/transformations/translate.py

+    """
+
+    pbc_arg = pbc
+    if isinstance(ag, AtomGroup):


Rather than explicitly test for the type of the argument, it would be better to try running the method and raise the correct error if it fails. This allows somebody to give as argument anything that roughly behave like an atom group. The python philosophy in that matter is that if it looks like a duck and quacks like a duck, then it is a duck.

try: ag_center = ag.center_of_geometry(...) except AttributeError: raise ValueError('{:r} is not an AtomGroup'.format(ag))

jbarnoud · 2018-06-20T07:14:38Z

package/MDAnalysis/transformations/translate.py

+        else:
+            boxcenter = np.sum(ts.triclinic_dimensions, axis=0) / 2
+        vector = boxcenter - ag_center
+        translate(vector)(ts)


While it is good to reuse existing function, I am worried this would create a big overhead for just a sum.

jbarnoud · 2018-06-20T07:16:56Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+    return reference, transformed
+
+def test_translate():
+    ref_u, trans_u = translate_universes()


The expected use of a fixture is to let pytest call the function. Here, it should be:

def test_translate(translate_universes): ref_u, trans_u = translate_universes ...

davidercruz · 2018-06-20T19:03:27Z

@jbarnoud Thank you.
I added some more tests. One of the tests made me realize that there was no distinction between failing because the input arg was not an AtomGroup object, or because it did not have masses. Now it's a bit clearer - but I think it can be better.

One thing regarding the box argument when pbc=True: currently the custom box is not taken into account when calculating the center of mass/geometry with PBC. I could change these two methods to correct this. What do you think?

jbarnoud · 2018-06-21T06:33:02Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+    with pytest.raises(ValueError):
+        translate(vector)(ts)
+
+def test_center_in_box_args(translate_universes):


This should be one test function per error. The way it is written, the function will stop at the first failing assertion. However, knowing if the other ones pass or fail is a valuable information to figure out what happens.

jbarnoud · 2018-06-21T06:34:37Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+    with pytest.raises(ValueError): 
+        center_in_box(ag, pbc=True)(ts)
+    # what if a wrong center type name is passed?
+    bad_center = "charge"


We may implement one day a center of charge. Having something obviously wrong such as "invalid" makes the test future proof and clearer (though, not so much clearer as the variable name is already good).

jbarnoud · 2018-06-21T06:37:43Z

testsuite/MDAnalysisTests/transformations/test_translate.py

-

-def test_center():
+def test_translate_vector(translate_universes):


This makes me realize: how is the dtype of ts.positions affected by the dtype of the vector?

I had an error in #1902 where if the vector was not np.float32 some floating point errors would pop up, and the transformed positions would be slightly different from the expected

jbarnoud · 2018-06-21T06:38:32Z

testsuite/MDAnalysisTests/transformations/test_translate.py

    reference = make_Universe(trajectory=True)
-    transformed = make_Universe(trajectory=True)
-    transformed.trajectory.ts.dimensions = np.array([372.,373.,374.,90,90,90])
+    transformed = make_Universe(['masses'], trajectory=True)


You may want to assign masses so that the center of mass and the center of geometry are different.

using the ['masses'} extra in make_Universe generates masses for the atoms. The centers are different for the transformed universe

jbarnoud · 2018-06-21T06:43:16Z

One thing regarding the box argument when pbc=True: currently the custom box is not taken into account when calculating the center of mass/geometry with PBC. I could change these two methods to correct this. What do you think?

What would be the use case for a custom box in your transformation? In what case would you want to use a box size different from the one of the time step?

I do not see a use case, which does not mean there isn't one. If there is indeed no use case, then the box option of the transformation should probably be removed.

davidercruz · 2018-06-21T12:45:08Z

@jbarnoud one use case would be to center one the atomgroup of a trajectory in the unit cell of another from a different universe, although the user could simply change the dimensions of the first to match those of the latter.

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

Then, for sure, you do not want to use the custom box for the center of mass or center of geometry as it will mess up your coordinates. It may also be less confusing to give a point to center on rather than a full custom box.

jbarnoud · 2018-06-23T08:46:59Z

package/MDAnalysis/transformations/translate.py

+    point: list, optional
+        overrides the unit cell center - the coordinates of the Timestep are translated so
+        that the center of mass/geometry of the given AtomGroup is aligned to this position
+        instead. Defined as a list of size 3. Example: [1, 2, 3]


"array like" rather than "list"

jbarnoud · 2018-06-23T08:47:57Z

package/MDAnalysis/transformations/translate.py

-        if box:
-            boxcenter = np.sum(triclinic_vectors(box), axis=0) / 2
+        if point:
+            if len(point)==3:


This could be done once in the parent function rather than every frame here.

jbarnoud

Some more changes. But it is getting there.

jbarnoud · 2018-06-23T08:51:08Z

package/CHANGELOG

  * 0.18.1

 Enhancements
+  * Added a box centering trajectory transformation


Refer to the PR. At the end of the day, it may be worth grouping all the transformations in a single changelog entry.

jbarnoud · 2018-06-23T08:52:09Z

package/MDAnalysis/transformations/translate.py

@@ -47,9 +50,7 @@ def translate(vector):
    ----------
    vector: list


This is not supposed to be a list, but an array of anything similar. I think "array-like" is the term.

jbarnoud · 2018-06-23T08:52:31Z

package/MDAnalysis/transformations/translate.py

+    center: str, optional
+        used to choose the method of centering on the given atom group. Can be 'geometry'
+        or 'mass'
+    point: list, optional


Same as above, "array-like".

jbarnoud · 2018-06-23T08:56:29Z

package/MDAnalysis/transformations/translate.py

+    pbc_arg = pbc
+    try:
+        if center == 'geometry':
+            ag_center = ag.center_of_geometry(pbc=pbc_arg)


This cannot be done in the parent function. The way it is written now means that the center is not updated every frame. What you can do, though, is to choose the method in the parent method:

# In the parent function if center =='geometry': center_method = functools.partial(ag.center_of_geometry, pbc=pbc_arg) elif ... # In the wrapped function ag_center = center_method()

jbarnoud · 2018-06-23T15:45:38Z

package/MDAnalysis/transformations/translate.py

+    def wrapped(ts):
+        ag_center = center_method()
+        if point:
+            if len(point)==3:


The comment disappeared, but this remains to change: you are testing and converting the point to an array for each frame while it could be done once in the parent function.

…nalysis into center-transform

jbarnoud · 2018-06-23T19:01:59Z

package/MDAnalysis/transformations/translate.py


    pbc_arg = pbc
+    point = np.array(point)
+    if point.any() and len(point)!=3:


if point is [0, 0, 0], then point.any() will be False. Also, I would recommend avoiding len on arrays: I got weird behaviours at times. Anyway, the logic you had before where point could be an array or None was misplaced but safer.

jbarnoud · 2018-06-23T19:02:28Z

package/MDAnalysis/transformations/translate.py

    """

    pbc_arg = pbc
+    point = np.array(point)


Shouldn't you specify the dtype? You had it specified, earlier.

jbarnoud · 2018-06-23T19:03:05Z

package/MDAnalysis/transformations/translate.py

-                boxcenter = np.float32(point)
-            else:
-                raise ValueError('{} is not a valid point'.format(point))
+        if point.any():


Same as before: point.any() is not safe here.

jbarnoud · 2018-06-23T19:04:42Z

Also, look at the doc failure on travis.

jbarnoud

One last thing :)

jbarnoud · 2018-06-24T10:36:50Z

package/MDAnalysis/transformations/translate.py

+
+    return wrapped
+
+def center_in_box(ag, center='geometry', point=None, pbc=None):


What does it mean to have pbc=None? Is there a difference from having the default being True or False? If there is such a difference, it should be explicitly stated in the docstring; if not, then the default should be the correct boolean. At the moment, I do not now if PBC is taken into account by default or not.

I tried it and it really makes no difference. However, the default in both center_of_mass and center_of_geometry methods is None, so I tried to keep things consistent. I added the default to the docs, however, to make things clearer for the user.

Then you should figure out what the default means for center_of_mass and center_of_geometry. I could already see that None was the default. The question remains: what am I going to get if I do not specify pbc?

Oh you were talking about the case when flags['use_pbc']

edit: I saw @richardjgowers 's comment. Nevermind

davidercruz · 2018-06-24T20:37:47Z

@richardjgowers Oh! Sorry.

Ok then, I will change None to False, since there's really only two intended behaviours

jbarnoud

I think this may be my last round of comments.

jbarnoud · 2018-07-02T14:56:32Z

package/MDAnalysis/transformations/translate.py

-
-.. autofunction:: translate   
-    
+The vector can either be user defined, using the function `translate`


It should be

:func:`translate`

in order to create a link to the function. Similar treatment should be applied bellow for center_in_box.

jbarnoud · 2018-07-02T14:58:45Z

package/MDAnalysis/transformations/translate.py

+    Translates the coordinates of a given :class:`~MDAnalysis.coordinates.base.Timestep`
+    instance so that the center of geometry/mass of the given :class:`~MDAnalysis.core.groups.AtomGroup`
+    is centered on the unit cell. The unit cell dimensions are taken from the input Timestep object or
+    can be defined using the `box` argument.


The last bit is out of date.

jbarnoud · 2018-07-02T14:59:32Z

package/MDAnalysis/transformations/translate.py

+    point: array-like, optional
+        overrides the unit cell center - the coordinates of the Timestep are translated so
+        that the center of mass/geometry of the given AtomGroup is aligned to this position
+        instead. Defined as an array of size 3. Example: [1, 2, 3]


The example is confusing, just remove it.

jbarnoud · 2018-07-02T15:02:30Z

package/MDAnalysis/transformations/translate.py

+        if point is None:
+            boxcenter = np.sum(ts.triclinic_dimensions, axis=0) / 2
+        else:
+            boxcenter = np.float32(point)


At this point, point already has the correct dtype, no need to convert it again.

jbarnoud · 2018-07-02T15:03:18Z

package/MDAnalysis/transformations/translate.py

+        return ts
+
+    return wrapped
+


Remove the trailing lines.

jbarnoud · 2018-07-02T15:17:27Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+def test_center_in_box_bad_ag(translate_universes):
+    # this universe as a box size zero
+    ts = translate_universes[0].trajectory.ts
+    ag = translate_universes[0].residues[0].atoms


You do not need ag in that test.

jbarnoud · 2018-07-02T15:21:48Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+    with pytest.raises(ValueError): 
+        center_in_box(bad_ag)(ts)
+
+def test_center_in_box_bad_point(translate_universes):


You should parametrize that test as described above for the bad vector.

jbarnoud · 2018-07-02T15:23:40Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+        center_in_box(ag, center=bad_center)(ts)
+
+def test_center_in_box_no_masses(translate_universes):   
+    # this universe as a box size zero


This is not the comment you care about. What matters here is that the universe does not have masses.

jbarnoud · 2018-07-02T15:24:32Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+@pytest.fixture()
+def translate_universes():
+    # create the Universe objects for the tests
+    reference = make_Universe(trajectory=True)


Some tests depend on that universe not having dimensions and masses. There should be a comment here warning about these requirements.

jbarnoud · 2018-07-02T15:25:41Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+
+def test_center_in_box_coords_no_options(translate_universes):
+    # what happens when we center the coordinates arround the center of geometry of a residue?
+    ref_u = translate_universes[0]


This is fine, but you can also write:

ref_u, trans_u = translate_universes

richardjgowers · 2018-07-02T21:28:59Z

package/MDAnalysis/transformations/translate.py

+from functools import partial
+
+from ..lib.mdamath import triclinic_vectors
+from ..core.groups import AtomGroup


I'm not sure you need this AtomGroup import here, maybe it's old.

In general these transformations should use duck typing, ie not explicitly checking the type of things coming in, but simply expecting them to do all the things an AtomGroup does.

Indeed, the import is not used anymore.

richardjgowers · 2018-07-02T21:32:01Z

package/MDAnalysis/transformations/translate.py

+        overrides the unit cell center - the coordinates of the Timestep are translated so
+        that the center of mass/geometry of the given AtomGroup is aligned to this position
+        instead. Defined as an array of size 3.
+    pbc: bool, optional


So you've been consistent, but imho pbc is ambiguous. Could we call this wrap instead?

richardjgowers · 2018-07-02T21:39:29Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+#
+


These tests are very good, but it would be good to test that an AtomGroup's positions are correctly modified (ie more of an integration test rather than solely relying on unit tests). Something like:

u = mda.Universe() u.trajectory.add_transformation() ref = .... assert u.atoms.positions = ref

It seems to be tested now.

…ext ; DummyReader didn't have a _transformations variable

jbarnoud

Tiny details to fix.

jbarnoud · 2018-07-04T17:26:19Z

package/MDAnalysis/transformations/translate.py

+from functools import partial
+
+from ..lib.mdamath import triclinic_vectors
+from ..core.groups import AtomGroup


Indeed, the import is not used anymore.

jbarnoud · 2018-07-04T17:27:55Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+#
+


It seems to be tested now.

jbarnoud · 2018-07-04T17:28:55Z

testsuite/MDAnalysisTests/transformations/test_translate.py

+    transformed = make_Universe(['masses'], trajectory=True)
+    transformed.trajectory.ts.dimensions = np.array([372., 373., 374., 90, 90, 90])
+    return reference, transformed
+


PEP8 says there should be 2 new lines between functions. Please fix in the whole PR.

davidercruz · 2018-07-04T18:00:13Z

@jbarnoud Regarding the center in plane transformations: what kind of planes do I use? xy xz and yz? Or also other custom planes?

jbarnoud · 2018-07-04T18:13:20Z

@davidercruz I'll be happy with the 3 main planes. In practice, the xy plane is going to be the most used. You can implement custom planes, but be careful that it should not impair the performances of the base feature.

davidercruz added 3 commits June 19, 2018 19:03

added a box centering transformation; added tests ; updated changelog

28c8072

cleanup

60d7195

improved docs

1ccb380

davidercruz changed the title ~~Center transform~~ Box centering transformation Jun 19, 2018

import division from the future

d9134fa

jbarnoud requested changes Jun 20, 2018

View reviewed changes

clarified exceptions ; moar tests

e95b5db

jbarnoud reviewed Jun 21, 2018

View reviewed changes

split test on args

cff6c00

replaced box with point; updated tests; fixed some typos

bf2faa9

jbarnoud reviewed Jun 23, 2018

View reviewed changes

jbarnoud requested changes Jun 23, 2018

View reviewed changes

davidercruz and others added 2 commits June 23, 2018 14:46

atom group center is now calculated inside wrapped ; updated changelog

51709be

Merge branch 'develop' into center-transform

9123478

jbarnoud reviewed Jun 23, 2018

View reviewed changes

davidercruz added 4 commits June 23, 2018 18:08

point check and conversion now done in parent function

bf62c1e

cleaned point checking

1775462

Merge branch 'center-transform' of https://github.com/davidercruz/mda…

55ef308

…nalysis into center-transform

cleaned point checking

1cac08e

jbarnoud reviewed Jun 23, 2018

View reviewed changes

code cleanup

cbe9049

jbarnoud requested changes Jun 24, 2018

View reviewed changes

defult pbc flag now in docs

d32abc3

jbarnoud mentioned this pull request Jun 24, 2018

Rotation transformation #1937

Merged

4 tasks

pbc description now includes reference to flags

c7b6622

davidercruz added 2 commits June 24, 2018 21:46

removed flags reference;pbc arg defaults to False

aff6615

no longer checking for lengths

0f6613f

jbarnoud requested changes Jul 2, 2018

View reviewed changes

davidercruz added 2 commits July 2, 2018 18:45

code cleanup; added parameters to tests

aecda76

forgot :func:

7cc6cbc

richardjgowers reviewed Jul 2, 2018

View reviewed changes

pbc argument is now wrap; added tests in the transformations API cont…

acbbf43

…ext ; DummyReader didn't have a _transformations variable

jbarnoud requested changes Jul 4, 2018

View reviewed changes

style cleanup

94f0196

jbarnoud approved these changes Jul 4, 2018

View reviewed changes

jbarnoud merged commit b718350 into MDAnalysis:develop Jul 4, 2018

jbarnoud self-assigned this Jul 4, 2018

jbarnoud added this to the 0.18.x milestone Jul 4, 2018

davidercruz mentioned this pull request Jul 22, 2018

corrected rotateby transformation calculation #2000

Merged

4 tasks



		def test_center():
		def test_translate_vector(translate_universes):

		@@ -47,9 +50,7 @@ def translate(vector):
		----------
		vector: list


		return wrapped

		def center_in_box(ag, center='geometry', point=None, pbc=None):


		.. autofunction:: translate

		The vector can either be user defined, using the function `translate`

Conversation

davidercruz commented Jun 19, 2018

PR Checklist

Uh oh!

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

Uh oh!

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 20, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

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 Jun 21, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

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 21, 2018

Uh oh!

davidercruz commented Jun 21, 2018

Uh oh!

jbarnoud commented Jun 22, 2018

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!

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 23, 2018

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 Jun 24, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

davidercruz commented Jun 24, 2018

Uh oh!

jbarnoud left a comment

coveralls commented Jun 19, 2018 •

edited

Loading

davidercruz commented Jun 20, 2018 •

edited

Loading

davidercruz Jun 21, 2018 •

edited

Loading

davidercruz Jun 24, 2018 •

edited

Loading