Add coordinate wrapping

[ojeda-e]

In this PR I added PBC conditions. Changes in this PR fix #36 and possibly #22

• self.ag.wrap() added to__init__ and _single_frame.
• Added tests for coordinate wrapping, including negative coordinates in x, y, and z, as suggested in Add test for the generation of z_ref (core_fast_leaflet) on toy system #22
• Since PBC has an effect on the overall results, respective tests for results.z_surface and results.mean were updated.
• I am not sure mentors would be interested in keeping previous tests, so I didn't delete but marked them as xfail(reason="PBC conditions not applied.") (See line 276 is tests_membrane_curvature.py.

@orbeckst @lilyminium @IAlibay are there any other relevant tests I am missing here?
Thanks

[pep8speaks]

Hello @ojeda-e! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

There are currently no PEP 8 issues detected in this Pull Request. Cheers! 🍻

Comment last updated at 2021-08-02 19:09:45 UTC

[codecov]

Codecov Report

Merging #48 (4aa405c) into main (8df8b78) will increase coverage by 2.43%.
The diff coverage is 100.00%.

[lilyminium]

I might be missing something, but could you please explain how the coordinate wrapping treats PBC? As far as I can tell you're still not treating the boundary rows and columns of the surface grid?

I'm also not convinced that you can validly wrap atom coordinates when they have been rotated and translated as I believe you preprocess your trajectories to do. Perhaps you could instead get the translation, or translation + rotation matrix (to center on the protein if any) for each frame, and use that to transform your surface?

This is indirectly related to @orbeckst's project 5 of GSoC (https://github.com/MDAnalysis/mdanalysis/wiki/Project-Ideas-2021#project-5-general-unit-cell-representation), by the way. I'm just pointing that out of interest and as more motivation for the project to go forward rather than suggesting you implement it (although you are welcome to, of course!) Just storing the matrix for each frame of this analysis would be a great starting point.

Edit: sorry if I'm a bit short here, just trying to catch a few bars of signal when I can find it -- which is less than I originally expected

[ojeda-e]

I pushed changes to use wrap instead of PBC. I didn't push changes for the keyword that will discriminate between raw and processed trajectories, as we discussed in this comment.
To clarify, should I add such changes in this PR as well @lilyminium @orbeckst @IAlibay
I think it is good enough with wrap. Instead, I can have clear documentation that says "wrap=True only if your Universe has a raw trajectory." or similar. We don't really need an extra keyword to discriminate the systems.

[ojeda-e]

@orbeckst @lilyminium @IAlibay @richardjgowers do I have a green light here? Thanks.

[richardjgowers]

@ojeda-e looks good but I'd clarify what wrapping is doing here. PBC treatment is confusing as it means one of two very different things

[orbeckst]

I would also add some documentation explaining when you should use wrap=True and explain in words what it achieves.

[ojeda-e]

I would also add some documentation explaining when you should use wrap=True and explain in words what it achieves.

Thanks, @orbeckst, initially I was planning to add the explanation of when to use wrap=True in the Usage (Issue #58)

Now that you point it out I'll add a brief description here in base.py and a more detailed one in Usage.

[ojeda-e]

Can I have an approval here, please? @richardjgowers @orbeckst @lilyminium

[lilyminium]

@ojeda-e just a few suggestions :-)

[ojeda-e]

I deleted the trajectory for now since I would need to create a new trajectory test with u.dimensions[2]=30. Skipping the trajectory the tests pass.

[ojeda-e]

I discover a "little" bug in my test_po4_small.gro (see comment). Given that the size of the box in z was 30 A, ag.wrap was returning z=0 (and therefore everything was zero, surface, mean, Gaussian). I also deleted the trajectory test for now, as I mentioned in this comment. I don't think is strictly necessary, but if you think is best to add the traj, I will add it. I hope we all agree this change doesn't represent any harm or has a negative impact in the tests.

Thanks to @lilyminium for asking me to add tests with wrap=False and wrap=True because otherwise, I'd have never seen that bug. :)

[ojeda-e]

I added two fixtures for universes:

• universe_dummy_wrap with atoms out of bounds in x,
• universe_dummy_wrap_xy with atoms out of bounds in x and y.

Intentionally, I made the two fixtures match with the derived z_surface when wrap=True is applied. Then, the third fixture added dummy_surface. For each one of the universes with atoms out of bounds, I added a test for

• .results.average_surface
• .results.average_mean
• .results.average_gaussian.

For each test that involves curvature, I calculated mean_curvature and gaussian_curvature directly from dummy_surface.

I also discover another bug in surface.py (see comment). But fixed it with an if statement. Hope you are ok with this change as well.

Thanks again @lilyminium

[lilyminium]

Thank you for adding the new tests @ojeda-e!

[ojeda-e]

The pushed changes include:

• Corrected indexes in diagram:

    #   +-----------+          +-----------+
    #   |   | 8 | 9 | 7        | 7 | 8 | 9 |
    #   +-----------+          +-----------+
    # 6 | 4 | 5 |   |   --->   | 4 | 5 | 6 |
    #   +-----------+          +-----------+
    #   |   | 2 | 3 | 1        | 1 | 2 | 3 |
    #   +-----------+          +-----------+
    #

• Added fixtures that returns membranecurvature.run() for each universe with unwrapped coordinates:

  • curvature_unwrapped_universe(self, universe_dummy_wrap)
  • curvature_unwrapped_universe_xy(self, universe_dummy_wrap_xy)

• Removed previous tests for x and y coordinates as suggested, and added an extra one with negative z coordinates with 4 bins.

• In surface.py, prints were replaced by warnings as suggested. However I have a question here:

  • I understand the purpose of logger in base.py. But for the warning in surface.py, I am not sure I should add a logger. Would it make sense to add a logger there?

Thanks @lilyminium

[lilyminium]

You could use the curvature_unwrapped_universe_xy fixture here too, right? (the line after)

[ojeda-e]

Oh, yeah, I missed that, thanks!

[lilyminium]

Thanks @ojeda-e, the new indices are much easier to understand! I found another place you could use a fixture, but otherwise, this looks good to me :D

Would it make sense to add a logger there?

This might come down to personal preference. I don't ever log warnings because users can just logging.captureWarnings(True) in their scripts if they really want it to log. @orbeckst might have a different opinion.

[ojeda-e]

Thanks for checking this @lilyminium! This was a very buggy PR 🐛 but glad it's all fixed now! :D