Conversation
Codecov Report
@@ Coverage Diff @@
## develop #2189 +/- ##
===========================================
+ Coverage 89.48% 89.53% +0.04%
===========================================
Files 157 157
Lines 18819 18883 +64
Branches 2723 2741 +18
===========================================
+ Hits 16841 16906 +65
Misses 1377 1377
+ Partials 601 600 -1
Continue to review full report at Codecov.
|
|
You can try a DNA as an example. In GROMACS the two DNA strands are modeled as individual fragments. An unwrap algorithm should however still keep both strands together. With Gromacs that required several passes over the trajectory. I will check if I can share an example configuration. |
|
We still have the problem that using Another issue is that |
|
I just thought about having a "lazy" |
|
@kain88-de Regarding the DNA example: Right now my implementation does not keep groups of individual compounds together, since this would additionally require making the COMs or COGs of the compounds whole. Such a functionality should either go into a separate method or it should be accessible via an optional kwarg. Furthermore, this would indeed require a method which is able to compute the necessary shifts from the (naturally non-bonded) COMs or COGs of the respective compounds. |
|
@jbarnoud That's a bit hard to judge since #2012 still has some TODOs and as far as I could see several lines of code which need fixing. Nevertheless, I'm pretty sure my implementation already covers most of what #2012 is aiming at. I didn't fully get what the |
richardjgowers
left a comment
richardjgowers
left a comment
There was a problem hiding this comment.
I like this, but a few suggested tweaks
|
|
||
| Note | ||
| ---- | ||
| * This is the inverse transformation of packing atoms into the primary |
There was a problem hiding this comment.
There's some ascii in make_whole which you can steal for this. Maybe also a seealso for pack_into_box (which maybe we can rename wrap?)
There was a problem hiding this comment.
*Group.wrap() already exists and extends the functionality of *Group.pack_into_box(). I'm in favor of deprecating pack_into_box().
There was a problem hiding this comment.
I just realized that pack_into_box() can do out-of-place packing, while wrap() always acts inplace. That should be changed. Furthermore, wrap() potentially acts on atoms outside of the calling group. That's dangerous and should be changed as well.
| if ortho: | ||
| positions = unique_atoms.positions | ||
| # Only unwrap if necessary: | ||
| if np.any(np.abs(positions - positions[0]) > hbox): |
There was a problem hiding this comment.
this optimisation is cute, but maybe deserves to live inside of make_whole
There was a problem hiding this comment.
Very true! That way, the optimization will always be applied.
There was a problem hiding this comment.
also if it's done in the Cython level you can exit the any() loop early, I think this will do the whole calculation then evaluate the boolean values
There was a problem hiding this comment.
For an atomgroup consisting of 500 fragments with 25 atoms per fragment, i measured the following runtimes:
- without optimization: 154 ms,
- optimization in
make_whole(): 65 ms - optimization in Python: 66.5 ms
Turns out the speed-up is not that great anymore. This simple loop runs in 22 ms on the same system:
def group_make_whole(ag):
box = ag.dimensions
fragments = ag.fragments
if np.any(box[3:] != 90.0):
for f in fragments:
mda.lib.mdamath.make_whole(f)
else:
hbox = 0.5 * box[:3]
for f in fragments:
pos = f.positions
if np.any(np.abs(pos - pos[0]) > hbox):
mda.lib.mdamath.make_whole(f)Seems to be the rest of the machinery that slows down group.unwrap().
| try: | ||
| compound_indices = unique_atoms.molnums | ||
| except AttributeError: | ||
| raise NoDataError("Cannot use compound='molecules': " |
There was a problem hiding this comment.
the existing attributeerror should suffice rather than adding another layer here
There was a problem hiding this comment.
I put this here intentionally since the keyword is called 'molecules' but the AttributeError would be raised for molnums, which might be confusing. Same as in *Group.center(), btw.
|
@zemanj |
Where else is it used? I can't find it... |
|
@zemanj this is doing the image shift: https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/core/groups.py#L1251 But more generally I think |
|
@richardjgowers Requests are addressed, I'll write the tests tomorrow. |
|
@zemanj ok cool, this is looking really good. Once we've merged this it should be fairly easy to
|
|
Even though unwrapping requires
One can only accomplish both by either adding a custom |
|
I changed the dtype of |
|
There are at least two minor issues (unnecessary import in |
| .. versionadded:: 0.20.0 | ||
| """ | ||
| atoms = self.atoms | ||
| # bail out early if no bonds in topology: |
There was a problem hiding this comment.
there's a requires decorator that can do this type of thing (and create uniform error messages)
There was a problem hiding this comment.
@richardjgowers I know. The @requires() decorator only works on AtomGroups, that's why it's not used here.
There was a problem hiding this comment.
One could extend the @requires decorator to check for GroupBase instances instead of AtomGroup and then check the group's .atoms member for the given attribute(s). However, I don't like this solution since it
- doesn't allow to check for
ResiduGrouporSegmentGroupattributes and - more importantly, isn't lightweight anymore since it involves object instantiation when accessing the
.atomsmember.
The check I implemented doesn't instantiate additional objects and its impact on performance should therefore be negligible.
There was a problem hiding this comment.
Hmm, really we could have something like assert 'bonds' in u.topology which works... I'll raise an issue. I just like error messages being uniform, so hitting eg 'no masses' always looks similar
|
Ok all checks passed. I'll just add changelog entries and this is good to go. |
|
I restarted the timed out Mac entry in the Travis CI matrix. The test suite ran to completion and then the job seemed to hang. Not sure if that's been happening here lately, but 100 % patch diff coverage is very nice to see--definitely comforting if reviewers don't really have the bandwidth to go over everything. Looks like Richard has checked more thoroughly than I. Does this mean MDAnalysis will soon be able to process my virus trajectories (rhombic dodecahedron with superstructure split all over the place requiring all kinds of trjconv processing & translations)? |
If you mean that your system's periodic unit cell is a rhombic dodecahedron, then the answer is no. |
|
@zemanj I think you can represent a rhombic dodec as a triclinic cell (I think I read somewhere a triclinic cell can represent any 3d tiling shape...) @tylerjereddy I think this will make your unsplitting dreams come true... but if you could check before we merge that would be great. Am I right that you can use a triclinic representation? |
| is_unwrapped = True | ||
| for i in range(1, natoms): | ||
| for j in range(3): | ||
| if fabs(oldpos[i, j] - oldpos[0, j]) >= half_box[j]: |
| from MDAnalysis.lib.mdamath import triclinic_box | ||
|
|
||
|
|
||
| def unwrap_test_universe(have_bonds=True, have_masses=True, have_molnums=True, |
There was a problem hiding this comment.
Yeah writing tests took longer than implementing the feature itself... The cool thing about this is that the minimal test universe makes the test lightning fast and (hopefully) covers all relevant situations. I also designed it with cluster unwrapping in mind.
Not sure about any space-filling 3d shape (sounds plausible, though), but definitely true for a rhombic dodecahedron (gromacs is doing it that way). |
|
I probably won't get a chance to test on my system before the merge, but sounds like these are solid steps in that direction. I guess when I do get around to it, I can open an issue / PR if there's something special that needs tweaking or adding in. |
|
@richardjgowers thanks! 👍 |
5 participants
Fixes (in part) #1004, #1185
Changes made in this Pull Request:
*Group.unwrap(self, box=None, compound='fragments', reference='com', inplace=True)Since I'd really like to have this feature and there wasn't too much progress going on lately in PRs #2012 and #2038, I decided to go ahead and get this going. Using the method in other methods such as
*Group.center*()and the like should be straightforward but is not yet covered here.Since in most cases, boxes are orthorhombic and the majority of molecules are not split accross boundaries, I implemented a short-cut so that only those fragments that are potentially broken are made whole. This yielded a speed-up of a factor between 4 and 5 compared to a simple loop blindly calling
make_whole(fragment)on every fragment in the group (as it's implemented in PR #2038 at the moment).Tests and changelog entry are yet to come.
PR Checklist