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