Issue 156 in mdanalysis: Periodic Boundary Treatment in AtomGroup

[mdana...@googlecode.com]

Status: New
Owner: ----
Labels: Type-Defect Priority-Medium

New issue 156 by White.D....@gmail.com: Periodic Boundary Treatment in
AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

I've noticed there's no accounting for periodic boundary conditions for the
AtomGroup methods. Many of the methods that rely on atom distances don't
use the minimum image convention. I'd like to put this as an issue, because
these implementations are not correct in cases when AtomGroups extend
across boundaries. I know it's very time consuming to fix these
implementations, but perhaps there is a way of warning the user or noting
in the documentation that they do not account for periodic boundaries.

Here's an example:

def centerOfMass(self):
"""Center of mass of the selection."""
return
numpy.sum(self.coordinates()*self.masses()[:,numpy.newaxis],axis=0)/self.totalMass()

This will not give the correct center of mass.

This code here, will not give the correct atom-atom distance if the atoms
are not in the same image:

def bond(self):
"""Returns the distance between atoms in a 2-atom group.

Distance between atoms 0 and 1::

0---1

.. Note::

Only makes sense for a :class:`AtomGroup` with exactly 2
:class:`Atom`; anything else will raise a
:exc:`ValueError`.

.. versionadded:: 0.7.3
"""
if len(self) != 2:
raise ValueError("distance computation only makes sense for
a group with exactly 2 atoms")
return numpy.linalg.norm(self[0].pos - self[1].pos)



--
You received this message because this project is configured to send all
issue notifications to this address.
You may adjust your notification preferences at:
https://code.google.com/hosting/settings

[mdana...@googlecode.com]

Updates:
Status: Accepted
Blockedon: mdanalysis:136

Comment #1 on issue 156 by orbeckst: Periodic Boundary Treatment in
AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

Technically these methods do the correct calculation: they take the
positions provided in the trajectory and calculate distances based on these
positions and nothing else. However, I agree with you that what you
describe is the behavior many users might expect.

Given that this might trip up users we should at least note this in the
documentation.

In my experience, for many people doing simulations of a single protein in
solution or in a membrane that's often not an issue because the
trajectories are pre-processed to be centered and fitted on the protein in
the simulation cell. Of course, that's no good if you are interested in
multiple proteins or a solution of small molecules...

If we had a general facility for calculating minimum images (see Issue 136)
then we could at least add code that takes PBC into account if a flag is
set. (Incidentally, such a flag would be a strong argument in favor of
leaving many methods such as centerOfMass() actual methods instead of
turning them into managed attributes; see also the thread
https://groups.google.com/d/msg/mdnalysis-devel/jxEns_KC0xY/kDdlDe0W_R4J ).

[Richard Gowers]

I made an AtomGroup method that moves an atom group into the simulation box which would fix your example of center of mass calculation
http://code.google.com/p/mdanalysis/issues/detail?id=153

This doesn't fix minimum image convention though.  I could write a general function in analysis that takes a 3d list of separations and then applies minimum image convention to this.  Other methods could then call this through some sort of optional flag [pbc= False] which I could add in to existing methods?  This would keep the current behaviour while adding an option for this new case?

[Oliver Beckstein]

On 27 Oct, 2013, at 11:44, Richard Gowers wrote:

> I made an AtomGroup method that moves an atom group into the simulation box which would fix your example of center of mass calculation
> http://code.google.com/p/mdanalysis/issues/detail?id=153
>
> This doesn't fix minimum image convention though. I could write a general function in analysis that takes a 3d list of separations and then applies minimum image convention to this. Other methods could then call this through some sort of optional flag [pbc= False] which I could add in to existing methods? This would keep the current behaviour while adding an option for this new case?

The cleanest approach would be to have a general minimum image calculation working http://code.google.com/p/mdanalysis/issues/detail?id=136 and then use a flag as you suggest, e.g.

ag.centerOfMass(pbc=False) # old behavior (and faster), default
ag.centerOfMass(pbc=True) # new behavior, uses minimum image distances

Oliver

> You received this message because you are subscribed to the Google Groups "MDnalysis-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to mdnalysis-dev...@googlegroups.com.
> To post to this group, send email to mdnalys...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/mdnalysis-devel/9a751686-885a-4495-9d45-053078391c9a%40googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.

--
Oliver Beckstein * orbe...@gmx.net
skype: orbeckst * orbe...@gmail.com

[Richard Gowers]

https://code.google.com/r/richardjgowers-applypbc/source/detail?r=d5cbfecff60d8dbf7e999cb21ca860680ce7e22e

This seems to do the job.  I added 2 functions to core.util (adding them to analysis.distances caused an error in loading the MDAnalysis package for me?), applyPBC makes coordinates always fall within the box, minimumImage makes vectors between atoms conform to MIC. 

I then made the AtomGroup methods use these functions if flagged to do so.

If this is what you were looking for then I can add tests and make it a bit more robust (and add to any other methods that might benefit from it?)

Richard

[Oliver Beckstein]

Hi Richard,

Thanks for looking at this.

On 28 Oct, 2013, at 09:33, Richard Gowers wrote:

> https://code.google.com/r/richardjgowers-applypbc/source/detail?r=d5cbfecff60d8dbf7e999cb21ca860680ce7e22e

I had a quick look.

I don't think that your solution with applyPBC and centerOfMass() gives the expected answer: the center of mass might still end up outside of the user's box. We can't assume that the user wants all atoms shifted. You'd need to (1) calculate the COM for the group in the primary unit cell and (2) then shift the COM back into the correct spot in the user's unit cell.

The bond distances + MI looks ok but I would really like a general and fast solution for MI instead of the simple orthorhombic boxes – MI calculations can really slow down analysis so I'd like this to be done in the best possible way if we can. If you look at Issue 136 http://code.google.com/p/mdanalysis/issues/detail?id=136 then you'll see that Rob McGibbon actually has highly optimized code in his own package MDtraj to do these calculations. He was willing to share this code, it's just a question of integrating it into MDAnalysis.

>
> This seems to do the job. I added 2 functions to core.util (adding them to analysis.distances caused an error in loading the MDAnalysis package for me?), applyPBC makes coordinates always fall within the box, minimumImage makes vectors between atoms conform to MIC.
>
> I then made the AtomGroup methods use these functions if flagged to do so.
>
> If this is what you were looking for then I can add tests and make it a bit more robust (and add to any other methods that might benefit from it?)

I'd rather have Issue 136 solved for once and for all because then we can start building on a solid foundation instead of putting on band-aids ;-). Once we have robust MI and distance calculations we can then add them in various parts of the code that will benefit. (On that note: PBS-aware KD-trees Issue 137 is also of interest...)

Thanks,
Oliver

[mdana...@googlecode.com]

Comment #2 on issue 156 by richardj...@gmail.com: Periodic Boundary

Treatment in AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

Ok I think I've fixed this issue.

There's now a function called applyPBC(coords, box) in core.distances that
applies periodic boundary conditions. (Triclinic too)

Lots of AtomGroup methods now have a pbc flag to move atoms inside the box
before doing their thing. This is false by default so default behaviour
isn't affected.

The bond method can now use the minimage=True flag to apply minimum image.
This flag does not have any relevance to small wizards.

packIntoBox now has an inplace flag, default False.

Finally, calc_bonds from the topology work now handles triclinic minimum
image too.

https://code.google.com/r/richardjgowers-topologylists/source/detail?r=ae1b14ceda80d99008d59f49be44eac7cbd9148d&name=issue156

[mdana...@googlecode.com]

Comment #3 on issue 156 by orbeckst: Periodic Boundary Treatment in
AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

Richard, sounds good. Do you also have unit tests?

How hard would it be to add a flag in core.flags that determines the
default behavior to use or not to use the PBC code?

Oliver

[mdana...@googlecode.com]

Comment #4 on issue 156 by richardj...@gmail.com: Periodic Boundary

Treatment in AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

I think I'll be able to get flags working for this later today.

Tests are there for the applyPBC function, but not for each AtomGroup
method which uses the flag. I might add some tests to check the methods
and flags are working properly together,

Would it be a good idea to add the 'pbc' flag to ag.get_positions so
that .coordinates() is also affected by this?
Ie.
mda.core.flags['use_pbc'] = True
ag.coordinates() # this now always returns coordinates shifted inside the
box

[mdana...@googlecode.com]

Comment #5 on issue 156 by richardj...@gmail.com: Periodic Boundary

Treatment in AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

Ok, flag added 'use_pbc' default is False. Setting to True is equivalent
to using pbc=True for all the methods.

There's also now tests for everything

https://code.google.com/r/richardjgowers-topologylists/source/detail?r=517bfa44edfbd03c5fee1680f06973792f3a99f8&name=issue156

[mdana...@googlecode.com]

Updates:
Status: Fixed
Owner: richardj...@gmail.com

Comment #6 on issue 156 by orbeckst: Periodic Boundary Treatment in
AtomGroup
http://code.google.com/p/mdanalysis/issues/detail?id=156

I merged Richard's code and added a few cosmetic fixes:

- I changed the keyword for AtomGroup.bond() from "minimage" to "pbc" for
consistency with all the other methods.
for notes on using math)

Please check that everything works as it should.