Search code examples
mdanalysis

MDAnalysis selects atoms from the PBC box but does not shift the coordinates

MDAnalysis distance selection commands like 'around' and 'sphzere' selects atoms from periodic image (I am using a rectangular box). 

universe.select_atoms("name OW and around 4 (resid 20 and name O2)")

However, the coordinates of the atoms from the PBC box reside on the other side of the box. In other words, I have to manually translate the atoms to ensure that they actually are withing the 4 Angstrom distance.

Is there a selection feature to achieve this using the select_atoms function?

Solution

  • If I well understand, you would like to get the atoms around a given selection in the image that is the closest to that selection.

    universe.select_atoms does not modify the coordinates, and I am not aware of a function that gives you what you want. The following function could work for an orthorhombic box like yours:

    def pack_around(atom_group, center):
    """
    Translate atoms to their periodic image the closest to a given point.

    The function assumes that the center is in the main periodic image.
    """
    # Get the box for the current frame
    box = atom_group.universe.dimensions
    # The next steps assume that all the atoms are in the same
    # periodic image, so let's make sure it is the case
    atom_group.pack_into_box()
    # AtomGroup.positions is a property rather than a simple attribute.
    # It does not always propagate changes very well so let's work with
    # a copy of the coordinates for now.
    positions = atom_group.positions.copy()
    # Identify the *coordinates* to translate.
    sub = positions - center
    culprits = numpy.where(numpy.sqrt(sub**2) > box[:3] / 2)
    # Actually translate the coordinates.
    positions[culprits] -= (u.dimensions[culprits[1]]
                            * numpy.sign(sub[culprits]))
    # Propagate the new coordinates.
    atom_group.positions = positions

    Using that function, I got the expected behavior on one of MDAnalysis test files. You need MDAnalysisTests to be installed to run the following piece of code:

    import numpy
import MDAnalysis as mda
from MDAnalysisTests.datafiles import PDB_sub_sol

u = mda.Universe(PDB_sub_sol)
selection = u.select_atoms('around 15 resid 32')
center = u.select_atoms('resid 32').center_of_mass()

# Save the initial file for latter comparison
u.atoms.write('original.pdb')
selection.write('selection_original.pdb')

# Translate the coordinates
pack_around(selection, center)

# Save the new coordinates
u.atoms.write('modified.pdb')
selection.write('selection_modified.pdb')