diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index f7af77f7b8b..c75ba4aa536 100644 --- a/package/MDAnalysis/core/groups.py +++ b/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` in the group. @@ -666,6 +666,9 @@ def center(self, weights, pbc=None, compound='group'): be calculated without moving any :class:`Atoms` 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]) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index dec11c68ea7..f15349a772c 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/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) \ No newline at end of file + return sqrt(result)