Add feature somewhere to remove periodic boundary conditions for trajectory

[sryckbos]

Not quite sure how to implement this, but having a feature analogous to that of gromacs, where you use trjconv -pbc whole. Since you have to specify a topology anyway when doing md.load, all the information to do this should be available. The problem is just figuring out a way to do it and implementing it.

[rmcgibbo]

Yes, +1 for adding this feature. I think it's also on the issue tracker as #490.

[kyleabeauchamp]

PS: I think VMD has a nice TCL library called pbctools which could serve as a template for various ways for us to implement this. The key ingredient is md.compute_displacements().

[rmcgibbo]

For the center of mass in PBCs, there's also http://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions

[kyleabeauchamp]

At this point, the easiest way to implement this might be via a round-trip journey to pytraj and back.

[gph82]

Pinging the mdtraj crew...can you check markovmodel/PyEMMA#732 (in particular markovmodel/PyEMMA#732 (comment)) and give us your thoughts? I think this issues are somehow related: It is about the "problem of using an unwrapped trajectory with pbc information as an input? " Thanks

[rmcgibbo]

It would be easy to add a bool periodic=True kwarg to compute_contacts. I
think otherwise this is solved by #1058
?

On Thu, Mar 17, 2016 at 10:10 AM, Guillermo Pérez-Hernández <
notifications@github.com> wrote:

Pinging the mdtraj crew...can you check markovmodel/PyEMMA#732
markovmodel/PyEMMA#732 (in particular markovmodel/PyEMMA#732
(comment)
markovmodel/PyEMMA#732 (comment))
and give us your thoughts? I think this issues are somehow related: It is
about the "problem of using an unwrapped trajectory with pbc information as
an input? " Thanks

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[franknoe]

May I'm missing something here, but I think this is a separate issue. #1058 appears to implement a function for re-centering molecules in a box and applying pbc's (wrapping). markovmodel/PyEMMA#732 has two aspects:

1. add bool periodic=True to compute_contacts as Robert suggests (yes, that would help)
2. deal with inconsistent or misleading trajectory input. The problem arises when you have an already unwrapped trajectory (e.g. by using VMD), but the trajectory still has pbc information. Passing such a trajectory to compute_distances with the default setting periodic=True or to comput_contacts (which currently has no periodic), will give unphysical distances. So a possible solution would be to either catch loading of a trajectory that has pbc information but atoms outside the box and warn the user (if this case is considered to be inconsistent input), or alternatively deal with this at the level of compute_distances / compute_contacts, and set back to periodic=False while issuing a warning.

[rmcgibbo]

Do you have any suggestions on how to efficiently find and remedy
"inconsistent or misleading trajectory input"?

For PDBs, we now have an option (which is on by default) that checks the
density when loading PDB files with PBCs, and ignores the PBCs if the
density is crazy-unphysically high. The reason for this is mostly that a
lot of NMR models in the PDB contain CRYST1 records with a 1Ax1Ax1A box
which is supposed to be ignored. That's a pretty good use case, but
generally I'm a little skeptical of trying to be clever about deducing what
the user meant to do -- one person's "inconsistent or misleading"
trajectory might be precisely what someone else is interested in, no?

Adding bool periodic=True to compute_contacts is an easy fix. If someone
wants to submit a PR for that, it should only take a 2-line change and I'd
be happy to merge it immediately.

-Robert

On Thu, Mar 17, 2016 at 11:41 AM, Frank Noe notifications@github.com
wrote:

May I'm missing something here, but I think this is a separate issue.
#1058 #1058 appears to implement a
function for re-centering molecules in a box and applying pbc's (wrapping).
markovmodel/PyEMMA#732 markovmodel/PyEMMA#732
has two aspects:

1. add bool periodic=True to compute_contacts as Robert suggests (yes,
   that would help)
2. deal with inconsistent or misleading trajectory input. The problem
   arises when you have an already unwrapped trajectory (e.g. by using VMD),
   but the trajectory still has pbc information. Passing such a trajectory to
   compute_distances with the default setting periodic=True or to
   comput_contacts (which currently has no periodic), will give unphysical
   distances. So a possible solution would be to either catch loading of a
   trajectory that has pbc information but atoms outside the box and warn the
   user (if this case is considered to be inconsistent input), or
   alternatively deal with this at the level of compute_distances /
   compute_contacts, and set back to periodic=False while issuing a warning.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[franknoe]

