*Group.unwrap

package/CHANGELOG

@@ -27,6 +27,11 @@ Enhancements
     molecule or per fragment (PR #2182)
   * Atoms (AtomGroups) have a new fragindex (fragindices) property corresponding
     to the respective fragment (fragments) (PR #2182)
+  * Atom-, Residue-, and SegmentGroups have a new unwrap() method allowing to
+    make compounds (group/residues/segments/fragments/molecules) that are broken
+    across periodic boundaries whole (Issues #1004 and #1185, PR #2189)
+  * lib.mdamath.make_whole() now returns the resulting coordinates and is able
+    to operate out-of-place, i.e., without modifying atom coordinates (PR #2189)
 
 Changes
   * added official support for Python 3.7 (PR #1963)

package/MDAnalysis/core/groups.py

@@ -107,6 +107,7 @@
 from ..lib.util import cached, warn_if_not_unique, unique_int_1d
 from ..lib import distances
 from ..lib import transformations
+from ..lib import mdamath
 from ..selections import get_writer as get_selection_writer_for
 from . import selection
 from . import flags
@@ -654,7 +655,9 @@ def center(self, weights, pbc=None, compound='group'):
         Computes the weighted center of :class:`Atoms<Atom>` in the group.
         Weighted centers per :class:`Residue`, :class:`Segment`, molecule, or
         fragment can be obtained by setting the `compound` parameter
-        accordingly.
+        accordingly. If the weights of a compound sum up to zero, the
+        coordinates of that compound's weighted center will be ``nan`` (not a
+        number).
 
         Parameters
         ----------
@@ -1211,6 +1214,12 @@ def wrap(self, compound="atoms", center="com", box=None):
            :meth:`wrap` with all default keywords is identical to
            :meth:`pack_into_box`
 
+        See Also
+        --------
+        :meth:`pack_into_box`
+        :meth:`unwrap`
+        :meth:`~MDanalysis.lib.distances.apply_PBC`
+
 
         .. versionadded:: 0.9.2
         """
@@ -1253,6 +1262,160 @@ def wrap(self, compound="atoms", center="com", box=None):
             if not all(s == 0.0):
                 o.atoms.translate(s)
 
+    def unwrap(self, compound='fragments', reference='com', inplace=True):
+        """Move atoms of this group so that bonds within the
+        group's compounds aren't split across periodic boundaries.
+
+        This function is most useful when atoms have been packed into the
+        primary unit cell, causing breaks mid-molecule, with the molecule then
+        appearing on either side of the unit cell. This is problematic for
+        operations such as calculating the center of mass of the molecule. ::
+
+           +-----------+       +-----------+
+           |           |       |           |
+           | 6       3 |       |         3 | 6
+           | !       ! |       |         ! | !
+           |-5-8   1-2-|  ==>  |       1-2-|-5-8
+           | !       ! |       |         ! | !
+           | 7       4 |       |         4 | 7
+           |           |       |           |
+           +-----------+       +-----------+
+
+        Parameters
+        ----------
+        compound : {'group', 'segments', 'residues', 'molecules', \
+                    'fragments'}, optional
+            Which type of component to unwrap. Note that, in any case, all
+            atoms within each compound must be interconnected by bonds, i.e.,
+            compounds must correspond to (parts of) molecules.
+        reference : {'com', 'cog', None}, optional
+            If ``'com'`` (center of mass) or ``'cog'`` (center of geometry), the
+            unwrapped compounds will be shifted so that their individual
+            reference point lies within the primary unit cell.
+            If ``None``, no such shift is performed.
+        inplace : bool, optional
+            If ``True``, coordinates are modified in place.
+
+        Returns
+        -------
+        coords : numpy.ndarray
+            Unwrapped atom coordinate array of shape ``(n, 3)``.
+
+        Raises
+        ------
+        NoDataError
+            If `compound` is ``'molecules'`` but the underlying topology does
+            not contain molecule information, or if `reference` is ``'com'``
+            but the topology does not contain masses.
+        ValueError
+            If `reference` is not one of ``'com'``, ``'cog'``, or ``None``, or
+            if `reference` is ``'com'`` and the total mass of any `compound` is
+            zero.
+
+        Note
+        ----
+        Be aware of the fact that only atoms *belonging to the group* will
+        be unwrapped! If you want entire molecules to be unwrapped, make sure
+        that all atoms of these molecules are part of the group.\n
+        An AtomGroup containing all atoms of all fragments in the group ``ag``
+        can be created with::
+
+          all_frag_atoms = sum(ag.fragments)
+
+
+        See Also
+        --------
+        :func:`~MDAnalysis.lib.mdamath.make_whole`,
+        :meth:`wrap`,
+        :meth:`pack_into_box`,
+        :func:`~MDanalysis.lib.distances.apply_PBC`
+
+
+        .. versionadded:: 0.20.0
+        """
+        atoms = self.atoms
+        # bail out early if no bonds in topology:
+        if not hasattr(atoms, 'bonds'):
+            raise NoDataError("{}.unwrap() not available; this requires Bonds"
+                              "".format(self.__class__.__name__))
+        unique_atoms = atoms.unique
+        if reference is not None:
+            ref = reference.lower()
+            if ref  == 'com':
+                if not hasattr(self, 'masses'):
+                    raise NoDataError("Cannot perform unwrap with "
+                                      "reference='com', this requires masses.")
+            elif ref != 'cog':
+                raise ValueError("Unrecognized reference '{}'. Please use one "
+                                 "of 'com', 'cog', or None.".format(reference))
+        comp = compound.lower()
+        if comp not in ('fragments', 'group', 'residues', 'segments',
+                        'molecules'):
+            raise ValueError("Unrecognized compound definition: {}\nPlease use"
+                             " one of 'group', 'residues', 'segments', "
+                             "'molecules', or 'fragments'.".format(compound))
+        # The 'group' needs no splitting:
+        if comp == 'group':
+            positions = mdamath.make_whole(unique_atoms, inplace=False)
+            # Apply reference shift if required:
+            if reference is not None and len(positions) > 0:
+                if ref == 'com':
+                    masses = unique_atoms.masses
+                    total_mass = masses.sum()
+                    if np.isclose(total_mass, 0.0):
+                        raise ValueError("Cannot perform unwrap with "
+                                         "reference='com' because the total "
+                                         "mass of the group is zero.")
+                    refpos = np.sum(positions * masses[:, None], axis=0)
+                    refpos /= total_mass
+                else:  # ref == 'cog'
+                    refpos = positions.mean(axis=0)
+                target = distances.apply_PBC(refpos, self.dimensions)
+                positions += target - refpos
+        # We need to split the group into compounds:
+        else:
+            if comp == 'fragments':
+                compound_indices = unique_atoms.fragindices
+            elif comp == 'residues':
+                compound_indices = unique_atoms.resindices
+            elif comp == 'segments':
+                compound_indices = unique_atoms.segindices
+            else:  # comp == 'molecules'
+                try:
+                    compound_indices = unique_atoms.molnums
+                except AttributeError:
+                    raise NoDataError("Cannot use compound='molecules': "
+                                      "No molecule information in topology.")
+            # Now process every compound:
+            unique_compound_indices = unique_int_1d(compound_indices)
+            positions = unique_atoms.positions
+            for i in unique_compound_indices:
+                mask = np.where(compound_indices == i)
+                c = unique_atoms[mask]
+                positions[mask] = mdamath.make_whole(c, inplace=False)
+                # Apply reference shift if required:
+                if reference is not None:
+                    if ref == 'com':
+                        masses = c.masses
+                        total_mass = masses.sum()
+                        if np.isclose(total_mass, 0.0):
+                            raise ValueError("Cannot perform unwrap with "
+                                             "reference='com' because the "
+                                             "total mass of at least one of "
+                                             "its {} is zero.".format(comp))
+                        refpos = np.sum(positions[mask] * masses[:, None],
+                                        axis=0)
+                        refpos /= total_mass
+                    else:  # ref == 'cog'
+                        refpos = positions[mask].mean(axis=0)
+                    target = distances.apply_PBC(refpos, self.dimensions)
+                    positions[mask] += target - refpos
+        if inplace:
+            unique_atoms.positions = positions
+        if not atoms.isunique:
+            positions = positions[atoms._unique_restore_mask]
+        return positions
+
     def copy(self):
         """Get another group identical to this one.
 
@@ -1842,11 +2005,11 @@ def __getattr__(self, attr):
         #
         # is this a known attribute failure?
         # TODO: Generalise this to cover many attributes
-        if attr in ('fragments', 'fragindices', 'n_fragments'):
+        if attr in ('fragments', 'fragindices', 'n_fragments', 'unwrap'):
             # eg:
             # if attr in _ATTR_ERRORS:
             # raise NDE(_ATTR_ERRORS[attr])
-            raise NoDataError("AtomGroup has no {}; this requires Bonds"
+            raise NoDataError("AtomGroup.{} not available; this requires Bonds"
                               "".format(attr))
         elif hasattr(self.universe._topology, 'names'):
             # Ugly hack to make multiple __getattr__s work

package/MDAnalysis/core/topologyattrs.py

@@ -775,7 +775,9 @@ def center_of_mass(group, pbc=None, compound='group'):
         Computes the center of mass of :class:`Atoms<Atom>` in the group.
         Centers of mass per :class:`Residue`, :class:`Segment`, molecule, or
         fragment can be obtained by setting the `compound` parameter
-        accordingly.
+        accordingly. If the masses of a compound sum up to zero, the
+        center of mass coordinates of that compound will be ``nan`` (not a
+        number).
 
         Parameters
         ----------
@@ -1457,7 +1459,7 @@ class Molnums(ResidueAttr):
     attrname = 'molnums'
     singular = 'molnum'
     target_classes = [AtomGroup, ResidueGroup, Atom, Residue]
-    dtype = int
+    dtype = np.int64
 
 # segment attributes
 

package/MDAnalysis/lib/_cutil.pyx

@@ -25,7 +25,7 @@
 import cython
 import numpy as np
 cimport numpy as np
-from libc.math cimport sqrt
+from libc.math cimport sqrt, fabs
 
 from MDAnalysis import NoDataError
 
@@ -103,10 +103,9 @@ cdef intset difference(intset a, intset b):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def make_whole(atomgroup, reference_atom=None):
-    """Move all atoms in a single molecule so that bonds don't split over images
-
-    Atom positions are modified in place.
+def make_whole(atomgroup, reference_atom=None, inplace=True):
+    """Move all atoms in a single molecule so that bonds don't split over
+    images.
 
     This function is most useful when atoms have been packed into the primary
     unit cell, causing breaks mid molecule, with the molecule then appearing
@@ -133,6 +132,13 @@ def make_whole(atomgroup, reference_atom=None):
     reference_atom : :class:`~MDAnalysis.core.groups.Atom`
         The atom around which all other atoms will be moved.
         Defaults to atom 0 in the atomgroup.
+    inplace : bool, optional
+        If ``True``, coordinates are modified in place.
+
+    Returns
+    -------
+    coords : numpy.ndarray
+        The unwrapped atom coordinates.
 
     Raises
     ------
@@ -155,12 +161,12 @@ def make_whole(atomgroup, reference_atom=None):
         # This algorithm requires bonds, these can be guessed!
         u = mda.Universe(......, guess_bonds=True)
 
-        # MDAnalysis can split molecules into their fragments
+        # MDAnalysis can split AtomGroups into their fragments
         # based on bonding information.
         # Note that this function will only handle a single fragment
         # at a time, necessitating a loop.
-        for frag in u.fragments:
-          make_whole(frag)
+        for frag in u.atoms.fragments:
+            make_whole(frag)
 
     Alternatively, to keep a single atom in place as the anchor::
 
@@ -169,24 +175,41 @@ def make_whole(atomgroup, reference_atom=None):
         make_whole(atomgroup, reference_atom=atomgroup[10])
 
 
+    See Also
+    --------
+    :meth:`MDAnalysis.core.groups.AtomGroup.unwrap`
+
+
     .. versionadded:: 0.11.0
+    .. versionchanged:: 0.20.0
+        Inplace-modification of atom positions is now optional, and positions
+        are returned as a numpy array.
     """
     cdef intset refpoints, todo, done
-    cdef int i, nloops, ref, atom, other, natoms
+    cdef int i, j, nloops, ref, atom, other, natoms
     cdef cmap[int, int] ix_to_rel
     cdef intmap bonding
     cdef int[:, :] bonds
     cdef float[:, :] oldpos, newpos
     cdef bint ortho
     cdef float[:] box
     cdef float tri_box[3][3]
+    cdef float half_box[3]
     cdef float inverse_box[3]
     cdef double vec[3]
     cdef ssize_t[:] ix_view
+    cdef bint is_unwrapped
 
     # map of global indices to local indices
     ix_view = atomgroup.ix[:]
     natoms = atomgroup.ix.shape[0]
+
+    oldpos = atomgroup.positions
+
+    # Nothing to do for less than 2 atoms
+    if natoms < 2:
+        return np.array(oldpos)
+
     for i in range(natoms):
         ix_to_rel[ix_view[i]] = i
 
@@ -198,9 +221,11 @@ def make_whole(atomgroup, reference_atom=None):
             raise ValueError("Reference atom not in atomgroup")
         ref = ix_to_rel[reference_atom.ix]
 
-    box = atomgroup.dimensions
+    box = atomgroup.dimensions.astype(np.float32)
+    # TODO: remove astype(np.float32) once all universes return float32 boxes
 
     for i in range(3):
+        half_box[i] = 0.5 * box[i]
         if box[i] == 0.0:
             raise ValueError("One or more dimensions was zero.  "
                              "You can set dimensions using 'atomgroup.dimensions='")
@@ -211,6 +236,17 @@ def make_whole(atomgroup, reference_atom=None):
             ortho = False
 
     if ortho:
+        # If atomgroup is already unwrapped, bail out
+        is_unwrapped = True
+        for i in range(1, natoms):
+            for j in range(3):
+                if fabs(oldpos[i, j] - oldpos[0, j]) >= half_box[j]:
+                    is_unwrapped = False
+                    break
+            if not is_unwrapped:
+                break
+        if is_unwrapped:
+            return np.array(oldpos)
         for i in range(3):
             inverse_box[i] = 1.0 / box[i]
     else:
@@ -233,7 +269,6 @@ def make_whole(atomgroup, reference_atom=None):
             bonding[atom].insert(other)
             bonding[other].insert(atom)
 
-    oldpos = atomgroup.positions
     newpos = np.zeros((oldpos.shape[0], 3), dtype=np.float32)
 
     refpoints = intset()  # Who is safe to use as reference point?
@@ -274,8 +309,9 @@ def make_whole(atomgroup, reference_atom=None):
 
     if refpoints.size() < natoms:
         raise ValueError("AtomGroup was not contiguous from bonds, process failed")
-    else:
+    if inplace:
         atomgroup.positions = newpos
+    return np.array(newpos)
 
 
 @cython.boundscheck(False)