diff --git a/package/.pylintrc b/package/.pylintrc index ab966147016..da780ce8ac3 100644 --- a/package/.pylintrc +++ b/package/.pylintrc @@ -200,7 +200,6 @@ enable=abstract-class-instantiated, xrange-builtin, zip-builtin-not-iterating, map-builtin-not-iterating, - range-builtin-not-iterating, # Things we'd like to try. # Procedure: diff --git a/package/CHANGELOG b/package/CHANGELOG index 05e9b078d1c..c491f708128 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,12 +14,20 @@ The rules for this file: ------------------------------------------------------------------------------ ??/??/18 tylerjereddy, richardjgowers, palnabarun, orbeckst, kain88-de, zemanj, - VOD555, davidercruz, jbarnoud, ayushsuhane, hfmull + VOD555, davidercruz, jbarnoud, ayushsuhane, hfmull, micaela-matta * 0.18.1 Enhancements + * Added wrap/unwrap transformations + * Modified around selections to work with KDTree and periodic boundary + conditions. Should reduce memory usage (#974 PR #2022) + * Modified topology.guessers.guess_bonds to automatically select the + fastest method for guessing bonds using + lib.distance.self_capped_distance (PR # 2006) + * Added lib.distances.self_capped_distance to internally select the + optimized method for distance evaluations of coordinates with itself. (PR # 2006) * Added augment functionality to create relevant images of particles in the vicinity of central cell to handle periodic boundary conditions (PR #1977) @@ -68,6 +76,7 @@ Fixes * PDBWriter now properly sets start value * Zero length TopologyGroup now return proper shape array via to_indices (Issue #1974) + * Added periodic boundary box support to the Tinker file reader (Issue #2001) Changes * TopologyAttrs are now statically typed (Issue #1876) diff --git a/package/MDAnalysis/coordinates/TXYZ.py b/package/MDAnalysis/coordinates/TXYZ.py index 260d11d6ead..aa604b96ad0 100644 --- a/package/MDAnalysis/coordinates/TXYZ.py +++ b/package/MDAnalysis/coordinates/TXYZ.py @@ -74,12 +74,19 @@ def __init__(self, filename, **kwargs): root, ext = os.path.splitext(self.filename) self.xyzfile = util.anyopen(self.filename) self._cache = dict() - + # Check if file has box information saved + with util.openany(self.filename) as inp: + inp.readline() + line = inp.readline() + # If second line has float at second position, we have box info + try: + float(line.split()[1]) + except ValueError: + self.periodic = False + else: + self.periodic = True self.ts = self._Timestep(self.n_atoms, **self._ts_kwargs) - # Haven't quite figured out where to start with all the self._reopen() - # etc. - # (Also cannot just use seek() or reset() because that would break - # with urllib2.urlopen() streams) + self._read_next_timestep() @property @@ -103,6 +110,8 @@ def _read_xyz_n_frames(self): # the number of lines in the XYZ file will be 1 greater than the # number of atoms linesPerFrame = self.n_atoms + 1 + if self.periodic: + linesPerFrame += 1 counter = 0 offsets = [] @@ -134,6 +143,8 @@ def _read_next_timestep(self, ts=None): try: # we assume that there is only one header line per frame f.readline() + if self.periodic: + ts.dimensions = f.readline().split() # convert all entries at the end once for optimal speed tmp_buf = [] for i in range(self.n_atoms): diff --git a/package/MDAnalysis/core/selection.py b/package/MDAnalysis/core/selection.py index bc4a256e20e..4532ac854a7 100644 --- a/package/MDAnalysis/core/selection.py +++ b/package/MDAnalysis/core/selection.py @@ -243,9 +243,6 @@ def __init__(self): self.apply = self._apply_distmat self.periodic = flags['use_periodic_selections'] - # KDTree doesn't support periodic - if self.periodic: - self.apply = self._apply_distmat def validate_dimensions(self, dimensions): r"""Check if the system is periodic in all three-dimensions. @@ -282,6 +279,9 @@ def _apply_KDTree(self, group): # All atoms in group that aren't in sel sys = group[~np.in1d(group.indices, sel.indices)] + if not sys: + return sys[[]] + box = self.validate_dimensions(group.dimensions) cut = self.cutoff if box is not None else None @@ -321,7 +321,8 @@ def _apply_KDTree(self, group): """ sel = self.sel.apply(group) box = self.validate_dimensions(group.dimensions) - ref = sel.center_of_geometry(pbc=self.periodic) + periodic = box is not None + ref = sel.center_of_geometry(pbc=periodic) kdtree = PeriodicKDTree(box=box) cutoff = self.exRadius if box is not None else None @@ -363,7 +364,8 @@ def _apply_KDTree(self, group): """ sel = self.sel.apply(group) box = self.validate_dimensions(group.dimensions) - ref = sel.center_of_geometry(pbc=self.periodic) + periodic = box is not None + ref = sel.center_of_geometry(pbc=periodic) cut = self.cutoff if box is not None else None kdtree = PeriodicKDTree(box=box) @@ -483,7 +485,7 @@ def __init__(self, parser, tokens): self.cutoff = float(tokens.popleft()) def _apply_KDTree(self, group): - box = group.dimensions if self.periodic else None + box = self.validate_dimensions(group.dimensions) kdtree = PeriodicKDTree(box=box) cut = self.cutoff if box is not None else None kdtree.set_coords(group.positions, cutoff=cut) @@ -496,7 +498,7 @@ def _apply_distmat(self, group): ref_coor = self.ref[np.newaxis, ...] ref_coor = np.asarray(ref_coor, dtype=np.float32) - box = group.dimensions if self.periodic else None + box = self.validate_dimensions(group.dimensions) dist = distances.distance_array(group.positions, ref_coor, box) mask = (dist <= self.cutoff).any(axis=1) diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py index 57a5625b074..1cb86147235 100644 --- a/package/MDAnalysis/lib/distances.py +++ b/package/MDAnalysis/lib/distances.py @@ -59,6 +59,8 @@ .. autofunction:: calc_angles(atom1, atom2, atom3 [,box [, result [, backend]]]) .. autofunction:: calc_dihedrals(atom1, atom2, atom3, atom4 [,box [, result [, backend]]]) .. autofunction:: apply_PBC(coordinates, box [, backend]) +.. autofunction:: capped_distance(reference, configuration, max_cutoff [, min_cutoff [, box [, method]]]) +.. autofunction:: self_capped_distance(reference, max_cutoff, [, min_cutoff [, box [, method]]]) .. autofunction:: transform_RtoS(coordinates, box [, backend]) .. autofunction:: transform_StoR(coordinates, box [,backend]) .. autofunction:: MDAnalysis.lib._augment.augment_coordinates(coordinates, box, radius) @@ -400,8 +402,9 @@ def self_distance_array(reference, box=None, result=None, backend="serial"): return distances -def capped_distance(reference, configuration, max_cutoff, min_cutoff=None, box=None, method=None): - """Calculates the pairs and distance within a specified distance +def capped_distance(reference, configuration, max_cutoff, min_cutoff=None, + box=None, method=None): + """Calculates the pairs and distances within a specified distance If a *box* is supplied, then a minimum image convention is used to evaluate the distances. @@ -460,7 +463,6 @@ def capped_distance(reference, configuration, max_cutoff, min_cutoff=None, box=N .. SeeAlso:: :func:'MDAnalysis.lib.distances.distance_array' .. SeeAlso:: :func:'MDAnalysis.lib.pkdtree.PeriodicKDTree' - """ if box is not None: if box.shape[0] != 6: @@ -483,12 +485,38 @@ def _determine_method(reference, configuration, max_cutoff, min_cutoff=None, All the rules to select the method based on the input can be incorporated here. + Parameters + ---------- + reference : array + reference coordinates array with shape ``reference.shape = (3,)`` + or ``reference.shape = (len(reference), 3)`` + configuration : array + Configuration coordinate array with shape ``reference.shape = (3,)`` + or ``reference.shape = (len(reference), 3)`` + max_cutoff : float + Maximum cutoff distance between the reference and configuration + min_cutoff : (optional) float + Minimum cutoff distance between reference and configuration [None] + box : (optional) array or None + The dimensions, if provided, must be provided in the same + The unitcell dimesions for this system format as returned + by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: + ``[lx,ly, lz, alpha, beta, gamma]``. Minimum image convention + is applied if the box is provided [None] + method : (optional) 'bruteforce' or 'pkdtree' or 'None' + Keyword to override the automatic guessing of method built-in + in the function [None] + Returns ------- - Function object based on the rules and specified method - Currently implemented methods are - bruteforce : returns ``_bruteforce_capped`` - PKDtree : return ``_pkdtree_capped` + Method : Function object + Returns function object based on the rules and specified method + + Note + ---- + Currently implemented methods are present in the ``methods`` dictionary + bruteforce : returns ``_bruteforce_capped`` + PKDtree : return ``_pkdtree_capped` """ methods = {'bruteforce': _bruteforce_capped, @@ -498,37 +526,49 @@ def _determine_method(reference, configuration, max_cutoff, min_cutoff=None, return methods[method] if len(reference) > 5000 and len(configuration) > 5000: - if box is None and reference.shape[0] != 3 and configuration.shape[0] != 3: + if box is None: min_dim = np.array([reference.min(axis=0), configuration.min(axis=0)]) max_dim = np.array([reference.max(axis=0), configuration.max(axis=0)]) size = max_dim.max(axis=0) - min_dim.min(axis=0) - elif box is not None: - if np.allclose(box[3:], 90): - size = box[:3] - else: - tribox = triclinic_vectors(box) - size = tribox.max(axis=0) - tribox.min(axis=0) - - if (np.any(size < 10.0*max_cutoff) and + elif np.allclose(box[3:], 90): + size = box[:3] + else: + tribox = triclinic_vectors(box) + size = tribox.max(axis=0) - tribox.min(axis=0) + + if ((np.any(size < 10.0*max_cutoff) and len(reference) > 100000 and - len(configuration) > 100000): + len(configuration) > 100000)): return methods['bruteforce'] else: return methods['pkdtree'] return methods['bruteforce'] -def _bruteforce_capped(reference, configuration, max_cutoff, min_cutoff=None, box=None): - """ - Using naive distance calulations, returns a list +def _bruteforce_capped(reference, configuration, max_cutoff, + min_cutoff=None, box=None): + """Internal method for bruteforce calculations + + Uses naive distance calulations and returns a list containing the indices with one from each reference and configuration arrays, such that the distance between them is less than the specified cutoff distance + + Returns + ------- + pairs : list + List of ``[(i, j)]`` pairs such that atom-index ``i`` is + from reference and ``j`` from configuration array + distance: list + Distance between ``reference[i]`` and ``configuration[j]`` + atom coordinate + """ pairs, distance = [], [] + reference = np.asarray(reference, dtype=np.float32) configuration = np.asarray(configuration, dtype=np.float32) @@ -543,18 +583,21 @@ def _bruteforce_capped(reference, configuration, max_cutoff, min_cutoff=None, bo for i, coords in enumerate(reference): dist = distance_array(coords[None, :], configuration, box=box)[0] if min_cutoff is not None: - idx = np.where((dist <= max_cutoff) & (dist > min_cutoff))[0] + idx = np.where((dist < max_cutoff) & (dist > min_cutoff))[0] else: - idx = np.where((dist <= max_cutoff))[0] + idx = np.where((dist < max_cutoff))[0] for j in idx: pairs.append((i, j)) distance.append(dist[j]) return pairs, distance -def _pkdtree_capped(reference, configuration, max_cutoff, min_cutoff=None, box=None): +def _pkdtree_capped(reference, configuration, max_cutoff, + min_cutoff=None, box=None): """ Capped Distance evaluations using KDtree. + Uses minimum image convention if *box* is specified + Returns: -------- pairs : list @@ -564,6 +607,7 @@ def _pkdtree_capped(reference, configuration, max_cutoff, min_cutoff=None, box=N distance : list Distance between two atoms corresponding to the (i, j) indices in pairs. + """ from .pkdtree import PeriodicKDTree @@ -600,6 +644,220 @@ def _pkdtree_capped(reference, configuration, max_cutoff, min_cutoff=None, box=N return pairs, distances +def self_capped_distance(reference, max_cutoff, min_cutoff=None, + box=None, method=None): + """Finds all the pairs and respective distances within a specified cutoff + for a configuration *reference* + + If a *box* is supplied, then a minimum image convention is used + to evaluate the distances. + + An automatic guessing of optimized method to calculate the distances is + included in the function. An optional keyword for the method is also + provided. Users can override the method with this functionality. + Currently pkdtree and bruteforce are implemented. + + Parameters + ----------- + reference : array + reference coordinates array with shape ``reference.shape = (3,)`` + or ``reference.shape = (len(reference), 3)`` + max_cutoff : float + Maximum cutoff distance to check the neighbors with itself + min_cutoff : (optional) float + Minimum cutoff distance [None] + box : (optional) array or None + The dimensions, if provided, must be provided in the same + The unitcell dimesions for this system format as returned + by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: + ``[lx,ly, lz, alpha, beta, gamma]``. Minimum image convention + is applied if the box is provided [None] + method : (optional) 'bruteforce' or 'pkdtree' or 'None' + Keyword to override the automatic guessing of method built-in + in the function [None] + + Returns + ------- + pairs : array + Pair of indices such that distance between them is + within the ``max_cutoff`` and ``min_cutoff`` + distances : array + Distances corresponding to each pair of indices. + d[k] corresponding to the pairs[i,j] gives the distance between + i-th and j-th coordinate in reference + + .. code-block:: python + + pairs, distances = self_capped_distances(reference, max_cutoff) + for indx, [a,b] in enumerate(pairs): + coord1, coords2 = reference[a], reference[b] + distance = distances[indx] + + Note + ----- + Currently only supports brute force and Periodic KDtree + + .. SeeAlso:: :func:'MDAnalysis.lib.distances.self_distance_array' + .. SeeAlso:: :func:'MDAnalysis.lib.pkdtree.PeriodicKDTree' + """ + if box is not None: + if box.shape[0] != 6: + raise ValueError('Box Argument is of incompatible type. The dimension' + 'should be either None or ' + 'of the type [lx, ly, lz, alpha, beta, gamma]') + method = _determine_method_self(reference, max_cutoff, + min_cutoff=min_cutoff, + box=box, method=method) + pairs, dist = method(reference, max_cutoff, + min_cutoff=min_cutoff, box=box) + + return np.asarray(pairs), np.asarray(dist) + + +def _determine_method_self(reference, max_cutoff, min_cutoff=None, + box=None, method=None): + """ + Switch between different methods based on the the optimized time. + All the rules to select the method based on the input can be + incorporated here. + + Parameters + ---------- + reference : array + reference coordinates array with shape ``reference.shape = (3,)`` + or ``reference.shape = (len(reference), 3)`` + max_cutoff : float + Maximum cutoff distance + min_cutoff : (optional) float + Minimum cutoff distance [None] + box : (optional) array or None + The dimensions, if provided, must be provided in the same + The unitcell dimesions for this system format as returned + by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`: + ``[lx,ly, lz, alpha, beta, gamma]``. Minimum image convention + is applied if the box is provided [None] + method : (optional) 'bruteforce' or 'pkdtree' or 'None' + Keyword to override the automatic guessing of method built-in + in the function [None] + + Returns + ------- + Method : Function object + Returns function object based on the rules and specified method + + Note + ---- + Currently implemented methods are present in the ``methods`` dictionary + bruteforce : returns ``_bruteforce_capped_self`` + PKDtree : return ``_pkdtree_capped_self`` + + """ + methods = {'bruteforce': _bruteforce_capped_self, + 'pkdtree': _pkdtree_capped_self} + + if method is not None: + return methods[method] + + if len(reference) > 5000: + if box is None: + min_dim = np.array([reference.min(axis=0)]) + max_dim = np.array([reference.max(axis=0)]) + size = max_dim.max(axis=0) - min_dim.min(axis=0) + elif np.allclose(box[3:], 90): + size = box[:3] + else: + tribox = triclinic_vectors(box) + size = tribox.max(axis=0) - tribox.min(axis=0) + + if ((np.any(size < 10.0*max_cutoff) and + (len(reference) > 100000))): + return methods['bruteforce'] + else: + return methods['pkdtree'] + return methods['bruteforce'] + + +def _bruteforce_capped_self(reference, max_cutoff, min_cutoff=None, + box=None): + """Finds all the pairs among the *reference* coordinates within + a fixed distance using brute force method + + Internal method using brute force method to evaluate all the pairs + of atoms within a fixed distance. + + Returns + ------- + pairs : array + Arrray of ``[i, j]`` pairs such that atom-index ``i`` + and ``j`` from reference array are within the fixed distance + distance: array + Distance between ``reference[i]`` and ``reference[j]`` + atom coordinate + + """ + pairs, distance = [], [] + + reference = np.asarray(reference, dtype=np.float32) + if reference.shape == (3, ): + reference = reference[None, :] + for i, coords in enumerate(reference): + # Each pair of atoms needs to be checked only once. + # Only calculate distance for atomA and atomB + # if atomidA < atomidB + dist = distance_array(coords[None, :], reference[i+1:], + box=box)[0] + + if min_cutoff is not None: + idx = np.where((dist < max_cutoff) & (dist > min_cutoff))[0] + else: + idx = np.where((dist < max_cutoff))[0] + for other_idx in idx: + # Actual atomid for atomB + # can be direclty obtained in this way + j = other_idx + 1 + i + pairs.append((i, j)) + distance.append(dist[other_idx]) + return np.asarray(pairs), np.asarray(distance) + + +def _pkdtree_capped_self(reference, max_cutoff, min_cutoff=None, + box=None): + """Finds all the pairs among the coordinates within a fixed distance + using PeriodicKDTree + + Internal method using PeriodicKDTree method to evaluate all the pairs + of atoms within a fixed distance. + + Returns + ------- + pairs : array + Array of ``[(i, j)]`` pairs such that atom-index ``i`` + and ``j`` from reference array are within the fixed distance + distance: array + Distance between ``reference[i]`` and ``reference[j]`` + atom coordinate + + """ + from .pkdtree import PeriodicKDTree + + reference = np.asarray(reference, dtype=np.float32) + if reference.shape == (3, ): + reference = reference[None, :] + + pairs, distance = [], [] + kdtree = PeriodicKDTree(box=box) + cut = max_cutoff if box is not None else None + kdtree.set_coords(reference, cutoff=cut) + pairs = kdtree.search_pairs(max_cutoff) + if pairs.size > 0: + refA, refB = pairs[:, 0], pairs[:, 1] + distance = calc_bonds(reference[refA], reference[refB], box=box) + if min_cutoff is not None: + mask = np.where(distance > min_cutoff)[0] + pairs, distance = pairs[mask], distance[mask] + return np.asarray(pairs), np.asarray(distance) + + def transform_RtoS(inputcoords, box, backend="serial"): """Transform an array of coordinates from real space to S space (aka lambda space) diff --git a/package/MDAnalysis/lib/formats/src/xdrfile.c b/package/MDAnalysis/lib/formats/src/xdrfile.c index 470bd26fe76..19a1c0b93b4 100644 --- a/package/MDAnalysis/lib/formats/src/xdrfile.c +++ b/package/MDAnalysis/lib/formats/src/xdrfile.c @@ -45,6 +45,7 @@ #endif #include +#define _FILE_OFFSET_BITS 64 #ifdef _WIN32 # include @@ -2518,8 +2519,8 @@ static int xdrstdio_getlong (XDR *, int32_t *); static int xdrstdio_putlong (XDR *, int32_t *); static int xdrstdio_getbytes (XDR *, char *, unsigned int); static int xdrstdio_putbytes (XDR *, char *, unsigned int); -static off_t xdrstdio_getpos (XDR *); -static int xdrstdio_setpos (XDR *, off_t, int); +static int64_t xdrstdio_getpos (XDR *); +static int xdrstdio_setpos (XDR *, int64_t, int); static void xdrstdio_destroy (XDR *); /* @@ -2601,7 +2602,7 @@ xdrstdio_putbytes (XDR *xdrs, char *addr, unsigned int len) } -static off_t +static int64_t xdrstdio_getpos (XDR *xdrs) { #ifdef _WIN32 @@ -2612,7 +2613,7 @@ xdrstdio_getpos (XDR *xdrs) } static int -xdrstdio_setpos (XDR *xdrs, off_t pos, int whence) +xdrstdio_setpos (XDR *xdrs, int64_t pos, int whence) { /* A reason for failure can be filesystem limits on allocation units, * before the actual off_t overflow (ext3, with a 4K clustersize, @@ -2640,7 +2641,7 @@ int xdr_seek(XDRFILE *xd, int64_t pos, int whence) /* Seeks to position in file */ { int result; - if ((result = xdrstdio_setpos(xd->xdr, (off_t) pos, whence)) != 0) + if ((result = xdrstdio_setpos(xd->xdr, (int64_t) pos, whence)) != 0) return result; return exdrOK; diff --git a/package/MDAnalysis/lib/pkdtree.py b/package/MDAnalysis/lib/pkdtree.py index 8ab562dbd45..5b1b21594f2 100644 --- a/package/MDAnalysis/lib/pkdtree.py +++ b/package/MDAnalysis/lib/pkdtree.py @@ -192,7 +192,7 @@ def search(self, centers, radius): radius)) self._indices = np.array(list( itertools.chain.from_iterable(indices)), - dtype=np.int) + dtype=np.int64) self._indices = np.asarray(unique_int_1d(self._indices)) return self._indices diff --git a/package/MDAnalysis/topology/TXYZParser.py b/package/MDAnalysis/topology/TXYZParser.py index 63698602d06..74bfdce4513 100644 --- a/package/MDAnalysis/topology/TXYZParser.py +++ b/package/MDAnalysis/topology/TXYZParser.py @@ -43,7 +43,10 @@ """ from __future__ import absolute_import + +import itertools import numpy as np +from six.moves import zip from . import guessers from ..lib.util import openany @@ -88,9 +91,22 @@ def parse(self, **kwargs): names = np.zeros(natoms, dtype=object) types = np.zeros(natoms, dtype=np.int) bonds = [] + # Find first atom line, maybe there's box information + fline = inf.readline() + try: + # If a box second value will be a float + # If an atom, second value will be a string + float(fline.split()[1]) + except ValueError: + # If float conversion failed, we have first atom line + pass + else: + # If previous try succeeded it was a box + # so read another line to find the first atom line + fline = inf.readline() # Can't infinitely read as XYZ files can be multiframe - for i in range(natoms): - line = inf.readline().split() + for i, line in zip(range(natoms), itertools.chain([fline], inf)): + line = line.split() atomids[i]= line[0] names[i] = line[1] types[i] = line[5] diff --git a/package/MDAnalysis/topology/guessers.py b/package/MDAnalysis/topology/guessers.py index 06689168079..bed9cca8b05 100644 --- a/package/MDAnalysis/topology/guessers.py +++ b/package/MDAnalysis/topology/guessers.py @@ -232,25 +232,14 @@ def guess_bonds(atoms, coords, box=None, **kwargs): bonds = [] - for i, atom in enumerate(atoms[:-1]): - vdw_i = vdwradii[atomtypes[i]] - max_d = (vdw_i + max_vdw) * fudge_factor - - # using self_distance_array scales O(n^2) - # 20,000 atoms = 1.6 Gb memory - dist = distances.distance_array(coords[i][None, :], coords[i + 1:], - box=box)[0] - idx = np.where((dist > lower_bound) & (dist <= max_d))[0] - - for a in idx: - j = i + 1 + a - atom_j = atoms[j] - - if dist[a] < (vdw_i + vdwradii[atomtypes[j]]) * fudge_factor: - # because of method used, same bond won't be seen twice, - # so don't need to worry about duplicates - bonds.append((atom.index, atom_j.index)) - + pairs, dist = distances.self_capped_distance(coords, + max_cutoff=2.0*max_vdw, + min_cutoff=lower_bound, + box=box) + for idx, (i, j) in enumerate(pairs): + d = (vdwradii[atomtypes[i]] + vdwradii[atomtypes[j]])*fudge_factor + if (dist[idx] < d): + bonds.append((atoms[i].index, atoms[j].index)) return tuple(bonds) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index bcefa903e2c..62488f7977a 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -94,3 +94,4 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby +from .wrap import wrap, unwrap diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index a9712d2a44a..ce6df7b8e26 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -21,12 +21,14 @@ # """\ -Rotate trajectory --- :mod:`MDAnalysis.transformations.translate` -================================================================= +Trajectory rotation --- :mod:`MDAnalysis.transformations.rotate` +================================================================ Rotates the coordinates by a given angle arround an axis formed by a direction and a point - + +.. autofunction:: rotateby + """ from __future__ import absolute_import @@ -35,21 +37,28 @@ from functools import partial from ..lib.transformations import rotation_matrix +from ..lib.util import get_weights -def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=None): +def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): ''' Rotates the trajectory by a given angle on a given axis. The axis is defined by the user, combining the direction vector and a point. This point can be the center - of geometry or the center of mass of a user defined AtomGroup, or a list defining custom - coordinates. - e.g. rotate the coordinates by 90 degrees on a x axis centered on a given atom group: + of geometry or the center of mass of a user defined AtomGroup, or an array defining + custom coordinates. + + Examples + -------- + + e.g. rotate the coordinates by 90 degrees on a axis formed by the [0,0,1] vector and + the center of geometry of a given AtomGroup: .. code-block:: python ts = u.trajectory.ts angle = 90 ag = u.atoms() - rotated = MDAnalysis.transformations.rotate(angle, ag)(ts) + d = [0,0,1] + rotated = MDAnalysis.transformations.rotate(angle, direction=d, ag=ag)(ts) e.g. rotate the coordinates by a custom axis: @@ -59,7 +68,7 @@ def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=Non angle = 90 p = [1,2,3] d = [0,0,1] - rotated = MDAnalysis.transformations.rotate(angle, point=point, direction=d)(ts) + rotated = MDAnalysis.transformations.rotate(angle, direction=d, point=point)(ts) Parameters ---------- @@ -67,24 +76,28 @@ def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=Non rotation angle in degrees direction: array-like vector that will define the direction of a custom axis of rotation from the - provided point. + provided point. Expected shapes are (3, ) or (1, 3). ag: AtomGroup, optional - use this to define the center of mass or geometry as the point from where the - rotation axis will be defined - center: str, optional - used to choose the method of centering on the given atom group. Can be 'geometry' - or 'mass' + use the weighted center of an AtomGroup as the point from where the rotation axis + will be defined. If no AtomGroup is given, the `point` argument becomes mandatory + point: array-like, optional + list of the coordinates of the point from where a custom axis of rotation will + be defined. Expected shapes are (3, ) or (1, 3). If no point is given, the + `ag` argument becomes mandatory. + weights: {"mass", ``None``} or array_like, optional + define the weights of the atoms when calculating the center of the AtomGroup. + With ``"mass"`` uses masses as weights; with ``None`` weigh each atom equally. + If a float array of the same length as `ag` is provided, use each element of + the `array_like` as a weight for the corresponding atom in `ag`. Default is + None. wrap: bool, optional If `True`, all the atoms from the given AtomGroup will be moved to the unit cell before calculating the center of mass or geometry. Default is `False`, no changes to the atom coordinates are done before calculating the center of the AtomGroup. - point: array-like, optional - list of the coordinates of the point from where a custom axis of rotation will - be defined. Returns ------- - :class:`~MDAnalysis.coordinates.base.Timestep` object + MDAnalysis.coordinates.base.Timestep Warning ------- @@ -92,25 +105,31 @@ def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=Non after rotating the trajectory. ''' - pbc_arg = wrap angle = np.deg2rad(angle) - if point: + try: + direction = np.asarray(direction, np.float32) + if direction.shape != (3, ) and direction.shape != (1, 3): + raise ValueError('{} is not a valid direction'.format(direction)) + direction = direction.reshape(3, ) + except ValueError: + raise ValueError('{} is not a valid direction'.format(direction)) + if point is not None: point = np.asarray(point, np.float32) if point.shape != (3, ) and point.shape != (1, 3): raise ValueError('{} is not a valid point'.format(point)) + point = point.reshape(3, ) elif ag: try: - if center == 'geometry': - center_method = partial(ag.center_of_geometry, pbc=pbc_arg) - elif center == 'mass': - center_method = partial(ag.center_of_mass, pbc=pbc_arg) - else: - raise ValueError('{} is not a valid argument for center'.format(center)) + atoms = ag.atoms except AttributeError: - if center == 'mass': - raise AttributeError('{} is not an AtomGroup object with masses'.format(ag)) - else: - raise ValueError('{} is not an AtomGroup object'.format(ag)) + raise ValueError('{} is not an AtomGroup object'.format(ag)) + else: + try: + weights = get_weights(atoms, weights=weights) + except (ValueError, TypeError): + raise TypeError("weights must be {'mass', None} or an iterable of the " + "same size as the atomgroup.") + center_method = partial(atoms.center, weights, pbc=wrap) else: raise ValueError('A point or an AtomGroup must be specified') @@ -119,8 +138,11 @@ def wrapped(ts): position = center_method() else: position = point - rotation = rotation_matrix(angle, direction, position)[:3, :3] + matrix = rotation_matrix(angle, direction, position) + rotation = matrix[:3, :3].T + translation = matrix[:3, 3] ts.positions= np.dot(ts.positions, rotation) + ts.positions += translation return ts diff --git a/package/MDAnalysis/transformations/wrap.py b/package/MDAnalysis/transformations/wrap.py new file mode 100644 index 00000000000..101591899f2 --- /dev/null +++ b/package/MDAnalysis/transformations/wrap.py @@ -0,0 +1,146 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# 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. +# +# 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 +# + +"""\ +Wrap/unwrap transformations --- :mod:`MDAnalysis.transformations.wrap` +====================================================================== + +Wrap/unwrap the atoms of a given AtomGroup in the unit cell. :func:`wrap` +translates the atoms back in the unit cell. :func:`unwrap` translates the +atoms of each molecule so that bons don't split over images. + +.. autofunction:: wrap + +.. autofunction:: unwrap + +""" + +from ..lib._cutil import make_whole + +def wrap(ag, compound='atoms'): + """ + Shift the contents of a given AtomGroup back into the unit cell. :: + + +-----------+ +-----------+ + | | | | + | 3 | 6 | 6 3 | + | ! | ! | ! ! | + | 1-2-|-5-8 -> |-5-8 1-2-| + | ! | ! | ! ! | + | 4 | 7 | 7 4 | + | | | | + +-----------+ +-----------+ + + Example + ------- + + .. code-block:: python + + ag = u.atoms + transform = mda.transformations.wrap(ag) + u.trajectory.add_transformations(transform) + + Parameters + ---------- + + ag: Atomgroup + Atomgroup to be wrapped in the unit cell + compound : {'atoms', 'group', 'residues', 'segments', 'fragments'}, optional + The group which will be kept together through the shifting process. + + Notes + ----- + When specifying a `compound`, the translation is calculated based on + each compound. The same translation is applied to all atoms + within this compound, meaning it will not be broken by the shift. + This might however mean that not all atoms from the compound are + inside the unit cell, but rather the center of the compound is. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + + """ + + def wrapped(ts): + ag.wrap(compound=compound) + + return ts + + return wrapped + + +def unwrap(ag): + """ + Move all atoms in an AtomGroup so that bonds don't split over images + + Atom positions are modified in place. + + 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 + | | | | + +-----------+ +-----------+ + + Example + ------- + + .. code-block:: python + + ag = u.atoms + transform = mda.transformations.unwrap(ag) + u.trajectory.add_transformations(transform) + + Parameters + ---------- + atomgroup : AtomGroup + The :class:`MDAnalysis.core.groups.AtomGroup` to work with. + The positions of this are modified in place. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + + """ + + try: + ag.fragments + except AttributeError: + raise AttributeError("{} has no fragments".format(ag)) + + def wrapped(ts): + for frag in ag.fragments: + make_whole(frag) + + return ts + + return wrapped + diff --git a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst index 71599ff44f9..42d16de3412 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -96,3 +96,7 @@ e.g. giving a workflow as a keyword argument when defining the universe: .. toctree:: ./transformations/translate + ./transformations/rotate + ./transformations/wrap + + diff --git a/package/doc/sphinx/source/documentation_pages/transformations/rotate.rst b/package/doc/sphinx/source/documentation_pages/transformations/rotate.rst new file mode 100644 index 00000000000..38bae418c86 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/transformations/rotate.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.transformations.rotate \ No newline at end of file diff --git a/package/doc/sphinx/source/documentation_pages/transformations/wrap.rst b/package/doc/sphinx/source/documentation_pages/transformations/wrap.rst new file mode 100644 index 00000000000..868288015a0 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/transformations/wrap.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.transformations.wrap \ No newline at end of file diff --git a/testsuite/MDAnalysisTests/coordinates/test_txyz.py b/testsuite/MDAnalysisTests/coordinates/test_txyz.py index 6334352645d..8dd49ac8d7a 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_txyz.py +++ b/testsuite/MDAnalysisTests/coordinates/test_txyz.py @@ -23,8 +23,9 @@ import pytest from numpy.testing import assert_almost_equal +from six.moves import zip -from MDAnalysisTests.datafiles import TXYZ, ARC +from MDAnalysisTests.datafiles import TXYZ, ARC, ARC_PBC import MDAnalysis as mda @@ -39,6 +40,11 @@ def ARC_U(): return mda.Universe(ARC) +@pytest.fixture +def ARC_PBC_U(): + return mda.Universe(ARC_PBC) + + def test_txyz_positions(TXYZ_U): assert_almost_equal(TXYZ_U.atoms.positions[0], [-6.553398, -1.854369, 0.000000]) @@ -56,5 +62,19 @@ def test_arc_positions_frame_2(ARC_U): [-0.231579, -0.350841, -0.037475]) -def test_arc_traj_legnth(ARC_U): +def test_arc_traj_length(ARC_U): assert len(ARC_U.trajectory) == 2 + + +def test_arcpbc_traj_length(ARC_PBC_U): + assert len(ARC_PBC_U.trajectory) == 3 + + +def test_pbc_boxsize(ARC_PBC_U): + ref_dimensions=[[ 38.860761, 38.860761, 38.860761, 90.000000, 90.000000, 90.000000], + [ 39.860761, 39.860761, 39.860761, 90.000000, 90.000000, 90.000000], + [ 40.860761, 40.860761, 40.860761, 90.000000, 90.000000, 90.000000]] + + for ref_box, ts in zip(ref_dimensions, ARC_PBC_U.trajectory): + assert_almost_equal(ref_box, ts.dimensions, decimal=5) + diff --git a/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc b/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc new file mode 100644 index 00000000000..6705ae636a9 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc @@ -0,0 +1,24 @@ + 6 Hexane (267 OPLS-UA in 38.8 Ang Box; JPC, 100, 14511 '96) + 38.860761 38.860761 38.860761 90.000000 90.000000 90.000000 + 1 CH3 -0.533809 8.893843 -17.718901 83 2 + 2 CH2 -0.492066 8.734009 -16.201869 86 1 3 + 3 CH2 0.609320 7.723360 -15.894925 86 2 4 + 4 CH2 0.723163 7.301395 -14.432851 86 3 5 + 5 CH2 1.923300 6.401842 -14.151512 86 4 6 + 6 CH3 1.576645 4.928146 -14.343148 83 5 + 6 Hexane (test file for MDA) + 39.860761 39.860761 39.860761 90.000000 90.000000 90.000000 + 1 CH3 0.039519 9.187867 -17.899888 83 2 + 2 CH2 -0.171280 8.486514 -16.561104 86 1 3 + 3 CH2 0.988507 7.632303 -16.057222 86 2 4 + 4 CH2 0.630146 7.086837 -14.677831 86 3 5 + 5 CH2 1.729866 6.240985 -14.042356 86 4 6 + 6 CH3 2.065117 4.899261 -14.687383 83 5 + 6 Hexane + 40.860761 40.860761 40.860761 90.000000 90.000000 90.000000 + 1 CH3 -0.533809 8.893843 -17.718901 83 2 + 2 CH2 -0.492066 8.734009 -16.201869 86 1 3 + 3 CH2 0.609320 7.723360 -15.894925 86 2 4 + 4 CH2 0.723163 7.301395 -14.432851 86 3 5 + 5 CH2 1.923300 6.401842 -14.151512 86 4 6 + 6 CH3 1.576645 4.928146 -14.343148 83 5 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 3109f100527..afc51555d0e 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -85,7 +85,7 @@ "XTC_sub_sol", "XYZ", "XYZ_psf", "XYZ_bz2", "XYZ_mini", "XYZ_five", # 3 and 5 atoms xyzs for an easy topology - "TXYZ", "ARC", # Tinker files + "TXYZ", "ARC", "ARC_PBC", # Tinker files "PRM", "TRJ", "TRJ_bz2", # Amber (no periodic box) "INPCRD", "PRMpbc", "TRJpbc_bz2", # Amber (periodic box) @@ -294,6 +294,7 @@ XYZ_five = resource_filename(__name__, 'data/five.xyz') TXYZ = resource_filename(__name__, 'data/coordinates/test.txyz') ARC = resource_filename(__name__, 'data/coordinates/test.arc') +ARC_PBC = resource_filename(__name__, 'data/coordinates/new_hexane.arc') PRM = resource_filename(__name__, 'data/Amber/ache.prmtop') TRJ = resource_filename(__name__, 'data/Amber/ache.mdcrd') diff --git a/testsuite/MDAnalysisTests/lib/test_distances.py b/testsuite/MDAnalysisTests/lib/test_distances.py index 37133cb5abc..38f634844e3 100644 --- a/testsuite/MDAnalysisTests/lib/test_distances.py +++ b/testsuite/MDAnalysisTests/lib/test_distances.py @@ -106,7 +106,56 @@ def test_capped_distance_checkbrute(npoints, box, query, method, min_cutoff): assert_equal(np.sort(found_pairs, axis=0), np.sort(indices[1], axis=0)) +@pytest.mark.parametrize('npoints', npoints_1) +@pytest.mark.parametrize('box', boxes_1) +@pytest.mark.parametrize('method', method_1) +@pytest.mark.parametrize('min_cutoff', min_cutoff_1) +def test_self_capped_distance(npoints, box, method, min_cutoff): + np.random.seed(90003) + points = (np.random.uniform(low=0, high=1.0, + size=(npoints, 3))*(boxes_1[0][:3])).astype(np.float32) + max_cutoff = 0.1 + pairs, distance = mda.lib.distances.self_capped_distance(points, + max_cutoff, + min_cutoff=min_cutoff, + box=box, + method=method) + found_pairs, found_distance = [], [] + for i, coord in enumerate(points): + dist = mda.lib.distances.distance_array(coord[None, :], + points[i+1:], + box=box) + if min_cutoff is not None: + idx = np.where((dist < max_cutoff) & (dist > min_cutoff))[0] + else: + idx = np.where((dist < max_cutoff))[0] + for other_idx in idx: + j = other_idx + 1 + i + found_pairs.append((i, j)) + found_distance.append(dist[other_idx]) + assert_equal(len(pairs), len(found_pairs)) + +@pytest.mark.parametrize('box', (None, + np.array([1, 1, 1, 90, 90, 90], dtype=np.float32), + np.array([1, 1, 1, 60, 75, 80], dtype=np.float32))) +@pytest.mark.parametrize('npoints,cutoff,meth', + [(1, 0.02, '_bruteforce_capped_self'), + (1, 0.2, '_bruteforce_capped_self'), + (6000, 0.02, '_pkdtree_capped_self'), + (6000, 0.2, '_pkdtree_capped_self'), + (200000, 0.02, '_pkdtree_capped_self'), + (200000, 0.2, '_bruteforce_capped_self')]) +def test_method_selfselection(box, npoints, cutoff, meth): + np.random.seed(90003) + points = (np.random.uniform(low=0, high=1.0, + size=(npoints, 3))).astype(np.float32) + method = mda.lib.distances._determine_method_self(points, cutoff, box=box) + assert_equal(method.__name__, meth) + +@pytest.mark.parametrize('box', (None, + np.array([1, 1, 1, 90, 90, 90], dtype=np.float32), + np.array([1, 1, 1, 60, 75, 80], dtype=np.float32))) @pytest.mark.parametrize('npoints,cutoff,meth', [(1, 0.02, '_bruteforce_capped'), (1, 0.2, '_bruteforce_capped'), @@ -114,20 +163,11 @@ def test_capped_distance_checkbrute(npoints, box, query, method, min_cutoff): (6000, 0.2, '_pkdtree_capped'), (200000, 0.02, '_pkdtree_capped'), (200000, 0.2, '_bruteforce_capped')]) -def test_method_selection(npoints, cutoff, meth): +def test_method_selection(box, npoints, cutoff, meth): np.random.seed(90003) - box = np.array([1, 1, 1, 90, 90, 90], dtype=np.float32) points = (np.random.uniform(low=0, high=1.0, - size=(npoints, 3)) * (box[:3])).astype(np.float32) - if box is not None: - boxtype = mda.lib.distances._box_check(box) - # Convert [A,B,C,alpha,beta,gamma] to [[A],[B],[C]] - if (boxtype == 'tri_box'): - box = triclinic_vectors(box) - if (boxtype == 'tri_vecs_bad'): - box = triclinic_vectors(triclinic_box(box[0], box[1], box[2])) - method = mda.lib.distances._determine_method(points, points, - cutoff, box=box) + size=(npoints, 3)).astype(np.float32)) + method = mda.lib.distances._determine_method(points, points, cutoff, box=box) assert_equal(method.__name__, meth) diff --git a/testsuite/MDAnalysisTests/topology/test_guessers.py b/testsuite/MDAnalysisTests/topology/test_guessers.py index 574d8bd7a45..7bbb585bf22 100644 --- a/testsuite/MDAnalysisTests/topology/test_guessers.py +++ b/testsuite/MDAnalysisTests/topology/test_guessers.py @@ -25,11 +25,13 @@ from numpy.testing import assert_equal import numpy as np +import MDAnalysis as mda from MDAnalysis.topology import guessers from MDAnalysis.core.topologyattrs import Angles from MDAnalysisTests import make_Universe from MDAnalysisTests.core.test_fragments import make_starshape +from MDAnalysisTests.datafiles import two_water_gro class TestGuessMasses(object): @@ -99,3 +101,12 @@ def test_guess_impropers(): vals = guessers.guess_improper_dihedrals(ag.angles) assert_equal(len(vals), 12) + + +def test_guess_bonds(): + u = mda.Universe(two_water_gro) + bonds = guessers.guess_bonds(u.atoms, u.atoms.positions, u.dimensions) + assert_equal(bonds, ((0, 1), + (0, 2), + (3, 4), + (3, 5))) diff --git a/testsuite/MDAnalysisTests/topology/test_mmtf.py b/testsuite/MDAnalysisTests/topology/test_mmtf.py index 1389e567745..0c4bfd66051 100644 --- a/testsuite/MDAnalysisTests/topology/test_mmtf.py +++ b/testsuite/MDAnalysisTests/topology/test_mmtf.py @@ -108,11 +108,12 @@ def u(self): class TestMMTFFetch(TestMMTFUniverse): @pytest.fixture() - @mock.patch('mmtf.fetch') - def u(self, mock_fetch): + def u(self): top = mmtf.parse(MMTF) - mock_fetch.return_value = top - return mda.fetch_mmtf('173D') # string is irrelevant + with mock.patch('mmtf.fetch') as mock_fetch: + mock_fetch.return_value = top + + return mda.fetch_mmtf('173D') # string is irrelevant class TestSelectModels(object): diff --git a/testsuite/MDAnalysisTests/topology/test_pqr.py b/testsuite/MDAnalysisTests/topology/test_pqr.py index 216a7a71c2d..c6c89ab5d2b 100644 --- a/testsuite/MDAnalysisTests/topology/test_pqr.py +++ b/testsuite/MDAnalysisTests/topology/test_pqr.py @@ -90,7 +90,7 @@ def test_gromacs_flavour(): assert u.atoms[0].type == 'O' assert u.atoms[0].segid == 'SYSTEM' assert not u._topology.types.is_guessed - assert u.atoms[0].radius == pytest.approx(1.48) - assert u.atoms[0].charge == pytest.approx(-0.67) + assert_almost_equal(u.atoms[0].radius, 1.48, decimal=5) + assert_almost_equal(u.atoms[0].charge, -0.67, decimal=5) # coordinatey things assert_almost_equal(u.atoms[0].position, [15.710, 17.670, 23.340], decimal=4) diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index f22e7771770..39384712e51 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -39,6 +39,7 @@ def rotate_universes(): transformed.trajectory.ts.dimensions = np.array([372., 373., 374., 90, 90, 90]) return reference, transformed + def test_rotation_matrix(): # test if the rotation_matrix function is working properly # simple angle and unit vector on origin @@ -60,21 +61,44 @@ def test_rotation_matrix(): matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] assert_array_almost_equal(matrix, ref_matrix, decimal=6) +@pytest.mark.parametrize('point', ( + np.asarray([0, 0, 0]), + np.asarray([[0, 0, 0]])) +) +def test_rotateby_custom_point(rotate_universes, point): + # what happens when we use a custom point for the axis of rotation? + ref_u = rotate_universes[0] + trans_u = rotate_universes[1] + trans = trans_u.trajectory.ts + ref = ref_u.trajectory.ts + vector = [1, 0, 0] + pos = point.reshape(3, ) + angle = 90 + matrix = rotation_matrix(np.deg2rad(angle), vector, pos) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, point=point)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + -def test_rotateby_custom_position(rotate_universes): +@pytest.mark.parametrize('vector', ( + np.asarray([1, 0, 0]), + np.asarray([[1, 0, 0]])) +) +def test_rotateby_vector(rotate_universes, vector): # what happens when we use a custom point for the axis of rotation? ref_u = rotate_universes[0] trans_u = rotate_universes[1] trans = trans_u.trajectory.ts ref = ref_u.trajectory.ts - vector = [1,0,0] - pos = [0,0,0] + point = [0, 0, 0] angle = 90 - matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(angle, vector, point=pos)(trans) + vec = vector.reshape(3, ) + matrix = rotation_matrix(np.deg2rad(angle), vec, point) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, point=point)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) - + + def test_rotateby_atomgroup_cog_nopbc(rotate_universes): # what happens when we rotate arround the center of geometry of a residue # without pbc? @@ -85,12 +109,13 @@ def test_rotateby_atomgroup_cog_nopbc(rotate_universes): center_pos = [6,7,8] vector = [1,0,0] angle = 90 - matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + ref_u.atoms.transform(matrix) selection = trans_u.residues[0].atoms - transformed = rotateby(angle, vector, ag=selection, center='geometry')(trans) + transformed = rotateby(angle, vector, ag=selection, weights=None)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + def test_rotateby_atomgroup_com_nopbc(rotate_universes): # what happens when we rotate arround the center of mass of a residue # without pbc? @@ -102,10 +127,11 @@ def test_rotateby_atomgroup_com_nopbc(rotate_universes): angle = 90 selection = trans_u.residues[0].atoms center_pos = selection.center_of_mass() - matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(angle, vector, ag=selection, center='mass')(trans) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, weights='mass')(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + def test_rotateby_atomgroup_cog_pbc(rotate_universes): # what happens when we rotate arround the center of geometry of a residue @@ -118,11 +144,12 @@ def test_rotateby_atomgroup_cog_pbc(rotate_universes): angle = 90 selection = trans_u.residues[0].atoms center_pos = selection.center_of_geometry(pbc=True) - matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(angle, vector, ag=selection, center='geometry', wrap=True)(trans) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, weights=None, wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + def test_rotateby_atomgroup_com_pbc(rotate_universes): # what happens when we rotate arround the center of mass of a residue # with pbc? @@ -134,12 +161,23 @@ def test_rotateby_atomgroup_com_pbc(rotate_universes): angle = 90 selection = trans_u.residues[0].atoms center_pos = selection.center_of_mass(pbc=True) - matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) - transformed = rotateby(angle, vector, ag=selection, center='mass', wrap=True)(trans) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, weights='mass', wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) - -def test_rotateby_bad_ag(rotate_universes): + + +@pytest.mark.parametrize('ag', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotanag', + 1) +) +def test_rotateby_bad_ag(rotate_universes, ag): # this universe as a box size zero ts = rotate_universes[0].trajectory.ts ag = rotate_universes[0].residues[0].atoms @@ -150,37 +188,80 @@ def test_rotateby_bad_ag(rotate_universes): with pytest.raises(ValueError): rotateby(angle, vector, ag = bad_ag)(ts) -def test_rotateby_bad_position(rotate_universes): + +@pytest.mark.parametrize('point', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotapoint', + 1) +) +def test_rotateby_bad_point(rotate_universes, point): # this universe as a box size zero ts = rotate_universes[0].trajectory.ts # what if the box is in the wrong format? angle = 90 vector = [0, 0, 1] - bad_position = [1] + bad_position = point with pytest.raises(ValueError): rotateby(angle, vector, point=bad_position)(ts) - + + +@pytest.mark.parametrize('direction', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotadirection', + 1) +) +def test_rotateby_bad_direction(rotate_universes, direction): + # this universe as a box size zero + ts = rotate_universes[0].trajectory.ts + # what if the box is in the wrong format? + angle = 90 + point = [0, 0, 0] + with pytest.raises(ValueError): + rotateby(angle, direction, point=point)(ts) + + def test_rotateby_bad_pbc(rotate_universes): # this universe as a box size zero ts = rotate_universes[0].trajectory.ts ag = rotate_universes[0].residues[0].atoms # is pbc passed to the center methods? # if yes it should raise an exception for boxes that are zero in size + vector = [1, 0, 0] angle = 90 - vector = [0, 0, 1] with pytest.raises(ValueError): rotateby(angle, vector, ag = ag, wrap=True)(ts) -def test_rotateby_bad_center(rotate_universes): + +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_rotateby_bad_weights(rotate_universes, weights): # this universe as a box size zero ts = rotate_universes[0].trajectory.ts ag = rotate_universes[0].residues[0].atoms # what if a wrong center type name is passed? angle = 90 vector = [0, 0, 1] - bad_center = " " - with pytest.raises(ValueError): - rotateby(angle, vector, ag = ag, center=bad_center)(ts) + bad_weights = " " + with pytest.raises(TypeError): + rotateby(angle, vector, ag = ag, weights=bad_weights)(ts) + def test_rotateby_no_masses(rotate_universes): # this universe as a box size zero @@ -190,8 +271,9 @@ def test_rotateby_no_masses(rotate_universes): angle = 90 vector = [0, 0, 1] bad_center = "mass" - with pytest.raises(AttributeError): - rotateby(angle, vector, ag = ag, center=bad_center)(ts) + with pytest.raises(TypeError): + rotateby(angle, vector, ag = ag, weights=bad_center)(ts) + def test_rotateby_no_args(rotate_universes): # this universe as a box size zero diff --git a/testsuite/MDAnalysisTests/transformations/test_wrap.py b/testsuite/MDAnalysisTests/transformations/test_wrap.py new file mode 100644 index 00000000000..09cdfe83868 --- /dev/null +++ b/testsuite/MDAnalysisTests/transformations/test_wrap.py @@ -0,0 +1,129 @@ +#-*- 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. +# +# 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 + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +import MDAnalysis as mda +from MDAnalysis.transformations import wrap, unwrap +from MDAnalysisTests.datafiles import fullerene + + +@pytest.fixture() +def wrap_universes(): + # create the Universe objects for the tests + # this universe is used for the unwrap testing cases + reference = mda.Universe(fullerene) + transformed = mda.Universe(fullerene) + transformed.dimensions = np.asarray([10, 10, 10, 90, 90, 90], np.float32) + transformed.atoms.wrap() + + return reference, transformed + + +@pytest.mark.parametrize('ag', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotanag', + 1) +) +def test_wrap_bad_ag(wrap_universes, ag): + # this universe has a box size zero + ts = wrap_universes[0].trajectory.ts + # what happens if something other than an AtomGroup is given? + bad_ag = ag + with pytest.raises(AttributeError): + wrap(bad_ag)(ts) + + +def test_wrap_no_options(wrap_universes): + # since were testing if the wrapping works + # the reference and the transformed are switched + trans, ref = wrap_universes + trans.dimensions = ref.dimensions + wrap(trans.atoms)(trans.trajectory.ts) + assert_array_almost_equal(trans.trajectory.ts.positions, ref.trajectory.ts.positions, decimal=6) + + +@pytest.mark.parametrize('compound', ( + "group", + "residues", + "segments", + "fragments") +) +def test_wrap_with_compounds(wrap_universes, compound): + trans, ref = wrap_universes + trans.dimensions = ref.dimensions + # reference is broken so let's rebuild it + unwrap(ref.atoms)(ref.trajectory.ts) + # lets shift the transformed universe atoms 2*boxlength before wrapping + # the compound keywords will shift the whole molecule back into the unit cell + trans.atoms.positions += np.float32([20,0,0]) + wrap(trans.atoms, compound=compound)(trans.trajectory.ts) + assert_array_almost_equal(trans.trajectory.ts.positions, ref.trajectory.ts.positions, decimal=6) + + +def test_wrap_api(wrap_universes): + trans, ref = wrap_universes + trans.dimensions = ref.dimensions + trans.trajectory.add_transformations(wrap(trans.atoms)) + assert_array_almost_equal(trans.trajectory.ts.positions, ref.trajectory.ts.positions, decimal=6) + + +@pytest.mark.parametrize('ag', ( + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + np.array([0, 1, 2, 3, 4]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]), + np.array([[0], [1], [2]]), + 'thisisnotanag', + 1) +) +def test_unwrap_bad_ag(wrap_universes, ag): + # this universe has a box size zero + ts = wrap_universes[0].trajectory.ts + # what happens if something other than an AtomGroup is given? + bad_ag = ag + with pytest.raises(AttributeError): + unwrap(bad_ag)(ts) + + +def test_unwrap(wrap_universes): + ref, trans = wrap_universes + # after rebuild the trans molecule it should match the reference + unwrap(trans.atoms)(trans.trajectory.ts) + assert_array_almost_equal(trans.trajectory.ts.positions, ref.trajectory.ts.positions, decimal=6) + + +def test_unwrap_api(wrap_universes): + ref, trans = wrap_universes + # after rebuild the trans molecule it should match the reference + trans.trajectory.add_transformations(unwrap(trans.atoms)) + assert_array_almost_equal(trans.trajectory.ts.positions, ref.trajectory.ts.positions, decimal=6)