Am 17/03/16 um 16:52 schrieb Robert T. McGibbon:

Do you have any suggestions on how to efficiently find and remedy
"inconsistent or misleading trajectory input"?

For PDBs, we now have an option (which is on by default) that checks the
density when loading PDB files with PBCs, and ignores the PBCs if the
density is crazy-unphysically high. The reason for this is mostly that a
lot of NMR models in the PDB contain CRYST1 records with a 1Ax1Ax1A box
which is supposed to be ignored. That's a pretty good use case, but
generally I'm a little skeptical of trying to be clever about deducing
what
the user meant to do -- one person's "inconsistent or misleading"
trajectory might be precisely what someone else is interested in, no?
I agree, it's difficult to be clever about this when loading the
trajectory. However,
I think at the level compute_distances or compute_contacts it seems more
easy to
decide. When you get distances of 0.1 Angstrom, something obviously went
wrong.
I think when distances are computed for data which is not wrapped inside
the box,
peridic=True should not be applied, or at least one could issue a warning.
If it's efficient to check if there are data outside the box boundaries,
that check
could be made in compute_distances or compute_contacts.

I'm happy to check this in PyEMMA and issue a warning there, but I have
the feeling
that mdtraj would be the better place for this.

Adding bool periodic=True to compute_contacts is an easy fix. If someone
wants to submit a PR for that, it should only take a 2-line change and I'd
be happy to merge it immediately.
ok. @marscher or @gph82, can you do that?

-Robert

On Thu, Mar 17, 2016 at 11:41 AM, Frank Noe notifications@github.com
wrote:

May I'm missing something here, but I think this is a separate issue.
#1058 #1058 appears to
implement a
function for re-centering molecules in a box and applying pbc's
(wrapping).
markovmodel/PyEMMA#732
markovmodel/PyEMMA#732
has two aspects:

1. add bool periodic=True to compute_contacts as Robert suggests (yes,
   that would help)
2. deal with inconsistent or misleading trajectory input. The problem
   arises when you have an already unwrapped trajectory (e.g. by using
   VMD),
   but the trajectory still has pbc information. Passing such a
   trajectory to
   compute_distances with the default setting periodic=True or to
   comput_contacts (which currently has no periodic), will give unphysical
   distances. So a possible solution would be to either catch loading of a
   trajectory that has pbc information but atoms outside the box and
   warn the
   user (if this case is considered to be inconsistent input), or
   alternatively deal with this at the level of compute_distances /
   compute_contacts, and set back to periodic=False while issuing a
   warning.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

---

Prof. Dr. Frank Noe
Head of Computational Molecular Biology group
Freie Universitaet Berlin

Phone: (+49) (0)30 838 75354
Web: research.franknoe.de

Mail: Arnimallee 6, 14195 Berlin, Germany

[gph82]

See PR #1072

[rmcgibbo]

I'm not sure that checking that distances come out greater than 0.1A (and
issuing a warning if they don't) is a great idea. You might have Drudes or
something. I think this kind of check is more appropriate in application
code than in mdtraj.

On Thu, Mar 17, 2016 at 12:24 PM, Guillermo Pérez-Hernández <
notifications@github.com> wrote:

See PR #1072 #1072

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[franknoe]

You're right. I would also prefer to check whether there are out-of-box
coordinates in conjunction with periodic=True. Wouldn't you agree that
this is a problematic case? I think periodic=True implies that your
coordinates
should adhere to the box boundaries.

We can do distance checks in PyEMMA, but this is probably less efficient
than checking for out-of-box coordinates at mdtraj level (assuming this
can be
done on the C/C++ level).

Is there an efficient way to check if any coordinates are ouf of the box
boundaries
/ are unwrapped in mdtraj? If that exists, we can do it on the analysis
level and prevent
this issue from happening.

Am 17/03/16 um 18:05 schrieb Robert T. McGibbon:

