Periodic Boundary Treatment in AtomGroup

[GoogleCodeExporter]

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)

Original issue reported on code.google.com by White.D....@gmail.com on 25 Oct 2013 at 7:47

• Blocked on: treatment of PBC in general triclinic boxes: general minimum image convention #136

[GoogleCodeExporter]

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 ).

Original comment by orbeckst on 25 Oct 2013 at 11:35

• Changed state: Accepted
• Now blocked on: treatment of PBC in general triclinic boxes: general minimum image convention #136

[GoogleCodeExporter]

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=ae1b14ced
a80d99008d59f49be44eac7cbd9148d&name=issue156

Original comment by richardjgowers on 26 Jan 2014 at 10:39

[GoogleCodeExporter]

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

Original comment by orbeckst on 26 Jan 2014 at 11:09

[GoogleCodeExporter]

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

Original comment by richardjgowers on 27 Jan 2014 at 12:43

[GoogleCodeExporter]

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=517bfa44e
dfbd03c5fee1680f06973792f3a99f8&name=issue156

Original comment by richardjgowers on 27 Jan 2014 at 9:47

[GoogleCodeExporter]

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.

Original comment by orbeckst on 3 Feb 2014 at 5:09

• Changed state: Fixed