Issue 1004 unwrap

package/MDAnalysis/core/groups.py

@@ -103,7 +103,7 @@
 
 from .. import _ANCHOR_UNIVERSES
 from ..lib import util
-from ..lib.util import cached, warn_if_not_unique, unique_int_1d
+from ..lib.util import cached, warn_if_not_unique, unique_int_1d, make_whole
 from ..lib import distances
 from ..lib import transformations
 from ..selections import get_writer as get_selection_writer_for
@@ -647,7 +647,7 @@ def isunique(self):
         return not np.count_nonzero(mask)
 
     @warn_if_not_unique
-    def center(self, weights, pbc=None, compound='group'):
+    def center(self, weights, pbc=None, unwrap=False, compound='group'):
         """Weighted center of (compounds of) the group
 
         Computes the weighted center of :class:`Atoms<Atom>` in the group.
@@ -666,6 +666,9 @@ def center(self, weights, pbc=None, compound='group'):
             be calculated without moving any :class:`Atoms<Atom>` to keep the
             compounds intact. Instead, the resulting position vectors will be
             moved to the primary unit cell after calculation.
+        unwrap : bool, optional
+            if True, atom positions are shifted so as to not break any bonds over
+            periodic boundaries.
         compound : {'group', 'segments', 'residues'}, optional
             If ``'group'``, the weighted center of all atoms in the group will
             be returned as a single position vector. Else, the weighted centers
@@ -1138,6 +1141,98 @@ def pack_into_box(self, box=None, inplace=True):
                 return ag.universe.trajectory.ts.positions[unique_ix[restore_mask]]
             return packed_coords[restore_mask]
 
+    def unwrap(self, compound='group', center=None, reference_point=None):
+        """Move coordinates to prevent bonds being split over periodic boundaries
+
+        This is the inverse transformation to packing atoms into the primary unit
+        cell. It should therefore stop bonds being "broken" over the periodic
+        boundaries in the system.
+
+        This method requires the Universe to have information about bonding.  This
+        can be guessed using the :meth:`guess_bonds` method if necessary.
+
+        Parameters
+        ----------
+        compound : str, optional {'group', 'residues', 'segments', 'fragments'}
+            which level of topology to keep together [``fragment``]
+        center : numpy.ndarray, optional
+            how to define center of group, 'com', 'cog'
+        reference_point : np.ndarray
+            position to try and center unwrapped molecules around.  If given
+            then for each group unwrapped, the COM of the atoms will be as close
+            as possible to this point.  Defaults to center of primary unit cell.
+
+        # TODO the function isn't actually here, merge the _cutil in to make docs work
+        .. seealso:: :func:`MDAnalysis.lib.util.make_whole`
+
+        .. versionadded:: 0.18.1
+        """
+        atomgroup = self.atoms.unique
+
+        if compound.lower() == 'group':
+            objects = [atomgroup]
+        elif compound.lower() == 'residues':
+            objects = atomgroup.residues
+        elif compound.lower() == 'segments':
+            objects = atomgroup.segments
+        elif compound.lower() == 'fragments':
+            objects = atomgroup.fragments
+        else:
+            raise ValueError("Unrecognized compound definition: {0}"
+                             "Please use one of 'group' 'residues' 'segments'"
+                             "or 'fragments'".format(compound))
+
+        if reference_point is None:
+            tri_box = mdamath.triclinic_vectors(self.dimensions)
+            center = np.diag(tri_box) / 2.0
+
+        for o in objects:
+            make_whole(o.atoms)
+            o.choose_image(reference_point, center=center, inplace=True)
+
+    def choose_image(self, reference_point, center='com', inplace=True):
+        """Shift coordinates to be as close as possible to reference point
+
+        Parameters
+        ----------
+        reference_point : coordinate
+          position to move group close to
+        center : str, optional
+          how to define the center of AtomGroup
+        inplace : bool, optional
+          if True, changes AtomGroup positions directly, otherwise returns
+          the positions [``True``]
+
+        Returns
+        -------
+        shifted_coords : np.ndarray
+          if not inplace, returns the coordinates of this Group which
+          minimise the distance from center of Group to reference_point
+
+        .. versionadded:: 0.19.0
+        """
+        ag = self.atoms.unique
+
+        if center.lower() in ('com', 'centerofmass'):
+            center = ag.center_of_mass()
+        elif center.lower() in ('cog', 'centroid', 'centerofgeometry'):
+            center = ag.center_of_geometry()
+        else:
+            raise ValueError("Unrecognized center definition: {0}"
+                             "Please use one of 'com' or 'cog'".format(center))
+
+        shift = reference_point - center
+        box_shifts = np.rint(shift)
+
+        if inplace:
+            if any(box_shifts):
+                ag.translate(s)
+        else:
+            if any(box_shifts):
+                return self.positions + shift
+            else:
+                return self.positions
+
     def wrap(self, compound="atoms", center="com", box=None):
         """Shift the contents of this group back into the unit cell.
 
@@ -1198,6 +1293,8 @@ def wrap(self, compound="atoms", center="com", box=None):
                              "Please use one of 'group' 'residues' 'segments'"
                              "or 'fragments'".format(compound))
 
+        for o in self.objects:
+
         # TODO: ADD TRY-EXCEPT FOR MASSES PRESENCE
         if center.lower() in ('com', 'centerofmass'):
             centers = np.vstack([o.atoms.center_of_mass() for o in objects])

package/MDAnalysis/lib/_cutil.pyx

@@ -98,7 +98,7 @@ cdef intset difference(intset a, intset b):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def make_whole(atomgroup, reference_atom=None):
+def make_whole(atomgroup, reference_atom=None, inplace=True):
     """Move all atoms in a single molecule so that bonds don't split over images
 
     Atom positions are modified in place.
@@ -125,9 +125,12 @@ def make_whole(atomgroup, reference_atom=None):
         The :class:`MDAnalysis.core.groups.AtomGroup` to work with.
         The positions of this are modified in place.  All these atoms
         must belong in the same molecule or fragment.
-    reference_atom : :class:`~MDAnalysis.core.groups.Atom`
+    reference_atom : :class:`~MDAnalysis.core.groups.Atom`, optional
         The atom around which all other atoms will be moved.
         Defaults to atom 0 in the atomgroup.
+    inplace : bool, optional
+        if True, directly modify the Atom's coordinates, if False the
+        unwrapped coordinates are instead returned [``True``]
 
     Raises
     ------
@@ -165,6 +168,10 @@ def make_whole(atomgroup, reference_atom=None):
 
 
     .. versionadded:: 0.11.0
+    .. versionchanged:: 0.18.1
+       Added support for triclinic boxes
+       Rewrote in Cython (better performance)
+       Added inplace keyword
     """
     cdef intset refpoints, todo, done
     cdef int i, nloops, ref, atom, other, natoms
@@ -270,7 +277,10 @@ def make_whole(atomgroup, reference_atom=None):
     if refpoints.size() < natoms:
         raise ValueError("AtomGroup was not contiguous from bonds, process failed")
     else:
-        atomgroup.positions = newpos
+        if inplace:
+            atomgroup.positions = newpos
+        else:
+            return newpos
 
 
 @cython.boundscheck(False)
@@ -310,4 +320,4 @@ cdef float _norm(float * a):
     result = 0.0
     for n in range(3):
         result += a[n]*a[n]
-    return sqrt(result)
+    return sqrt(result)