I'm not sure that checking that distances come out greater than 0.1A (and
issuing a warning if they don't) is a great idea. You might have Drudes or
something. I think this kind of check is more appropriate in application
code than in mdtraj.

On Thu, Mar 17, 2016 at 12:24 PM, Guillermo Pérez-Hernández <
notifications@github.com> wrote:

See PR #1072 #1072

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

---

Prof. Dr. Frank Noe
Head of Computational Molecular Biology group
Freie Universitaet Berlin

Phone: (+49) (0)30 838 75354
Web: research.franknoe.de

Mail: Arnimallee 6, 14195 Berlin, Germany

[rmcgibbo]

What exactly is an out-of-box coordinate under PBCs? For a cubic box, for
example, there's no reason to constrain the positions to be in
[0,boxLenth). Some people want the box to be centered at (0,0,0) or others
at (boxLength/2, boxLength/2, boxLength/2).

On Thu, Mar 17, 2016 at 1:35 PM, Frank Noe notifications@github.com wrote:

You're right. I would also prefer to check whether there are out-of-box
coordinates in conjunction with periodic=True. Wouldn't you agree that
this is a problematic case? I think periodic=True implies that your
coordinates
should adhere to the box boundaries.

We can do distance checks in PyEMMA, but this is probably less efficient
than checking for out-of-box coordinates at mdtraj level (assuming this
can be
done on the C/C++ level).

Is there an efficient way to check if any coordinates are ouf of the box
boundaries
/ are unwrapped in mdtraj? If that exists, we can do it on the analysis
level and prevent
this issue from happening.

Am 17/03/16 um 18:05 schrieb Robert T. McGibbon:

I'm not sure that checking that distances come out greater than 0.1A (and
issuing a warning if they don't) is a great idea. You might have Drudes
or
something. I think this kind of check is more appropriate in application
code than in mdtraj.

On Thu, Mar 17, 2016 at 12:24 PM, Guillermo Pérez-Hernández <
notifications@github.com> wrote:

See PR #1072 #1072

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

---

Prof. Dr. Frank Noe
Head of Computational Molecular Biology group
Freie Universitaet Berlin

Phone: (+49) (0)30 838 75354
Web: research.franknoe.de

Mail: Arnimallee 6, 14195 Berlin, Germany

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[nsplattner]

These different PBC-conventions are actually argument for not using periodic=True as a default. I think the PBC options and storage formats are not sufficiently standardized for this to work correctly in all cases and its much better if the user sets the option explicitly.

[rmcgibbo]

I don't really understand. If you have periodic boundary conditions, the
distances are invariant to translations of the unit cell. These different
conventions are not meaningful, but they serve to illustrate the problem
with simply looking at the value of the cartesian positions to determine
whether PBCs are in effect or not.

If you don't want to use periodic boundary conditions or the minimum image
convention, the best thing with MDTraj is simply to not set any periodic
box vectors. If your trajectory has periodic box vectors associated with
it, it's assumed that this is what you meant to do, and distances will be
computed taking them into account by default.

On Thu, Mar 17, 2016 at 1:48 PM, nsplattner notifications@github.com
wrote:

These different PBC-conventions are actually argument for not using
periodic=True as a default. I think the PBC options and storage formats are
not sufficiently standardized for this to work correctly in all cases and
its much better if the user sets the option explicitly.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[rmcgibbo]

(or my_trajectory.unitcell_vectors = None if you just want to remove
them from an already-loaded trajectory)

On Thu, Mar 17, 2016 at 1:52 PM, Robert T. McGibbon rmcgibbo@gmail.com
wrote:

I don't really understand. If you have periodic boundary conditions, the
distances are invariant to translations of the unit cell. These different
conventions are not meaningful, but they serve to illustrate the problem
with simply looking at the value of the cartesian positions to determine
whether PBCs are in effect or not.

If you don't want to use periodic boundary conditions or the minimum image
convention, the best thing with MDTraj is simply to not set any periodic
box vectors. If your trajectory has periodic box vectors associated with
it, it's assumed that this is what you meant to do, and distances will be
computed taking them into account by default.

