-
Notifications
You must be signed in to change notification settings - Fork 41
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Blog post for on-the-fly transformations #87
Conversation
def up_by_2(): | ||
|
||
def wrapped(ts): | ||
# here's where the magic happens |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
There is no need for a wrapped function here. The transformation could be written:
def up_by_2(ts):
ts.positions += np.asarray([0,0,20])
return ts
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
True, but maybe it's easiest to always define transformations in this form so it's not weird when you need to do things with arguments. So maybe this blogpost shouldn't have this simple case that can flatten out @davidercruz
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In an ideal world, this would be the perfect example. However, I would just add a footnote saying that in this particular case we could also write as https://github.com/MDAnalysis/MDAnalysis.github.io/pull/87/files#r209500911 says but the pattern described here is general and will also work for situations when the simple approach doesn't work.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for the blog post.
I like the comparison with trjconv and the large number of examples. I think you should introduce the usage API explicitly after the first example and only have one comparison to trjconv, in the first example. If you always show the trjconv example it is not clear why MDAnalysis is superior because trjconv can do the same thing. So more focus on the MDAnalysis transformations and show of examples that are not possible with other tools currently. If you want to sell it we have to show what new things people can do with this API.
|
||
# Using the on-the-fly trajectory transformations in MDAnalysis | ||
|
||
On-the-fly transformations have been introduced in version 0.18.1 of MDAnalysis. This long awaited feature brings to MDAnalysis a whole new level of functionality, allowing for more efficient workflows when analyzing and visualizing simulation trajectories. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's not really long awaited. It's OK to write that this is your GSoC project.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is work you did and you can say that.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we link to docs from MDAnalysis as well? How much of this is documented in the code yet?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Each PR that I have adds a topic to the transformations in the online docs. So only the merged PR's are actually in the docs
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What is already in the docs (you can link to the dev docs at https://www.mdanalysis.org/mdanalysis) and which PRs are still pending?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Rotations and some translations are already in the docs. #1973 and #2038 are yet to be merged. Both include docs.
|
||
## Using MDAnalysis transformations | ||
Now it's time to learn how to use the trajectory transformations in MDAnalysis. During the following steps, we will apply some transformations on a 1 ns trajectory of a simple 19-residue peptide embeded in a 128-DMPC membrane, showing the GROMACS `trjconv` command and the equivalent MDAnalysis code and output. To keep thinks lightweight, frames are were taken every 100 ps, and water molecules were removed. | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
So here it would be cool to show that MDAnalysis can remove water, remove pbc and align a trajectory all in one go. Or at least mention it. It might not be super obvious to new users
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you mean using a rot+fit and make_whole transformation on a trajectory with water and then writing out the trajectory without?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
thinks -> things
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
embeded -> embedded
|
||
### Example 1: making everything whole again | ||
When performing MD simulations using periodic boundary conditions, molecules will often cross the limits of the unit cell. When this happens, some atoms of the molecule will show up on the the opposing side of the unit cell and some molecular viewers will show stretched bonds and other visual artifacts depending on the visual representation of the system. | ||
This is the case of our system. Without any modifications, when we look at the trajectory of our system, things will become very ugly, very fast: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
remove the ugly. The simulation is totally fine! The bonds we use for drawing are just broken producing admittingly ugly stripes in the visualization but the simulation is fine!
transformation = mda.transformations.unwrap(ag) | ||
# we add the transformation to the trajectory in it will be applied everytime we read a frame | ||
u.trajectory.add_transformations(transformation) | ||
# we can do other things with the trajectory without having to generate a file with the transformed |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This should be it's own paragraph separate from this and not hidden in a comment.
|
||
And we choose `Protein` as the group to be centered. | ||
|
||
In MDAnalysis we use the `center_in_box` transformation. This function takes an AtomGroup as a mandatory argument. Optional arguments include `weights`, which is used to calculate the weighted center of the given AtomGroup (if weights='mass' then the center of mass is calculated), `center_to` which is used when the user needs to center the AtomGroup in a custom point instead of the center of the unit cell, and `wrap` which, if `True`, causes all the atoms of the AtomGroup to be moved to the unit cell before calculating the weighted center. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
where in the code is wrap
set to True
? That sentence doesn't make sense to me. Also please explain what the compound keyword does here. You should also explain what the atomgroup argument to center_in_box
does. Will it only center those atoms? Is this the reference around which all atoms are moved?
NGLWidget(count=11) | ||
|
||
|
||
You can see that the `transformations` workflow above has three steps: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
explain this before. Say the transformation now consists of three steps and explain them like you do here. Then show the code. The code doesn't come as a surprise then.
transformations = (mda.transformations.unwrap(ag), | ||
mda.transformations.center_in_box(prot), | ||
mda.transformations.wrap(ag, compound='fragments'), | ||
mda.transformations.fit_rot_trans(prot, reference, plane='xy', weights="mass")) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Here you can enter another wrap call. This means, in the end, all atoms are nicely within the unit-cell without any lipids wandering off. That is especially a transformation that is not possible with gromacs tools.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Wrapping after a rotation is not advised. The rotation breaks the relationship between the coordinate and the box so the wrapping will be erroneous.
There is a comment about this in at least one of the docstring. Note also that the issue is not specific to MDAnalysis but that gromacs has it as well.
There is a solution for that but it is not trivial. The box matrix could be rotated, but it is not a common thing to do and it will likely break assumptions about how triclinic boxes are handled both in MDAnalysis and other softwares.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Unrelated to the PR, but related to this topic: is it impossible to center the system on a protein, i.e. helix lying on the membrane, and wrap the lipids such that the membrane does not rotate in circles? This would be neat for visualization / density calculations of membrane parameters in relation to the protein.
As far as I'm aware it's impossible to do with GROMACS, but I was hoping that the related PR would enable the functionality in MDAnalysis.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Do you mean fit instead of center? Because if you just center the helix in the unit cell then there are no rotation operations performed on the coordinates. Thus you can still perform pbc transformations - wrap and unwrap. This can also be done in gromacs.
If you do rotational fitting though, there are no transformations currently in MDAnalysis that can do that. As you rotate the system, it becomes "disconnected" from the unit cell , making pbc transformations trickier. It can be done, but the resulting system would be different - in membrane systems it would be mostly fine if the membrane is pure. But for mixtures it would throw the lipid disposition off.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I was talking about the latter, rotational fitting. If the membrane would be wrapped after rotation, one could more easily do analysis of lipid accumulation, thickness/density calculations. Just thought it might be easy to implement / was already implemented in the recent PRs.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Possibly add a footnote to discuss this case.
I agree with @mimischi that I'd love to have the unwrap after rotation functionality. Maybe another GSoC project...
|
||
## The advantage of using MDAnalysis for trajectory transformations | ||
Each simulation package is often bundled with tools to transform and analyze trajectories, such as GROMACS' `trjconv`. However, most of the times, the user is required to apply all the intended transformations to the whole trajectory (or the portion of interest) prior to visualization and analysis. This often requires processing huge files, sometimes more than once. Moreover, some tools such as `trjconv` do not support frame indexing for the most popular trajectory formats, requiring iterating over frames that are not needed for that particular analysis. | ||
Trajectory transformations in MDAnalysis, on the other end, have one great advantage - they are performed on-the-fly. This means that, after loading the trajectory file, the user adds a transformation workflow and the transformations are applied to each frame that is read. There is no need to iterate over the whole trajectory before performing other analysis, and, with the frame indexing provided by MDAnalysis, only the frames of interest are processed. Moreover, the way the transformations API is implemented makes it really easy to add custom transformations. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This without a noun following it is bad writing style.
Trajectory transformations in MDAnalysis, on the other end, have one great advantage - they are performed on-the-fly for each frame that is read. Transformations are added to a universe as a transformation workflow containing one or more transformations. The API also makes it easy to add new transformations for your own projects.
A good system to test against would be a DNA. It's difficult to do with gromacs since the two strands are not a connected topology, this often leads gromacs to break the DNA into two pieces. I'll look if I can share a few ns from a DNA simulation I have. |
@davidercruz do you need more input to make progress here? It would be good to get this online before 0.19.0 is released so that people can immediately see one of the cool new features. I think this write-up will become a initial version for part of the MDA 1.0 paper. Think about it this way: it should highlight the most important parts of your work. Show examples that you find worthwhile. Mention examples that change the way people might want to do their analysis. Anything that enables something new is highly interesting. |
@jbarnoud can I please ask you to shepherd the write-up to completion? Thanks! |
@orbeckst sorry, I've been busy with work after I took my vacations. I have some corrections that I made with @kain88-de 's recomendations. I will commit them soon. |
Is the wrap transformation actually in 0.19.0 or is this PR MDAnalysis/mdanalysis#2038 which is still pending? In this case you need to first show examples with what's available and then show things that are coming soon. Looking at the docs https://www.mdanalysis.org/docs/documentation_pages/trajectory_transformations.html I see only translations (which is actually empty – I'll try to add this to PR MDAnalysis/mdanalysis#2105) and rotations. |
The PR MDAnalysis/mdanalysis#1973 contains the changes to the docs for the translations module. Regarding wrap/unwrap, I think they are the most important transformations and announcing the on-the-fly transformations feature without them would rob the feature some of its impact. I can help close the PR as soon as possible, if changes are needed. |
I agree. We should really try to get your open PRs in ASAP. However, the fact is that the features aren't in 0.19.0. My suggestion: I think we should therefore only briefly mention the transformations in the 0.19.0 blog post PR #95 along the lines that the foundations are there, it's experimental, but cool things are coming. You can then publish this blog post after the 0.19.0 release announcement and we try to get the PRs done for a 0.20.0 (or 1.0.0) release. (I'll post the same comment in the discussion for #94 and we continue the discussion there.) |
You add them under public/images and include them with <img
src="{{site.images}}/numfocus.png"
title="NumFOCUS Foundation" alt="NumFOCUS Foundation"
style="float: right; width: 10em;" /> or similar code. |
import MDAnalysis as mda | ||
import nglview as nv | ||
u = mda.Universe('pept_in_memb.tpr', 'pept_in_memb.xtc') | ||
nv.show_mdanalysis(u) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@lilyminium would there be a way to embed the nglview widget here with actual content?
|
||
```python | ||
nv.show_mdanalysis(u) | ||
``` |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@lilyminium would there be a way to embed the nglview widget here with actual content — it would be neat if users could see the changes.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you really wanted to, I think you could do this (@hainm would know more):
- Run the code in a Jupyter notebook and write view out to an HTML file
view = nv.show_mdanalysis(u)
view # you must visualise it in the browser
view.write_html('transformed.html', [view], frame_range=[0, len(u.trajectory), 1])
- Scrap everything outside the
<body>
tags intransformed.html
(optional if in iframe?) {% include transformed.html %}
in your post or put it in an iframe
Otherwise making it a gif and including the gif is probably easier.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I have created the html files. Do I simply paste the code inside the <body>
tags in the md file?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
You could, but I think that using the include tag would be functionally the same and neater.
I think you could also make an iframe.html
that looks like this, stick it in _includes/
<iframe src={{ include.nv_html }}> </iframe>
and in your .md file go:
{% include iframe.html nv_html="transformed.html" %}
which hopefully would remove the need to delete the stuff outside <body />
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If that doesn't render well, here is an example of using the movie maker to make a gif (scroll up a bit).
I am feeling like a gravedigger here, trying to exhume this PR... but I think we have all transformations now in 0.21/1.0 https://www.mdanalysis.org/mdanalysis/documentation_pages/trajectory_transformations.html so we should also have this post. @davidercruz could you try to run everything with the latest develop and check that it still works? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a cool feature!
|
||
## Using MDAnalysis transformations | ||
Now it's time to learn how to use the trajectory transformations in MDAnalysis. During the following steps, we will apply some transformations on a 1 ns trajectory of a simple 19-residue peptide embeded in a 128-DMPC membrane, showing the GROMACS `trjconv` command and the equivalent MDAnalysis code and output. To keep thinks lightweight, frames are were taken every 100 ps, and water molecules were removed. | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
thinks -> things
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hooray! I'll update the date and merge.
- changed 0.18.1 to 0.19.0 - link to gsoc post - small edits
…nt center_in_axis and center_in_plane
e0dd09e
to
4012363
Compare
Rebased & changed date. I might spend a few minutes to look if I can quickly embed some nglview visualizations/gifs. That would help a lot. |
- signed post - links separated
b0d3d62
to
da4c7c8
Compare
- GIF movies for peptide examples - standard transformations - custom peptide 2 up transformation - updated code with nglview commands - use MDAnalysisData for runnable example
31ee785
to
baeedec
Compare
- performance - in-memory vs out of core
I added movies (see https://gist.github.com/orbeckst/ddbe912c665340791752c6622575cc6a for how they were made) and final remarks. I'll merge on Monday unless there are other comments. |
Adresses issue #86
Here's is the first version of the post. It demonstrates the "pure magic" (@jbarnoud) of the on-the-fly transformations.
I have some questions though: