Molecules across periodic boundaries are not properly treated

[hsidky]

Expected behaviour

A residue or atomgroup containing atoms that straddle a periodic boundary should correctly unwrap the collection to one side or the other in order to properly compute the center of mass or moment of inertia.

Actual behaviour

Although pbc=True flags exist for a variety of methods that operate on residues and atom groups, this only maps atom coordinates back into the original periodic simulation box. However, if analyzing a trajectory, such as in my case a collection of molecules, some of which are located across a periodic boundary, this does nothing to solve the issue of correctly computing centers of mass and (importantly for me) the moment of inertia.

Code to reproduce the behaviour

To deal with this, I wrote a few helper functions in my analysis scripts, which I've reproduced below. This is simply one approach; if the bond/topology information is available it would potentially be easier to utilize that information.

It would be great to see something like this incorporated into MDAnalysis, in a much more efficient manner. For reference, H is the Parrinello-Rahman matrix, or triclinic_dimensions as referenced in MDAnalysis.

import MDAnalysis as mda

def roundf(x):
    if x >= 0:
        return np.floor( x + 0.5 )
    else:
        return np.ceil( x - 0.5 )

def apply_min_image(H,x):
    s = np.dot(la.inv(H),x)
    for i in range(0,3):
        s[i] -= isperiodic[i]*roundf(s[i])
    return np.dot(H,s)

def upwrap_molecule(A, H):
    for i in range(1, A.shape[0]):
        A[i,:] = apply_min_image(H, A[i,:] - A[i-1,:]) + A[i-1,:]

My approach above is simple: it linearly marches through the coordinates and computes the minimium image distance between an atom and the next one, which will reconstruct a molecule on the side of the first (reference) atom.

If something like this, or a better approach, can be incorporated into MDAnalysis in a more efficient manner, then great! All I need is really a way to properly take into account periodicity. Currently the COM and inertial tensor are incorrect for molecules whose atoms are across a periodic boundary.

Thanks!

Currently version of MDAnalysis:

0.15.0

[richardjgowers]

Heya,

Thanks for getting in touch. There is currently MDAnalysis.lib.mdamath.make_whole which unwraps molecules:

https://github.com/MDAnalysis/mdanalysis/blob/develop/package/MDAnalysis/lib/mdamath.py#L263

This works based upon the bonded structure of an AtomGroup, so you'd need to either use a topology with bonds or guess them (AtomGroup.guess_bonds). If the molecule you are unwrapping is very large, then it is safer to use bonded structures rather than just comparing all distances within a molecule. It doesn't do triclinic boxes currently, but that's an easy fix.

Thinking about it, something like atomgroup.center_of_mass(unwrap=True) would be best right?

[hsidky]

Thank you, I didn't find that when I searched the docs. Yes, it would certainly be nice to have something like what you suggest. However, it would also be nice to have it on other methods such as moment_of_inertia. Those are currently the two most important to me, but I imagine other people might benefit from similar behavior on other AtomGroup calls. It fact, it would be quite nice if there was a global flag similar to use_pbc for this. Thanks again!

[shuaijiang-ustc]

Hi, everyone, I am wondering if this molecules across boundaries issue is solved or not. I just test the moment_of_inertia method under mdanalysis version 1.1.1, the problem still exists. Thanks!