From 09d5774cd280f74cb48e4d0813d6bd1a352095ae Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 7 Jul 2018 11:23:00 -0500 Subject: [PATCH 1/5] added unwrap method to Groups --- package/MDAnalysis/core/groups.py | 65 ++++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index f7af77f7b8b..a51091c1df2 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 @@ -1138,6 +1138,69 @@ 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, box=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 + 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. + box : array_like + Box dimensions, can be either orthogonal or triclinic information. + Cell dimensions must be in an identical to format to those returned + by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`, + ``[lx, ly, lz, alpha, beta, gamma]``. If ``None``, uses these + timestep dimensions. + + # 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 + """ + if compound.lower() == 'group': + objects = [atomgroup.atoms] + 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)) + + for o in objects: + make_whole(o) + + if center is None: + tri_box = mdamath.triclinic_vectors(u.dimensions) + center = np.diag(tri_box) / 2.0 + + object_centers = np.vstack([o.center_of_mass(pbc=False) for o in objects]) + + if box is None: + box = self.dimensions + + dests = distances.apply_PBC(object_centers, box=box) + shifts = dests - object_centers + + for o, s in zip(objects, shifts): + # Save some needless shifts + if not all(s == 0.0): + o.atoms.translate(s) + def wrap(self, compound="atoms", center="com", box=None): """Shift the contents of this group back into the unit cell. From 27436f4f280d420a8b8da4db5f3e555775869ba0 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 7 Jul 2018 15:05:55 -0500 Subject: [PATCH 2/5] WIP of unwrap method --- package/MDAnalysis/core/groups.py | 8 +++++--- package/MDAnalysis/lib/_cutil.pyx | 11 +++++++++-- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index a51091c1df2..8fc6b9edf95 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1168,8 +1168,10 @@ def unwrap(self, compound='group', center=None, box=None): .. versionadded:: 0.18.1 """ + atomgroup = self.atoms.unique + if compound.lower() == 'group': - objects = [atomgroup.atoms] + objects = [atomgroup] elif compound.lower() == 'residues': objects = atomgroup.residues elif compound.lower() == 'segments': @@ -1182,13 +1184,13 @@ def unwrap(self, compound='group', center=None, box=None): "or 'fragments'".format(compound)) for o in objects: - make_whole(o) + make_whole(o.atoms) if center is None: tri_box = mdamath.triclinic_vectors(u.dimensions) center = np.diag(tri_box) / 2.0 - object_centers = np.vstack([o.center_of_mass(pbc=False) for o in objects]) + object_centers = np.vstack([o.center_of_geometry(pbc=False) for o in objects]) if box is None: box = self.dimensions diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index dec11c68ea7..d2aaadd0429 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 From 046abfde8b18b3452e9b652d0de8421767f4e051 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 28 Jul 2018 10:28:50 -0500 Subject: [PATCH 3/5] implemented make_whole inplace --- package/MDAnalysis/lib/_cutil.pyx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index d2aaadd0429..f15349a772c 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -277,7 +277,10 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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) @@ -317,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) From d30ebbd2d654475183369b845e151f340e1f9e7c Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sun, 29 Jul 2018 16:05:58 -0500 Subject: [PATCH 4/5] wip of unwrap in methods --- package/MDAnalysis/core/groups.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 8fc6b9edf95..1280868040e 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -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 From 408b084150434bd121d90a7df88f923862d09598 Mon Sep 17 00:00:00 2001 From: richardjgowers Date: Sat, 4 Aug 2018 14:15:15 -0500 Subject: [PATCH 5/5] more work on wrapping --- package/MDAnalysis/core/groups.py | 67 ++++++++++++++++++++++--------- 1 file changed, 48 insertions(+), 19 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 1280868040e..c75ba4aa536 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1141,7 +1141,7 @@ 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, box=None): + 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 @@ -1156,15 +1156,11 @@ def unwrap(self, compound='group', center=None, box=None): 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. - box : array_like - Box dimensions, can be either orthogonal or triclinic information. - Cell dimensions must be in an identical to format to those returned - by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`, - ``[lx, ly, lz, alpha, beta, gamma]``. If ``None``, uses these - timestep dimensions. # TODO the function isn't actually here, merge the _cutil in to make docs work .. seealso:: :func:`MDAnalysis.lib.util.make_whole` @@ -1186,25 +1182,56 @@ def unwrap(self, compound='group', center=None, box=None): "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) - if center is None: - tri_box = mdamath.triclinic_vectors(u.dimensions) - center = np.diag(tri_box) / 2.0 + def choose_image(self, reference_point, center='com', inplace=True): + """Shift coordinates to be as close as possible to reference point - object_centers = np.vstack([o.center_of_geometry(pbc=False) for o in objects]) + 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``] - if box is None: - box = self.dimensions + 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 - dests = distances.apply_PBC(object_centers, box=box) - shifts = dests - object_centers + .. versionadded:: 0.19.0 + """ + ag = self.atoms.unique - for o, s in zip(objects, shifts): - # Save some needless shifts - if not all(s == 0.0): - o.atoms.translate(s) + 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. @@ -1266,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])