On Thu, Mar 17, 2016 at 1:48 PM, nsplattner notifications@github.com
wrote:

These different PBC-conventions are actually argument for not using
periodic=True as a default. I think the PBC options and storage formats are
not sufficiently standardized for this to work correctly in all cases and
its much better if the user sets the option explicitly.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

-Robert

[nsplattner]

This is technically o.k., but the case where this becomes problematic is the following: the PBC information stored in a trajectory is not necessarily meaningful and I think most users are not aware that its provided and affects the analysis. In my case I was using VMD to wrap the trajectories and I think the PBC information was just copied from the original trajectory, but is obviously meaningless when used in combination with the already wrapped coordinates. I don't think thats the only problematic case since most people work with processed trajectories (e.g. waters removed, coordinates aligned). In these cases the program used for this needs to store the PBC information in a way that its not wrongly interpreted and I think there is no convention for this.

[rmcgibbo]

Why not just remove the periodic box vectors from the trajectory if you've post-processed the trajectory in such a way such that you no longer want to interpret the coordinates with periodic boundary conditions?

On Thu, Mar 17, 2016 at 2:04 PM, nsplattner notifications@github.com
wrote:

This is technically o.k., but the case where this becomes problematic is
the following: the PBC information stored in a trajectory is not
necessarily meaningful and I think most users are not aware that its
provided and affects the analysis. In my case I was using VMD to wrap the
trajectories and I think the PBC information was just copied from the
original trajectory, but is obviously meaningless when used in combination
with the already wrapped coordinates. I don't think thats the only
problematic case since most people work with processed trajectories (e.g.
waters removed, coordinates aligned). In these cases the program used for
this needs to store the PBC information in a way that its not wrongly
interpreted and I think there is no convention for this.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

-Robert

[mpharrigan]

I don't understand: even if you've wrapped your trajectory with trjconv or something, the box vectors still have meaning and distances still make sense under minimum image. If your minimum image distances are interacting with a copy of the system, doesn't that just mean your box is too small?

[rmcgibbo]

#1072 is merged, so I'm going to close this issue. If there is a remaining problem feel free to open another issue.

[nsplattner]

I'm not sure why the information stored in the trajectories wrapped by VMD is wrong or has the wrong effect, but its clearly not a problem of interacting copies in the system. And the problem is solved when using periodic=False. So clearly this option is doing something wrong.

[franknoe]

Matthew - I think you're right. We need to check this. Thanks.

Am 17/03/16 um 19:08 schrieb Matthew Harrigan:

I don't understand: even if you've wrapped your trajectory with
trjconv or something, the box vectors still have meaning and distances
still make sense under minimum image. If your minimum image distances
are interacting with a copy of the system, doesn't that just mean your
box is too small?

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

---

Prof. Dr. Frank Noe
Head of Computational Molecular Biology group
Freie Universitaet Berlin

Phone: (+49) (0)30 838 75354
Web: research.franknoe.de

Mail: Arnimallee 6, 14195 Berlin, Germany

[nsplattner]

O.k., so probably this is a VMD-issue and not an mdtraj-issue. Thanks for the comments.

[franknoe]

Well not necessarily. Matthew's point is: If the original box size is 4
nanometers, we unwrap the coordinates and store (for some reason) an
incorrect box size of 3 nanometers, then computing distances with
periodic=False would work while periodic=True wouldn't (because we
re-wrap the wrong way). If, however we store the correct box size in the
unwrapped coordinates, then the argument periodic should make no
difference. If it would, then it would be an mdtraj issue.

Am 17/03/16 um 19:12 schrieb nsplattner:

I'm not sure why the information stored in the trajectories wrapped by
VMD is wrong or has the wrong effect, but its clearly not a problem of
interacting copies in the system. And the problem is solved when using
periodic=False. So clearly this option is doing something wrong.

—
You are receiving this because you commented.
Reply to this email directly or view it on GitHub
#787 (comment)

---

Prof. Dr. Frank Noe
Head of Computational Molecular Biology group
Freie Universitaet Berlin

Phone: (+49) (0)30 838 75354
Web: research.franknoe.de

Mail: Arnimallee 6, 14195 Berlin, Germany