From e8319b22af9f1dcf66276e84b18fac4f399c8678 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 30 Jan 2019 04:06:48 +0100 Subject: [PATCH 01/21] made inplace position modification of make_whole() optional --- package/MDAnalysis/lib/_cutil.pyx | 30 ++++++++++++++++++++++++------ package/MDAnalysis/lib/util.py | 2 +- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index 8a823d20c4e..194c6e74d77 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -103,7 +103,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. @@ -133,6 +133,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 ------ @@ -150,17 +157,17 @@ def make_whole(atomgroup, reference_atom=None): ------- Make fragments whole:: - from MDAnalysis.lib.mdamath import make_whole + from MDAnalysis.lib.util import make_whole # 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:: @@ -170,6 +177,9 @@ def make_whole(atomgroup, reference_atom=None): .. 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 @@ -187,6 +197,10 @@ def make_whole(atomgroup, reference_atom=None): # map of global indices to local indices ix_view = atomgroup.ix[:] natoms = atomgroup.ix.shape[0] + + if natoms == 0: + return atomgroup.positions + for i in range(natoms): ix_to_rel[ix_view[i]] = i @@ -198,6 +212,9 @@ def make_whole(atomgroup, reference_atom=None): raise ValueError("Reference atom not in atomgroup") ref = ix_to_rel[reference_atom.ix] + if natoms == 1: + return atomgroup.positions + box = atomgroup.dimensions for i in range(3): @@ -274,8 +291,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) diff --git a/package/MDAnalysis/lib/util.py b/package/MDAnalysis/lib/util.py index 28c9457198c..893448d1ad4 100644 --- a/package/MDAnalysis/lib/util.py +++ b/package/MDAnalysis/lib/util.py @@ -208,7 +208,7 @@ import inspect from ..exceptions import StreamWarning, DuplicateWarning -from ._cutil import unique_int_1d +from ._cutil import unique_int_1d, make_whole # Python 3.0, 3.1 do not have the builtin callable() From a43b6bea114b016bbf465f650db7d4bb9437b669 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 30 Jan 2019 04:12:14 +0100 Subject: [PATCH 02/21] added *Group.unwrap() method --- package/MDAnalysis/core/groups.py | 4 +- package/MDAnalysis/core/topologyattrs.py | 134 ++++++++++++++++++++++- package/MDAnalysis/lib/util.py | 2 +- 3 files changed, 136 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index a3f0e47120b..1415bf79387 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1842,11 +1842,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 diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index b4dd1ce3311..13693098c7b 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -50,8 +50,9 @@ from . import flags from ..lib.util import (cached, convert_aa_code, iterable, warn_if_not_unique, - unique_int_1d) + unique_int_1d, check_box) from ..lib import transformations, mdamath +from ..lib.distances import apply_PBC from ..exceptions import NoDataError, SelectionError from .topologyobjects import TopologyGroup from . import selection @@ -1799,6 +1800,137 @@ def n_fragments(self): ('n_fragments', property(n_fragments, None, None, n_fragments.__doc__))) + def unwrap(self, box=None, compound='fragments', reference='com', + inplace=True): + """Move atoms in this group so that bonds within the + group's compounds aren't split accross periodic boundaries. + + Parameters + ---------- + box : array_like or None, optional + The unitcell dimensions of the system, which can be orthogonal or + triclinic and must be provided in the same format as returned by + :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n + ``[lx, ly, lz, alpha, beta, gamma]``. + If ``None``, uses the dimensions of the current timestep. + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + Which type of components 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)``. + + Note + ---- + * This is the inverse transformation of packing atoms into the primary + unit cell, preventing bonds from being "broken" accross periodic + boundaries. + * 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 belong the group. + + See Also + -------- + :func:`MDAnalysis.lib.mdamath.make_whole` + + + .. versionadded:: 0.20.0 + """ + if reference is not None: + if reference.lower() not in ['com', 'cog']: + raise ValueError("Unrecognized reference definition: {}\n" + "Please use one of 'com', 'cog', or None." + "".format(compound)) + ref = reference.lower() + + if box is None: + box = self.dimensions + boxtype, box = check_box(box) + ortho = (boxtype == 'ortho') + if ortho: + hbox = 0.5 * box[:3] + + atoms = self.atoms + unique_atoms = atoms.unique + + comp = compound.lower() + if comp == 'group': + if ortho: + positions = unique_atoms.positions + # Only unwrap if necessary: + if np.any(np.abs(positions - positions[0]) > hbox): + positions = mdamath.make_whole(unique_atoms, inplace=False) + else: + positions = mdamath.make_whole(unique_atoms, inplace=False) + # Apply reference shift if required: + if reference is not None: + if ref == 'com': + refpos = unique_atoms.center_of_mass(pbc=False, + compound='group') + elif ref == 'cog': + refpos = unique_atoms.center_of_geometry(pbc=False, + compound='group') + target = apply_PBC(refpos, self.dimensions) + positions += target - refpos + if inplace: + unique_atoms.positions = positions + if not atoms.isunique: + positions = positions[atoms._unique_restore_mask] + return positions + elif comp == 'fragments': + compound_indices = unique_atoms.fragindices + elif comp == 'residues': + compound_indices = unique_atoms.resindices + elif comp == 'segments': + compound_indices = unique_atoms.segindices + elif comp == 'molecules': + try: + compound_indices = unique_atoms.molnums + except AttributeError: + raise NoDataError("Cannot use compound='molecules': " + "No molecule information in topology.") + else: + raise ValueError("Unrecognized compound definition: {}\nPlease use" + " one of 'group', 'residues', 'segments', " + "'molecules', or 'fragments'.".format(compound)) + + # 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] + cpos = positions[mask] + # Only unwrap if necessary: + if (not ortho) or np.any(np.abs(cpos - cpos[0]) > hbox): + positions[mask] = mdamath.make_whole(c, inplace=False) + # Apply reference shift if required: + if reference is not None: + if ref == 'com': + refpos = c.center_of_mass(pbc=False, compound='group') + elif ref == 'cog': + refpos = c.center_of_geometry(pbc=False, compound='group') + target = 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 + + transplants[GroupBase].append( + ('unwrap', unwrap)) + class Angles(_Connection): """Angles between three atoms diff --git a/package/MDAnalysis/lib/util.py b/package/MDAnalysis/lib/util.py index 893448d1ad4..28c9457198c 100644 --- a/package/MDAnalysis/lib/util.py +++ b/package/MDAnalysis/lib/util.py @@ -208,7 +208,7 @@ import inspect from ..exceptions import StreamWarning, DuplicateWarning -from ._cutil import unique_int_1d, make_whole +from ._cutil import unique_int_1d # Python 3.0, 3.1 do not have the builtin callable() From 44cd7097ee38c4acab06373a365d79a85d56acbd Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 30 Jan 2019 04:47:51 +0100 Subject: [PATCH 03/21] fixed docs of make_whole() --- package/MDAnalysis/lib/_cutil.pyx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index 194c6e74d77..92cc684f548 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -104,9 +104,8 @@ cdef intset difference(intset a, intset b): @cython.boundscheck(False) @cython.wraparound(False) 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. + """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 From e89a89d6f71c7117e79b2002046c41ad29778205 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 30 Jan 2019 19:54:56 +0100 Subject: [PATCH 04/21] cleaned up unwrap() --- package/MDAnalysis/core/groups.py | 4 +- package/MDAnalysis/core/topologyattrs.py | 167 +++++++++++++---------- package/MDAnalysis/lib/_cutil.pyx | 27 +++- 3 files changed, 117 insertions(+), 81 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 1415bf79387..7477e5ea1cc 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -654,7 +654,9 @@ def center(self, weights, pbc=None, compound='group'): Computes the weighted center of :class:`Atoms` 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 ---------- diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 13693098c7b..66b513db5e0 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -50,7 +50,7 @@ from . import flags from ..lib.util import (cached, convert_aa_code, iterable, warn_if_not_unique, - unique_int_1d, check_box) + unique_int_1d) from ..lib import transformations, mdamath from ..lib.distances import apply_PBC from ..exceptions import NoDataError, SelectionError @@ -776,7 +776,9 @@ def center_of_mass(group, pbc=None, compound='group'): Computes the center of mass of :class:`Atoms` 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 ---------- @@ -1800,20 +1802,29 @@ def n_fragments(self): ('n_fragments', property(n_fragments, None, None, n_fragments.__doc__))) - def unwrap(self, box=None, compound='fragments', reference='com', - inplace=True): - """Move atoms in this group so that bonds within the + 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 accross 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 ---------- - box : array_like or None, optional - The unitcell dimensions of the system, which can be orthogonal or - triclinic and must be provided in the same format as returned by - :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:\n - ``[lx, ly, lz, alpha, beta, gamma]``. - If ``None``, uses the dimensions of the current timestep. - compound : {'group', 'segments', 'residues', 'molecules', 'fragments'}, optional + compound : {'group', 'segments', 'residues', 'molecules', \ + 'fragments'}, optional Which type of components 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. @@ -1830,51 +1841,66 @@ def unwrap(self, box=None, compound='fragments', reference='com', 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 ---- - * This is the inverse transformation of packing atoms into the primary - unit cell, preventing bonds from being "broken" accross periodic - boundaries. * 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 belong the group. + 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` + :func:`~MDAnalysis.lib.mdamath.make_whole`, + :meth:`~MDAnalysis.AtomGroup.wrap`, + :meth:`~MDAnalysis.AtomGroup.pack_into_box`, + :func:`~MDanalysis.lib.distances.apply_PBC` .. versionadded:: 0.20.0 """ if reference is not None: - if reference.lower() not in ['com', 'cog']: - raise ValueError("Unrecognized reference definition: {}\n" - "Please use one of 'com', 'cog', or None." - "".format(compound)) ref = reference.lower() - - if box is None: - box = self.dimensions - boxtype, box = check_box(box) - ortho = (boxtype == 'ortho') - if ortho: - hbox = 0.5 * box[:3] + 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)) atoms = self.atoms unique_atoms = atoms.unique 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': - if ortho: - positions = unique_atoms.positions - # Only unwrap if necessary: - if np.any(np.abs(positions - positions[0]) > hbox): - positions = mdamath.make_whole(unique_atoms, inplace=False) - else: - positions = mdamath.make_whole(unique_atoms, inplace=False) + positions = mdamath.make_whole(unique_atoms, inplace=False) # Apply reference shift if required: if reference is not None: if ref == 'com': + if np.isclose(unique_atoms.total_mass(), 0.0): + raise ValueError("Cannot perform unwrap with " + "reference='com' because the total " + "mass of the group is zero.") refpos = unique_atoms.center_of_mass(pbc=False, compound='group') elif ref == 'cog': @@ -1882,46 +1908,41 @@ def unwrap(self, box=None, compound='fragments', reference='com', compound='group') target = apply_PBC(refpos, self.dimensions) positions += target - refpos - if inplace: - unique_atoms.positions = positions - if not atoms.isunique: - positions = positions[atoms._unique_restore_mask] - return positions - elif comp == 'fragments': - compound_indices = unique_atoms.fragindices - elif comp == 'residues': - compound_indices = unique_atoms.resindices - elif comp == 'segments': - compound_indices = unique_atoms.segindices - elif comp == 'molecules': - try: - compound_indices = unique_atoms.molnums - except AttributeError: - raise NoDataError("Cannot use compound='molecules': " - "No molecule information in topology.") + # We need to split the group into compounds: else: - raise ValueError("Unrecognized compound definition: {}\nPlease use" - " one of 'group', 'residues', 'segments', " - "'molecules', or 'fragments'.".format(compound)) - - # 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] - cpos = positions[mask] - # Only unwrap if necessary: - if (not ortho) or np.any(np.abs(cpos - cpos[0]) > hbox): + if comp == 'fragments': + compound_indices = unique_atoms.fragindices + elif comp == 'residues': + compound_indices = unique_atoms.resindices + elif comp == 'segments': + compound_indices = unique_atoms.segindices + elif 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': - refpos = c.center_of_mass(pbc=False, compound='group') - elif ref == 'cog': - refpos = c.center_of_geometry(pbc=False, compound='group') - target = apply_PBC(refpos, self.dimensions) - positions[mask] += target - refpos + # Apply reference shift if required: + if reference is not None: + if ref == 'com': + if np.isclose(c.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 = c.center_of_mass(pbc=False, compound='group') + elif ref == 'cog': + refpos = c.center_of_geometry(pbc=False, + compound='group') + target = apply_PBC(refpos, self.dimensions) + positions[mask] += target - refpos if inplace: unique_atoms.positions = positions if not atoms.isunique: diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index 92cc684f548..8e0e5781eb1 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -181,7 +181,7 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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 @@ -189,16 +189,21 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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 - if natoms == 0: - return 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 @@ -211,12 +216,10 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): raise ValueError("Reference atom not in atomgroup") ref = ix_to_rel[reference_atom.ix] - if natoms == 1: - return atomgroup.positions - box = atomgroup.dimensions 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='") @@ -227,6 +230,17 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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 abs(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: @@ -249,7 +263,6 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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? From fc7cc384f65fec25781be5c309d203aa04bc1aaf Mon Sep 17 00:00:00 2001 From: zeman Date: Fri, 1 Feb 2019 16:06:20 +0100 Subject: [PATCH 05/21] doc fixes --- package/MDAnalysis/core/topologyattrs.py | 35 ++++++++++++------------ 1 file changed, 18 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 66b513db5e0..edc9e48e20d 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1807,25 +1807,25 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): group's compounds aren't split accross 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 + 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 - | | | | - +-----------+ +-----------+ + +-----------+ +-----------+ + | | | | + | 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 components to unwrap. Note that, in any case, all + 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 @@ -1854,13 +1854,14 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): 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 + 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) - >>> all_frag_atoms = sum(ag.fragments) See Also -------- From e04ec6e4af2382fd7b322396b7fec3d67b06cadd Mon Sep 17 00:00:00 2001 From: zeman Date: Fri, 1 Feb 2019 16:06:59 +0100 Subject: [PATCH 06/21] fix make_whole() for double precision boxes --- package/MDAnalysis/lib/_cutil.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index 8e0e5781eb1..e100c160fc0 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -216,7 +216,8 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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] From 737b4feff892fddb5cf0ce7547b8d019b07289a1 Mon Sep 17 00:00:00 2001 From: zeman Date: Fri, 1 Feb 2019 17:44:37 +0100 Subject: [PATCH 07/21] fixed group.unwrap() for reference="com" --- package/MDAnalysis/core/topologyattrs.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index edc9e48e20d..4a3fde8f513 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1898,15 +1898,16 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): # Apply reference shift if required: if reference is not None: if ref == 'com': - if np.isclose(unique_atoms.total_mass(), 0.0): + 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 = unique_atoms.center_of_mass(pbc=False, - compound='group') + refpos = np.sum(positions * masses[:, None], axis=0) + refpos /= total_mass elif ref == 'cog': - refpos = unique_atoms.center_of_geometry(pbc=False, - compound='group') + refpos = positions.mean(axis=0) target = apply_PBC(refpos, self.dimensions) positions += target - refpos # We need to split the group into compounds: @@ -1933,15 +1934,18 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): # Apply reference shift if required: if reference is not None: if ref == 'com': - if np.isclose(c.total_mass(), 0.0): + 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 = c.center_of_mass(pbc=False, compound='group') + refpos = np.sum(positions[mask] * masses[:, None], + axis=0) + refpos /= total_mass elif ref == 'cog': - refpos = c.center_of_geometry(pbc=False, - compound='group') + refpos = positions[mask].mean(axis=0) target = apply_PBC(refpos, self.dimensions) positions[mask] += target - refpos if inplace: From 781b0eefcbd2bea0f13f4eae07b11088b36e2322 Mon Sep 17 00:00:00 2001 From: zeman Date: Fri, 1 Feb 2019 19:28:15 +0100 Subject: [PATCH 08/21] fix group.unwrap() for empty group --- package/MDAnalysis/core/topologyattrs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 4a3fde8f513..11dad2e86ab 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1896,7 +1896,7 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): if comp == 'group': positions = mdamath.make_whole(unique_atoms, inplace=False) # Apply reference shift if required: - if reference is not None: + if reference is not None and len(positions) > 0: if ref == 'com': masses = unique_atoms.masses total_mass = masses.sum() From 19c4190908cc632601dd0a980c1b00491abfffb5 Mon Sep 17 00:00:00 2001 From: zeman Date: Mon, 4 Feb 2019 14:20:25 +0100 Subject: [PATCH 09/21] moved *group.unwrap() from core.topologyattrs to core.groups --- package/MDAnalysis/core/groups.py | 155 +++++++++++++++++++++++ package/MDAnalysis/core/topologyattrs.py | 155 ----------------------- 2 files changed, 155 insertions(+), 155 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 7477e5ea1cc..6f36fdb7c84 100644 --- a/package/MDAnalysis/core/groups.py +++ b/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 @@ -1255,6 +1256,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 accross 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:`~MDAnalysis.AtomGroup.wrap`, + :meth:`~MDAnalysis.AtomGroup.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 + elif 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 + elif 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 + elif 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. diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 11dad2e86ab..2f7032344d3 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1802,161 +1802,6 @@ def n_fragments(self): ('n_fragments', property(n_fragments, None, None, n_fragments.__doc__))) - 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 accross 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:`~MDAnalysis.AtomGroup.wrap`, - :meth:`~MDAnalysis.AtomGroup.pack_into_box`, - :func:`~MDanalysis.lib.distances.apply_PBC` - - - .. versionadded:: 0.20.0 - """ - 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)) - - atoms = self.atoms - unique_atoms = atoms.unique - - 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 - elif ref == 'cog': - refpos = positions.mean(axis=0) - target = 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 - elif 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 - elif ref == 'cog': - refpos = positions[mask].mean(axis=0) - target = 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 - - transplants[GroupBase].append( - ('unwrap', unwrap)) - class Angles(_Connection): """Angles between three atoms From 5bd767f3fd02be69c6e7913a05414de8d2988faf Mon Sep 17 00:00:00 2001 From: zeman Date: Mon, 4 Feb 2019 14:21:05 +0100 Subject: [PATCH 10/21] tests for *group.unwrap() --- testsuite/MDAnalysisTests/core/test_unwrap.py | 772 ++++++++++++++++++ 1 file changed, 772 insertions(+) create mode 100644 testsuite/MDAnalysisTests/core/test_unwrap.py diff --git a/testsuite/MDAnalysisTests/core/test_unwrap.py b/testsuite/MDAnalysisTests/core/test_unwrap.py new file mode 100644 index 00000000000..f53716c0450 --- /dev/null +++ b/testsuite/MDAnalysisTests/core/test_unwrap.py @@ -0,0 +1,772 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +from __future__ import absolute_import, division +from six.moves import range + +import numpy as np +from numpy.testing import ( + assert_equal, + assert_almost_equal, + assert_array_equal +) +import pytest + +import MDAnalysis as mda +from MDAnalysis.core import topologyattrs +from MDAnalysis import NoDataError +from MDAnalysis.lib.mdamath import triclinic_box + +from MDAnalysisTests import make_Universe +from MDAnalysisTests.datafiles import TPR, XTC + + +def unwrap_test_universe(have_bonds=True, have_masses=True, have_molnums=True, + is_triclinic=False): + """Construct a universe comprising small molecule clusters with molecules + (as well as clusters) broken or whole accross periodic boundaries: + + The universe comprises 47 atoms in 12 molecules: + * Three "molecules" (type A) containing a single atom: + 1. within the central image + 2. in a neighboring image + 3. in an image two boxes away from the central image + + * Four molecules (type B) containing three atoms: + 4. within the central image + 5. in the corner of the box, close to 4., broken accross one boundary + but whole accross another + 6. outside the box, broken diagonally accross two neighboring images, + close to 5. + 7. whole accross a boundary opposite to the location of 4., i.e., + close to 4. in terms of PBC + * Two cyclical molecules (type C) containing four atoms: + 8. broken accross the front/back but whole accross the top face + 9. within the central image close to the front and bottom face, close + to 8. in terms of PBC + * Three linear chain molecules (type D) containing 8 atoms, spanning + more than half a box length: + 10. within the central image relatively close to the top boundary + 11. close to 10, broken mid-molecule accross the left/right boundary + 12. close to 11. but in another image, whole accross the same + boundary as 11. but not mid-molecule + + There are 15 residues in the universe: + Molecules of type A, B, and C each have a single residue, while each of + the chain molecules of type D have two residues with 4 atoms per residue. + + Atoms can be selected by their residue's resname. For molecules of type + A, B, and C, the resnames are A, B, and C, respectively. Molecules of + type D have the resnames D1 and D2. + + Atoms can also be selected by their moltype attribute, which is identical + to the corresponding molecule type (A, B, C, D). + + The molecules/residues are contained in 6 segments, whereof the first + three segments contain all molecules of type A, B, or C, respectively. + Each of the remaining three segments contain a molecule of type D. + + A projection onto the x-z plane of box (orthorhombic case) looks roughly + like this:: + + : : + : : + 6 : 8 : : + - + - -:+-------+----------:- - - - - - - - - -: - + |5 8 | : + | | : + | (10) | : + | o-o-o-o-x-x-x-x | : + +x-x-x-x o-o-o-o+ : + | (11) | : + | | : + | 1 | 2 : 3 + | | : + | 9 | : + | 4 ! 7| : + | ! 9 !| : + | 4-4 7+7 : + | | : + 5+5 | : + - + - -:+------------------:- - - - - - - - - -: - + 6-6 : : : + : : + : : + : : + : (12) : + : x-x-x-x-o-o+o-o + : : + + Note that the cyclic structures of molecules 8 and 9 lie in the x-y plane + and are therefore not fully visible in the x-z plane projection. + + Parameters + ---------- + have_bonds : bool, optional + If ``True``, intramolecular bonds will be present in the topology. + If ``False``, there will be no bond information. + have_masses : bool, optional + If ``True``, each atom will be assigned a mass of 1 u. + If ``False``, masses won't be present in the topology. + have_molnums : bool, optional + If ``True``, the topology will have molecule information (molnums). + If ``False``, that information will be missing. + is_triclinic : bool, optional + If ``False``, the box will be a cube with an edge length of 10 Angstrom. + If ``True``, the cubic box will be sheared by 1 Angstrom in the x and + y directions. + + Returns + ------- + u : MDAnalysis.Universe + The test universe as specified above. + """ + # box: + a = 10.0 # edge length + tfac = 0.1 # factor for box vector shift of triclinic boxes (~84°) + if is_triclinic: + box = triclinic_box([a, 0.0, 0.0], + [a * tfac, a, 0.0], + [a * tfac, a * tfac, a]) + else: + box = np.array([a, a, a, 90.0, 90.0, 90.0], dtype=np.float32) + + # number of atoms, residues, and segments: + n_atoms = 47 + n_residues = 15 + n_segments = 6 + + # resindices: + residx = np.empty(n_atoms, dtype=np.int64) + # type A + rix = 0 + for i in range(0, 3, 1): + residx[i] = rix + rix += 1 + # type B + for i in range(3, 15, 3): + residx[i:i+3] = rix + rix += 1 + # type C & D + for i in range(15, 47, 4): + residx[i:i+4] = rix + rix += 1 + + # segindices: + segidx = np.empty(n_residues, dtype=np.int64) + segidx[0:3] = 0 + segidx[3:7] = 1 + segidx[7:9] = 2 + segidx[9:11] = 3 + segidx[11:13] = 4 + segidx[13:15] = 5 + + # universe: + u = mda.Universe.empty( + # topology things + n_atoms=n_atoms, + n_residues=n_residues, + n_segments=n_segments, + atom_resindex=residx, + residue_segindex=segidx, + # trajectory things + trajectory=True, + velocities=False, + forces=False, + ) + + # resnames: we always want those for selection purposes + resnames = ['A'] * 3 + resnames += ['B'] * 4 + resnames += ['C'] * 2 + resnames += ['D1', 'D2'] * 3 + u.add_TopologyAttr(topologyattrs.Resnames(resnames)) + + # moltypes: we always want those for selection purposes + moltypes = ['A'] * 3 + moltypes += ['B'] * 4 + moltypes += ['C'] * 2 + moltypes += ['D'] * 6 + u.add_TopologyAttr(topologyattrs.Moltypes(moltypes)) + + # trajectory: + ts = u.trajectory.ts + ts.frame = 0 + ts.dimensions = box + + # positions: + relpos = np.empty((n_atoms, 3), dtype=np.float32) + # type A + relpos[0:3, :] = np.array([[0.5, 0.5, 0.5], + [1.4, 0.5, 0.5], + [2.1, 0.5, 0.5]], dtype=np.float32) + # type B + relpos[3:15, :] = np.array([[0.1, 0.1, 0.2], + [0.1, 0.1, 0.1], + [0.2, 0.1, 0.1], + [-0.05, 0.2, 0.05], + [0.05, 0.2, 0.05], + [0.05, 0.2, 0.95], + [-0.2, -0.9, 1.05], + [-0.2, 0.1, -0.05], + [-0.1, 0.1, -0.05], + [0.95, 0.2, 0.25], + [0.95, 0.2, 0.15], + [1.05, 0.2, 0.15]], dtype=np.float32) + # type C + relpos[15:23, :] = np.array([[0.4, 0.96, 1.06], + [0.4, 0.96, 0.96], + [0.4, 0.06, 0.96], + [0.4, 0.06, 1.06], + [0.6, 0.06, 0.26], + [0.6, 0.06, 0.16], + [0.6, 0.16, 0.16], + [0.6, 0.16, 0.26]], dtype=np.float32) + # type D + relpos[23:47, :] = np.array([[0.2, 0.7, 0.8], + [0.3, 0.7, 0.8], + [0.4, 0.7, 0.8], + [0.5, 0.7, 0.8], + [0.6, 0.7, 0.8], + [0.7, 0.7, 0.8], + [0.8, 0.7, 0.8], + [0.9, 0.7, 0.8], + [0.66, 0.75, 0.7], + [0.76, 0.75, 0.7], + [0.86, 0.75, 0.7], + [0.96, 0.75, 0.7], + [0.06, 0.75, 0.7], + [0.16, 0.75, 0.7], + [0.26, 0.75, 0.7], + [0.36, 0.75, 0.7], + [1.14, 0.65, -0.4], + [1.04, 0.65, -0.4], + [0.94, 0.65, -0.4], + [0.84, 0.65, -0.4], + [0.74, 0.65, -0.4], + [0.64, 0.65, -0.4], + [0.54, 0.65, -0.4], + [0.44, 0.65, -0.4]], dtype=np.float32) + # apply y- and z-dependent shift of x and y coords for triclinic boxes: + if is_triclinic: + # x-coord shift depends on y- and z-coords + relpos[:, 0] += tfac * relpos[:, 1:].sum(axis=1) + # y-coord shift depends on z-coords only + relpos[:, 1] += tfac * relpos[:, 2] + # scale relative to absolute positions: + ts.positions = relpos * np.array([a, a, a], dtype=np.float32) + + # bonds: + if have_bonds: + bonds = [] + # type A has no bonds + #type B + for base in range(3, 15, 3): + for i in range(2): + bonds.append((base + i, base + i + 1)) + # type C + for base in range(15, 23, 4): + for i in range(3): + bonds.append((base + i, base + i + 1)) + bonds.append((0 + base, 3 + base)) + # type D + for base in range(23, 47, 8): + for i in range(7): + bonds.append((base + i, base + i + 1)) + u.add_TopologyAttr(topologyattrs.Bonds(bonds)) + + # masses: + if have_masses: + # masses are all set to 1 so that one can cross-check the results of + # reference='com' with reference='cog' unwrapping + masses = np.ones(n_atoms) + u.add_TopologyAttr(topologyattrs.Masses(masses)) + + # molnums: + if have_molnums: + molnums = [0, 1, 2] + molnums += [3, 4, 5, 6] + molnums += [7, 8] + molnums += [9, 9, 10, 10, 11, 11] + u.add_TopologyAttr(topologyattrs.Molnums(molnums)) + + return u + + +def unwrapped_coords(compound, reference, is_triclinic): + """Returns coordinates which correspond to the system returned + by :func:`unwrap_test_universe` but are unwrapped. + + Parameters + ---------- + compound : {'group', 'segments', 'residues', 'molecules', 'fragments'} + Which type of component is unwrapped. + reference : {'com', 'cog', None} + The reference point of each compound that is shifted into the primary + unit cell. + is_triclinic : bool + Is the box triclinic? + + Note + ---- + This function assumes that all atom masses are equal. Therefore, the + returned coordinates for ``reference='com'`` and ``reference='cog'`` are + identical. + """ + if reference not in ['com', 'cog', None]: + raise ValueError("Unknown unwrap reference: {}".format(reference)) + if compound not in ['group', 'segments', 'residues', 'molecules', + 'fragments']: + raise ValueError("Unknown unwrap compound: {}".format(compound)) + + n_atoms = 47 + + # box: + a = 10.0 # edge length + tfac = 0.1 # factor for box vector shift of triclinic boxes (~84°) + if is_triclinic: + box = triclinic_box([a, 0.0, 0.0], + [a * tfac, a, 0.0], + [a * tfac, a * tfac, a]) + else: + box = np.array([a, a, a, 90.0, 90.0, 90.0], dtype=np.float32) + + # unwrapped positions: + relpos = np.empty((n_atoms, 3), dtype=np.float32) + # type A + relpos[0:3, :] = np.array([[0.5, 0.5, 0.5], + [1.4, 0.5, 0.5], + [2.1, 0.5, 0.5]], dtype=np.float32) + # type B + relpos[3:15, :] = np.array([[0.1, 0.1, 0.2], + [0.1, 0.1, 0.1], + [0.2, 0.1, 0.1], + [-0.05, 0.2, 0.05], + [0.05, 0.2, 0.05], + [0.05, 0.2, -0.05], + [-0.2, -0.9, 1.05], + [-0.2, -0.9, 0.95], + [-0.1, -0.9, 0.95], + [0.95, 0.2, 0.25], + [0.95, 0.2, 0.15], + [1.05, 0.2, 0.15]], dtype=np.float32) + # type C + relpos[15:23, :] = np.array([[0.4, 0.96, 1.06], + [0.4, 0.96, 0.96], + [0.4, 1.06, 0.96], + [0.4, 1.06, 1.06], + [0.6, 0.06, 0.26], + [0.6, 0.06, 0.16], + [0.6, 0.16, 0.16], + [0.6, 0.16, 0.26]], dtype=np.float32) + # type D + relpos[23:47, :] = np.array([[0.2, 0.7, 0.8], + [0.3, 0.7, 0.8], + [0.4, 0.7, 0.8], + [0.5, 0.7, 0.8], + [0.6, 0.7, 0.8], + [0.7, 0.7, 0.8], + [0.8, 0.7, 0.8], + [0.9, 0.7, 0.8], + [0.66, 0.75, 0.7], + [0.76, 0.75, 0.7], + [0.86, 0.75, 0.7], + [0.96, 0.75, 0.7], + [1.06, 0.75, 0.7], + [1.16, 0.75, 0.7], + [1.26, 0.75, 0.7], + [1.36, 0.75, 0.7], + [1.14, 0.65, -0.4], + [1.04, 0.65, -0.4], + [0.94, 0.65, -0.4], + [0.84, 0.65, -0.4], + [0.74, 0.65, -0.4], + [0.64, 0.65, -0.4], + [0.54, 0.65, -0.4], + [0.44, 0.65, -0.4]], dtype=np.float32) + # apply image shifts if necessary: + if reference is None: + if compound == 'residues': + #second residue of molecule 11 + relpos[35:39, :] = np.array([[0.06, 0.75, 0.7], + [0.16, 0.75, 0.7], + [0.26, 0.75, 0.7], + [0.36, 0.75, 0.7]], dtype=np.float32) + else: + # molecule 2 & 3 + relpos[1:3, :] = np.array([[0.4, 0.5, 0.5], + [0.1, 0.5, 0.5]], dtype=np.float32) + # molecule 6 + relpos[9:12, :] = np.array([[0.8, 0.1, 1.05], + [0.8, 0.1, 0.95], + [0.9, 0.1, 0.95]], dtype=np.float32) + #molecule 8 + relpos[15:19, :] = np.array([[0.4, -0.04, 0.06], + [0.4, -0.04, -0.04], + [0.4, 0.06, -0.04], + [0.4, 0.06, 0.06]], dtype=np.float32) + if compound == 'residues': + #molecule 11 & 12 + relpos[31:47, :] = np.array([[0.66, 0.75, 0.7], + [0.76, 0.75, 0.7], + [0.86, 0.75, 0.7], + [0.96, 0.75, 0.7], + [0.06, 0.75, 0.7], + [0.16, 0.75, 0.7], + [0.26, 0.75, 0.7], + [0.36, 0.75, 0.7], + [1.14, 0.65, 0.6], + [1.04, 0.65, 0.6], + [0.94, 0.65, 0.6], + [0.84, 0.65, 0.6], + [0.74, 0.65, 0.6], + [0.64, 0.65, 0.6], + [0.54, 0.65, 0.6], + [0.44, 0.65, 0.6]], dtype=np.float32) + else: + #molecule 11 & 12 + relpos[31:47, :] = np.array([[-0.34, 0.75, 0.7], + [-0.24, 0.75, 0.7], + [-0.14, 0.75, 0.7], + [-0.04, 0.75, 0.7], + [0.06, 0.75, 0.7], + [0.16, 0.75, 0.7], + [0.26, 0.75, 0.7], + [0.36, 0.75, 0.7], + [1.14, 0.65, 0.6], + [1.04, 0.65, 0.6], + [0.94, 0.65, 0.6], + [0.84, 0.65, 0.6], + [0.74, 0.65, 0.6], + [0.64, 0.65, 0.6], + [0.54, 0.65, 0.6], + [0.44, 0.65, 0.6]], dtype=np.float32) + + # apply y- and z-dependent shift of x and y coords for triclinic boxes: + if is_triclinic: + # x-coord shift depends on y- and z-coords + relpos[:, 0] += tfac * relpos[:, 1:].sum(axis=1) + # y-coord shift depends on z-coords only + relpos[:, 1] += tfac * relpos[:, 2] + + # scale relative to absolute positions: + positions = relpos * np.array([a, a, a], dtype=np.float32) + + return positions + + +class TestUnwrap(object): + """Tests the functionality of *Group.unwrap() using the test universe + returned by unwrap_test_universe() + """ + precision = 5 + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_unwrap_pass(self, level, compound, reference, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + # get the expected result: + ref_unwrapped_pos = unwrapped_coords(compound, reference, + is_triclinic=is_triclinic) + if compound == 'group': + ref_unwrapped_pos = ref_unwrapped_pos[39:47] # molecule 12 + elif compound == 'segments': + ref_unwrapped_pos = ref_unwrapped_pos[23:47] # molecules 10, 11, 12 + # first, do the unwrapping out-of-place: + unwrapped_pos = group.unwrap(compound=compound, reference=reference, + inplace=False) + # check for correct result: + assert_almost_equal(unwrapped_pos, ref_unwrapped_pos, + decimal=self.precision) + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + # now, do the unwrapping inplace: + unwrapped_pos2 = group.unwrap(compound=compound, reference=reference, + inplace=True) + # check that result is the same as for out-of-place computation: + assert_array_equal(unwrapped_pos, unwrapped_pos2) + # check that unwrapped positions are applied: + assert_array_equal(group.atoms.positions, unwrapped_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_wrap_unwrap_cycle(self, level, compound, reference, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # wrap: + group.wrap() + # store original wrapped positions: + orig_wrapped_pos = group.atoms.positions + # unwrap: + group.unwrap(compound=compound, reference=reference, inplace=True) + # wrap again: + group.wrap() + # make sure wrapped atom positions are as before: + assert_almost_equal(group.atoms.positions, orig_wrapped_pos, + decimal=self.precision) + + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_unwrap_partial_frags(self, compound, reference, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + # select group with one atom missing + group = u.atoms[39:46] # molecule 12 without its last atom + # store original position of last atom of molecule 12: + orig_pos = u.atoms[46].position + # get the expected result: + ref_unwrapped_pos = unwrapped_coords(compound, reference, + is_triclinic=is_triclinic)[39:46] + # first, do the unwrapping out-of-place: + group.unwrap(compound=compound, reference=reference, inplace=True) + # check for correct result: + assert_almost_equal(group.positions, ref_unwrapped_pos, + decimal=self.precision) + # make sure the position of molecule 12's last atom is unchanged: + assert_array_equal(u.atoms[46].position, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_unwrap_empty_group(self, level, compound, reference, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + if level == 'atoms': + group = mda.AtomGroup([], u) + elif level == 'residues': + group = mda.ResidueGroup([], u) + elif level == 'segments': + group = mda.SegmentGroup([], u) + group.unwrap(compound=compound, reference=reference, inplace=True) + # check for correct (empty) result: + assert_array_equal(group.atoms.positions, + np.empty((0, 3), dtype=np.float32)) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_unwrap_duplicates(self, level, compound, reference, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + # select molecule 6: + group = u.atoms[39:47] + # select the rest of the universe's atoms: + rest = u.atoms[:39] + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # duplicate the group: + group += group + # store original positions of the rest: + orig_rest_pos = rest.positions + # get the expected result with duplicates: + ref_unwrapped_pos = unwrapped_coords(compound, reference, + is_triclinic=is_triclinic)[39:47] + ref_unwrapped_pos = np.vstack((ref_unwrapped_pos, ref_unwrapped_pos)) + # unwrap: + group.unwrap(compound=compound, reference=reference, inplace=True) + # check for correct result: + assert_almost_equal(group.atoms.positions, ref_unwrapped_pos, + decimal=self.precision) + # check that the rest of the atoms are kept unmodified: + assert_array_equal(rest.positions, orig_rest_pos) + + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('is_triclinic', (False, True)) + def test_unwrap_com_cog_difference(self, compound, is_triclinic): + # get a pristine test universe: + u = unwrap_test_universe(is_triclinic=is_triclinic) + # select molecule 5: + group = u.atoms[6:9] + # make first atom of molecule 5 much more heavy than the other two. + # That way, the whole molecule's center of geometry will still lie + # inside the first unit cell but its center of mass will lie outside + # the first unit cell in negative x-direction. + group.masses = [100.0, 1.0, 1.0] + # unwrap with center of geometry as reference: + unwrapped_pos_cog = group.unwrap(compound=compound, reference='cog', + inplace=False) + # get expected result: + ref_unwrapped_pos = unwrapped_coords(compound, reference='cog', + is_triclinic=is_triclinic)[6:9] + # check for correctness: + assert_almost_equal(unwrapped_pos_cog, ref_unwrapped_pos, + decimal=self.precision) + # unwrap with center of mass as reference: + unwrapped_pos_com = group.unwrap(compound=compound, reference='com', + inplace=False) + # assert that the com result is shifted with respect to the cog result + # by one box length in the x-direction: + shift = np.array([10.0, 0.0, 0.0], dtype=np.float32) + assert_almost_equal(unwrapped_pos_cog, unwrapped_pos_com - shift, + decimal=self.precision) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + def test_unwrap_zero_mass_exception_safety(self, level, compound): + # get a pristine test universe: + u = unwrap_test_universe() + # set masses of molecule 12 to zero: + u.atoms[39:47].masses = 0.0 + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + # try to unwrap: + with pytest.raises(ValueError): + group.unwrap(compound=compound, reference='com', + inplace=True) + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + def test_unwrap_wrong_reference_exception_safety(self, level, compound): + # get a pristine test universe: + u = unwrap_test_universe() + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + # try to unwrap: + with pytest.raises(ValueError): + group.unwrap(compound=compound, reference='wrong', inplace=True) + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + def test_unwrap_wrong_compound_exception_safety(self, level, reference): + # get a pristine test universe: + u = unwrap_test_universe() + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + # try to unwrap: + with pytest.raises(ValueError): + group.unwrap(compound='wrong', reference=reference, inplace=True) + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + def test_unwrap_no_bonds_exception_safety(self, level): + # universe without bonds: + u = unwrap_test_universe(have_bonds=False) + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + with pytest.raises(NoDataError): + group.unwrap() + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + def test_unwrap_no_molnums_exception_safety(self, level, reference): + # universe without molnums: + u = unwrap_test_universe(have_molnums=False) + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + with pytest.raises(NoDataError): + group.unwrap(compound='molecules', reference=reference, + inplace=True) + assert_array_equal(group.atoms.positions, orig_pos) From b7a0760dc2c96b003f7bb20fd6957133306fb5bd Mon Sep 17 00:00:00 2001 From: zeman Date: Mon, 4 Feb 2019 16:16:59 +0100 Subject: [PATCH 11/21] tests for modifications in make_whole() --- testsuite/MDAnalysisTests/lib/test_util.py | 23 +++++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/lib/test_util.py b/testsuite/MDAnalysisTests/lib/test_util.py index 197f68f0de6..5764422f88b 100644 --- a/testsuite/MDAnalysisTests/lib/test_util.py +++ b/testsuite/MDAnalysisTests/lib/test_util.py @@ -254,6 +254,14 @@ def universe(self): universe.add_TopologyAttr(Bonds(bondlist)) return universe + def test_return_value(self, universe): + ag = universe.residues[0].atoms + orig_pos = ag.positions.copy() + retval = mdamath.make_whole(ag) + assert retval.dtype == np.float32 + assert_array_equal(ag.positions, retval) + assert np.any(ag.positions != orig_pos) + def test_single_atom_no_bonds(self): # Call make_whole on single atom with no bonds, shouldn't move u = mda.Universe(Make_Whole) @@ -265,7 +273,13 @@ def test_single_atom_no_bonds(self): refpos = ag.positions.copy() mdamath.make_whole(ag) - assert_array_almost_equal(ag.positions, refpos) + assert_array_equal(ag.positions, refpos) # must be untouched + + def test_empty_ag(self, universe): + ag = mda.AtomGroup([], universe) + retval = mdamath.make_whole(ag) + assert retval.dtype == np.float32 + assert_array_equal(retval, np.empty((0, 3), dtype=np.float32)) def test_scrambled_ag(self, universe): # if order of atomgroup is mixed @@ -277,6 +291,13 @@ def test_scrambled_ag(self, universe): # largest bond should be 20A assert ag.bonds.values().max() < 20.1 + def test_out_of_place(self, universe): + ag = universe.residues[0].atoms + orig_pos = ag.positions.copy() + mdamath.make_whole(ag, inplace=False) + # positions must be untouched: + assert_array_equal(ag.positions, orig_pos) + @staticmethod @pytest.fixture() def ag(universe): From 415c44377ba81a1d63ba874b6ae67f1018e6aa57 Mon Sep 17 00:00:00 2001 From: zeman Date: Mon, 4 Feb 2019 16:36:31 +0100 Subject: [PATCH 12/21] make_whole() test for double precision box --- testsuite/MDAnalysisTests/lib/test_util.py | 23 ++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/testsuite/MDAnalysisTests/lib/test_util.py b/testsuite/MDAnalysisTests/lib/test_util.py index 5764422f88b..605d80a0dca 100644 --- a/testsuite/MDAnalysisTests/lib/test_util.py +++ b/testsuite/MDAnalysisTests/lib/test_util.py @@ -298,6 +298,29 @@ def test_out_of_place(self, universe): # positions must be untouched: assert_array_equal(ag.positions, orig_pos) + def test_double_precision_box(self): + # universe with double precision box containing a 2-atom molecule + # broken accross a corner: + u = mda.Universe.empty( + n_atoms=2, + n_residues=1, + n_segments=1, + atom_resindex=[0, 0], + residue_segindex=[0], + trajectory=True, + velocities=False, + forces=False) + ts = u.trajectory.ts + ts.frame = 0 + ts.dimensions = [10, 10, 10, 90, 90, 90] + assert ts.dimensions.dtype == np.float64 + ts.positions = np.array([[1, 1, 1,], [9, 9, 9]], dtype=np.float32) + u.add_TopologyAttr(Bonds([(0, 1)])) + mdamath.make_whole(u.atoms) + assert_array_almost_equal(u.atoms.positions, + np.array([[1, 1, 1,], [-1, -1, -1]], + dtype=np.float32)) + @staticmethod @pytest.fixture() def ag(universe): From 8d8dfe5dfada216cc631160d44f498367f37e356 Mon Sep 17 00:00:00 2001 From: zeman Date: Mon, 4 Feb 2019 17:43:29 +0100 Subject: [PATCH 13/21] changed dtype of molnums from int to np.int64 --- package/MDAnalysis/core/topologyattrs.py | 2 +- testsuite/MDAnalysisTests/topology/test_tprparser.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 2f7032344d3..a9f451b6c79 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -1460,7 +1460,7 @@ class Molnums(ResidueAttr): attrname = 'molnums' singular = 'molnum' target_classes = [AtomGroup, ResidueGroup, Atom, Residue] - dtype = int + dtype = np.int64 # segment attributes diff --git a/testsuite/MDAnalysisTests/topology/test_tprparser.py b/testsuite/MDAnalysisTests/topology/test_tprparser.py index 47805a435c6..9fb46ea81e3 100644 --- a/testsuite/MDAnalysisTests/topology/test_tprparser.py +++ b/testsuite/MDAnalysisTests/topology/test_tprparser.py @@ -54,6 +54,7 @@ def test_moltypes(self, top): def test_molnums(self, top): molnums = top.molnums.values assert_equal(molnums, self.ref_molnums) + assert molnums.dtype == np.int64 class TestTPR(TPRAttrs): From de53c87841d33ba0c8b94819d66a72111a9fcabc Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 02:30:17 +0100 Subject: [PATCH 14/21] fixed typo and removed unneccessary "elif"s to avoid CodeCov partials --- package/MDAnalysis/core/groups.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 6f36fdb7c84..1ec445bcaa2 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1258,7 +1258,7 @@ def wrap(self, compound="atoms", center="com", box=None): 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 accross periodic boundaries. + 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 @@ -1362,7 +1362,7 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): "mass of the group is zero.") refpos = np.sum(positions * masses[:, None], axis=0) refpos /= total_mass - elif ref == 'cog': + else: # ref == 'cog' refpos = positions.mean(axis=0) target = distances.apply_PBC(refpos, self.dimensions) positions += target - refpos @@ -1400,7 +1400,7 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): refpos = np.sum(positions[mask] * masses[:, None], axis=0) refpos /= total_mass - elif ref == 'cog': + else: # ref == 'cog' refpos = positions[mask].mean(axis=0) target = distances.apply_PBC(refpos, self.dimensions) positions[mask] += target - refpos From 5bce7c047215deaeb4942c4c3030b8cd4fcecdb1 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 02:30:52 +0100 Subject: [PATCH 15/21] removed unnecessary import --- package/MDAnalysis/core/topologyattrs.py | 1 - 1 file changed, 1 deletion(-) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index a9f451b6c79..21045cd3088 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -52,7 +52,6 @@ from ..lib.util import (cached, convert_aa_code, iterable, warn_if_not_unique, unique_int_1d) from ..lib import transformations, mdamath -from ..lib.distances import apply_PBC from ..exceptions import NoDataError, SelectionError from .topologyobjects import TopologyGroup from . import selection From 9599c5edfb5012fcfa9c0c6503c18fb15f50ead6 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 02:31:53 +0100 Subject: [PATCH 16/21] use libc.math.fabs instead of abs in make_whole() --- package/MDAnalysis/lib/_cutil.pyx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index e100c160fc0..c861d8c8d33 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/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 @@ -235,7 +235,7 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): is_unwrapped = True for i in range(1, natoms): for j in range(3): - if abs(oldpos[i, j] - oldpos[0, j]) >= half_box[j]: + if fabs(oldpos[i, j] - oldpos[0, j]) >= half_box[j]: is_unwrapped = False break if not is_unwrapped: From a7e8fcced2f9192240ed86c1a9848efb6fa3fa35 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 02:35:06 +0100 Subject: [PATCH 17/21] unwrap test for universe without masses --- testsuite/MDAnalysisTests/core/test_unwrap.py | 50 +++++++++++++++---- 1 file changed, 39 insertions(+), 11 deletions(-) diff --git a/testsuite/MDAnalysisTests/core/test_unwrap.py b/testsuite/MDAnalysisTests/core/test_unwrap.py index f53716c0450..9327cb1d621 100644 --- a/testsuite/MDAnalysisTests/core/test_unwrap.py +++ b/testsuite/MDAnalysisTests/core/test_unwrap.py @@ -24,11 +24,7 @@ from six.moves import range import numpy as np -from numpy.testing import ( - assert_equal, - assert_almost_equal, - assert_array_equal -) +from numpy.testing import assert_almost_equal, assert_array_equal import pytest import MDAnalysis as mda @@ -36,9 +32,6 @@ from MDAnalysis import NoDataError from MDAnalysis.lib.mdamath import triclinic_box -from MDAnalysisTests import make_Universe -from MDAnalysisTests.datafiles import TPR, XTC - def unwrap_test_universe(have_bonds=True, have_masses=True, have_molnums=True, is_triclinic=False): @@ -737,10 +730,45 @@ def test_unwrap_wrong_compound_exception_safety(self, level, reference): assert_array_equal(group.atoms.positions, orig_pos) @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) - def test_unwrap_no_bonds_exception_safety(self, level): + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + def test_unwrap_no_masses_exception_safety(self, level, compound): + # universe without masses: + u = unwrap_test_universe(have_masses=False) + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms + # select topology level: + if level == 'residues': + group = group.residues + elif level == 'segments': + group = group.segments + # store original positions: + orig_pos = group.atoms.positions + # try to unwrap: + with pytest.raises(NoDataError): + group.unwrap(compound=compound, reference='com', inplace=True) + # make sure atom positions are unchanged: + assert_array_equal(group.atoms.positions, orig_pos) + + @pytest.mark.parametrize('level', ('atoms', 'residues', 'segments')) + @pytest.mark.parametrize('compound', ('fragments', 'molecules', 'residues', + 'group', 'segments')) + @pytest.mark.parametrize('reference', ('com', 'cog', None)) + def test_unwrap_no_bonds_exception_safety(self, level, compound, reference): # universe without bonds: u = unwrap_test_universe(have_bonds=False) - group = u.atoms + # select group appropriate for compound: + if compound == 'group': + group = u.atoms[39:47] # molecule 12 + elif compound == 'segments': + group = u.atoms[23:47] # molecules 10, 11, 12 + else: + group = u.atoms # select topology level: if level == 'residues': group = group.residues @@ -749,7 +777,7 @@ def test_unwrap_no_bonds_exception_safety(self, level): # store original positions: orig_pos = group.atoms.positions with pytest.raises(NoDataError): - group.unwrap() + group.unwrap(compound=compound, reference=reference, inplace=True) # make sure atom positions are unchanged: assert_array_equal(group.atoms.positions, orig_pos) From 407879fba4bde1c21e8bd3e59a7d56c1a80bbfcf Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 03:11:18 +0100 Subject: [PATCH 18/21] added changelog entry for *Group.unwrap() --- package/CHANGELOG | 3 +++ 1 file changed, 3 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index e9f0ed60d2b..d7e72b1204d 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -27,6 +27,9 @@ 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) Changes * added official support for Python 3.7 (PR #1963) From f01a87892e0094dc34eb336f0ae0577dde4e9e03 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 03:17:10 +0100 Subject: [PATCH 19/21] added changelog entry for lib.mdamath.make_whole() --- package/CHANGELOG | 2 ++ 1 file changed, 2 insertions(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index d7e72b1204d..71a16c96617 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -30,6 +30,8 @@ Enhancements * 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) From 1c83cf324ddc81132b6554e2d7534d9fe484ed17 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 04:18:42 +0100 Subject: [PATCH 20/21] docs candy --- package/MDAnalysis/core/groups.py | 10 ++++++++-- package/MDAnalysis/lib/_cutil.pyx | 7 ++++++- 2 files changed, 14 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 1ec445bcaa2..164dadc45a9 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1214,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 """ @@ -1320,8 +1326,8 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): See Also -------- :func:`~MDAnalysis.lib.mdamath.make_whole`, - :meth:`~MDAnalysis.AtomGroup.wrap`, - :meth:`~MDAnalysis.AtomGroup.pack_into_box`, + :meth:`wrap`, + :meth:`pack_into_box`, :func:`~MDanalysis.lib.distances.apply_PBC` diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx index c861d8c8d33..5b5c99e3faf 100644 --- a/package/MDAnalysis/lib/_cutil.pyx +++ b/package/MDAnalysis/lib/_cutil.pyx @@ -156,7 +156,7 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): ------- Make fragments whole:: - from MDAnalysis.lib.util import make_whole + from MDAnalysis.lib.mdamath import make_whole # This algorithm requires bonds, these can be guessed! u = mda.Universe(......, guess_bonds=True) @@ -175,6 +175,11 @@ def make_whole(atomgroup, reference_atom=None, inplace=True): 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 From 986e59fea84394adeb66825d1417772df445cf77 Mon Sep 17 00:00:00 2001 From: zeman Date: Wed, 6 Feb 2019 04:53:42 +0100 Subject: [PATCH 21/21] eliminated another useless elif --- package/MDAnalysis/core/groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/core/groups.py b/package/MDAnalysis/core/groups.py index 164dadc45a9..03bed317372 100644 --- a/package/MDAnalysis/core/groups.py +++ b/package/MDAnalysis/core/groups.py @@ -1380,7 +1380,7 @@ def unwrap(self, compound='fragments', reference='com', inplace=True): compound_indices = unique_atoms.resindices elif comp == 'segments': compound_indices = unique_atoms.segindices - elif comp == 'molecules': + else: # comp == 'molecules' try: compound_indices = unique_atoms.molnums except AttributeError: