diff --git a/package/CHANGELOG b/package/CHANGELOG index 57adc0f7b81..e1eecbb330f 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -95,7 +95,7 @@ Enhancements * Added lib.distances.capped_distance function to quickly calculate all distances up to a given maximum distance (PR #1941) * Added a rotation coordinate transformation (PR #1937) - * Added a box centering trajectory transformation (PR #1946) + * Added box, plane and axis centering trajectory transformations (PR #1946 and #1973) * Added a on-the-fly trajectory transformations API and a coordinate translation function (Issue #786) * Added various duecredit stubs diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index 0e0794fbdb6..1c797c4a5b0 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -93,5 +93,5 @@ def wrapped(ts): from __future__ import absolute_import -from .translate import translate, center_in_box +from .translate import translate, center_in_box, center_in_plane, center_in_axis from .rotate import rotateby diff --git a/package/MDAnalysis/transformations/translate.py b/package/MDAnalysis/transformations/translate.py index bdcc16870f0..ef0b81fc938 100644 --- a/package/MDAnalysis/transformations/translate.py +++ b/package/MDAnalysis/transformations/translate.py @@ -25,15 +25,19 @@ Trajectory translation --- :mod:`MDAnalysis.transformations.translate` ====================================================================== -Translate the coordinates of a given trajectory by a given vector. -The vector can either be user defined, using the function :func:`translate` -or defined by centering an AtomGroup in the unit cell using the function -:func:`center_in_box` +Translate the coordinates of a given trajectory. Coordinates may be +translated by a vector, using the function :func:`translate` or defined +by centering an AtomGroup in the unit cell, a plane or axis using the functions +:func:`center_in_box`, :func:`center_in_plane` or :func:`center_in_axis`, +respectively. .. autofunction:: translate .. autofunction:: center_in_box +.. autofunction:: center_in_plane + +.. autofunction:: center_in_axis """ from __future__ import absolute_import, division @@ -42,6 +46,8 @@ from functools import partial from ..lib.mdamath import triclinic_vectors +from ..lib.util import get_weights +from ..lib._cutil import make_whole def translate(vector): """ @@ -50,8 +56,14 @@ def translate(vector): Example ------- - ts = MDAnalysis.transformations.translate([1,2,3])(ts) + Translate the coordinates of the system by the [1, 2, 3] vector: + + .. code-block:: python + + transform = mda.transformations.translate([1,2,3]) + u.trajectory.add_transformations(transform) + Parameters ---------- vector: array-like @@ -59,7 +71,7 @@ def translate(vector): Returns ------- - :class:`~MDAnalysis.coordinates.base.Timestep` object + MDAnalysis.coordinates.base.Timestep """ if len(vector)>2: @@ -75,7 +87,7 @@ def wrapped(ts): return wrapped -def center_in_box(ag, center='geometry', point=None, wrap=False): +def center_in_box(ag, weights=None, center_to=None, wrap=False, unwrap=False): """ Translates the coordinates of a given :class:`~MDAnalysis.coordinates.base.Timestep` instance so that the center of geometry/mass of the given :class:`~MDAnalysis.core.groups.AtomGroup` @@ -84,58 +96,74 @@ def center_in_box(ag, center='geometry', point=None, wrap=False): Example ------- + + Translate the center of mass of of the second residue of the universe u to the center of the unit + cell: .. code-block:: python ag = u.residues[1].atoms - ts = MDAnalysis.transformations.center(ag,center='mass')(ts) + transform = mda.transformations.center(ag, weights='mass') + u.trajectory.add_transformations(transform) Parameters ---------- ag: AtomGroup atom group to be centered on the unit cell. - center: str, optional - used to choose the method of centering on the given atom group. Can be 'geometry' - or 'mass' - point: array-like, optional + 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. + center_to: array-like, optional overrides the unit cell center - the coordinates of the Timestep are translated so that the center of mass/geometry of the given AtomGroup is aligned to this position instead. Defined as an array of size 3. 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. - + before calculating the weighted center. Default is `False`, no changes to the atom + coordinates are done before calculating the center of the AtomGroup. + unwrap: bool, optional + If `True`, all the atoms from the given AtomGroup will be moved so as to not break + any bonds over periodic boundaries before calculating the weighted center. Default + is `False`, no changes to the atom coordinates are done before calculating the center + of the AtomGroup. + Returns ------- - :class:`~MDAnalysis.coordinates.base.Timestep` object - + MDAnalysis.coordinates.base.Timestep + """ pbc_arg = wrap - if point: - point = np.asarray(point, np.float32) - if point.shape != (3, ) and point.shape != (1, 3): - raise ValueError('{} is not a valid point'.format(point)) + if center_to: + center_to = np.asarray(center_to, np.float32) + if center_to.shape != (3, ) and center_to.shape != (1, 3): + raise ValueError('{} is not a valid point'.format(center_to)) + else: + center_to = center_to.reshape(3, ) 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.") + if unwrap and wrap: + raise ValueError("wrap and unwrap can't be both True") + center_method = partial(atoms.center, weights, pbc=wrap) def wrapped(ts): - if point is None: + if center_to is None: boxcenter = np.sum(ts.triclinic_dimensions, axis=0) / 2 else: - boxcenter = point - + boxcenter = center_to + if unwrap: + make_whole(ag) ag_center = center_method() vector = boxcenter - ag_center @@ -145,3 +173,195 @@ def wrapped(ts): return wrapped + +def center_in_plane(ag, plane, center_to="center", weights=None, wrap=False, unwrap=False): + """ + Translates the coordinates of a given :class:`~MDAnalysis.coordinates.base.Timestep` + instance so that the center of geometry/mass of the given :class:`~MDAnalysis.core.groups.AtomGroup` + is centered on a plane. This plane can be the yz, xz or xy planes. + + Example + ------- + + Translate the center of mass of the second residue of the universe u to the center of the + xy plane: + + .. code-block:: python + + ag = u.residues[1].atoms + transform = mda.transformations.center_in_plane(ag, xy, weights='mass') + u.trajectory.add_transformations(transform) + + Parameters + ---------- + ag: AtomGroup + atom group to be centered on the unit cell. + plane: str + used to define the plane on which the given AtomGroup will be centered. Suported + planes are yz, xz and xy planes. + center_to: array-like or str + coordinates to which the axes are centered. Can be an array of three coordinate + values or `center` which centers the AtomGroup coordinates on the center of the + unit cell. Default is `center`. + 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 weighted center. Default is `False`, no changes to the atom + coordinates are done before calculating the center of the AtomGroup. + unwrap: bool, optional + If `True`, all the atoms from the given AtomGroup will be moved so as to not break + any bonds over periodic boundaries before calculating the weighted center. Default + is `False`, no changes to the atom coordinates are done before calculating the center + of the AtomGroup. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + + """ + try: + axes = {'yz' : 0, 'xz' : 1, 'xy' : 2} + plane = axes[plane] + except (KeyError, TypeError): + raise ValueError('{} is not a valid plane'.format(plane)) + try: + if center_to != "center": + center_to = np.asarray(center_to, np.float32) + if center_to.shape != (3, ) and center_to.shape != (1, 3): + raise ValueError('{} is not a valid "center_to"'.format(center_to)) + center_to = center_to.reshape(3, ) + except ValueError: + raise ValueError('{} is not a valid "center_to"'.format(center_to)) + try: + atoms = ag.atoms + except AttributeError: + 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.") + if unwrap and wrap: + raise ValueError("wrap and unwrap can't be both True") + center_method = partial(atoms.center, weights, pbc=wrap) + + def wrapped(ts): + boxcenter = ts.triclinic_dimensions.sum(axis=0) / 2.0 + if center_to == 'center': + _origin = boxcenter + else: + _origin = center_to + if unwrap: + make_whole(ag) + position = center_method() + position[plane] = _origin[plane] + vector = np.asarray(position, np.float32) - center_method() + ts.positions += vector + + return ts + + return wrapped + + +def center_in_axis(ag, axis, center_to="center", weights=None, wrap=False, unwrap=False): + """ + Translates the coordinates of a given :class:`~MDAnalysis.coordinates.base.Timestep` + instance so that the center of geometry/mass of the given :class:`~MDAnalysis.core.groups.AtomGroup` + is centered on an axis. This axis can only be parallel to the x, y or z planes. + + Example + ------- + + Translate the center of mass of of the second residue of the universe u to the center of the line + parallel to the x axis that passes through the [0, 0, 1] point: + + .. code-block:: python + + ag = u.residues[1].atoms + transform = mda.transformations.center_in_axis(ag, 'x', [0,0,1], + weights='mass') + u.trajectory.add_transformation(transform) + + Parameters + ---------- + ag: AtomGroup + atom group to be centered on the unit cell. + axis: str + used to define the plane on which the given AtomGroup will be centered. Defined as + a string or an array-like of the plane and the value of its translation. Suported + planes are x, y and z. + center_to: array-like or str + coordinate on which the axes are centered. Can be an array of three coordinates or + `center` to shift the origin to the center of the unit cell. Default is `None` + which sets [0, 0, 0] as the origin of the axis. + 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 weighted center. Default is `False`, no changes to the atom + coordinates are done before calculating the center of the AtomGroup. + unwrap: bool, optional + If `True`, all the atoms from the given AtomGroup will be moved so as to not break + any bonds over periodic boundaries before calculating the weighted center. Default + is `False`, no changes to the atom coordinates are done before calculating the center + of the AtomGroup. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + + """ + try: + axes = {'x' : 0, 'y': 1, 'z' : 2} + axis = axes[axis] + except (KeyError, TypeError): + raise ValueError('{} is not a valid axis'.format(axis)) + try: + if center_to != "center": + center_to = np.asarray(center_to, np.float32) + if center_to.shape != (3, ) and center_to.shape != (1, 3): + raise ValueError('{} is not a valid "center_to"'.format(center_to)) + center_to = center_to.reshape(3, ) + except ValueError: + raise ValueError('{} is not a valid "center_to"'.format(center_to)) + try: + atoms = ag.atoms + except AttributeError: + 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.") + if unwrap and wrap: + raise ValueError("wrap and unwrap can't be both True") + center_method = partial(ag.atoms.center, weights, pbc=wrap) + + def wrapped(ts): + if center_to == 'center': + _origin = ts.triclinic_dimensions.sum(axis=0) / 2.0 + else: + _origin = center_to + if unwrap: + make_whole(ag) + ag_center = center_method() + center = _origin + center[axis] = ag_center[axis] + vector = center - ag_center + ts.positions += vector + + return ts + + return wrapped + diff --git a/testsuite/MDAnalysisTests/transformations/test_translate.py b/testsuite/MDAnalysisTests/transformations/test_translate.py index 413c3338bce..250ef82c30e 100644 --- a/testsuite/MDAnalysisTests/transformations/test_translate.py +++ b/testsuite/MDAnalysisTests/transformations/test_translate.py @@ -28,8 +28,9 @@ from numpy.testing import assert_array_almost_equal import MDAnalysis as mda -from MDAnalysis.transformations import translate, center_in_box +from MDAnalysis.transformations import translate, center_in_box, center_in_plane, center_in_axis from MDAnalysisTests import make_Universe +from MDAnalysisTests.datafiles import fullerene @pytest.fixture() @@ -43,6 +44,18 @@ def translate_universes(): return reference, transformed +@pytest.fixture() +def translate_unwrap_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 + + def test_translate_coords(translate_universes): ref_u, trans_u = translate_universes ref = ref_u.trajectory.ts @@ -66,7 +79,15 @@ def test_translate_vector(translate_universes, vector): with pytest.raises(ValueError): translate(vector)(ts) - + +def test_translate_coords_dtype(translate_universes): + # is the dtype of the coordinates correct? + trans_u = translate_universes[1] + vector = np.float32([1, 2, 3]) + trans = translate(vector)(trans_u.trajectory.ts) + assert(trans.positions.dtype == np.float32) + + def test_translate_transformations_api(translate_universes): # test if the translate transformation works when using the # on-the-fly transformations API @@ -78,49 +99,88 @@ def test_translate_transformations_api(translate_universes): assert_array_almost_equal(trans_u.trajectory.ts.positions, ref.positions, decimal=6) -def test_center_in_box_bad_ag(translate_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_center_in_box_bad_ag(translate_universes, ag): # this universe has a box size zero ts = translate_universes[0].trajectory.ts # what happens if something other than an AtomGroup is given? - bad_ag = 1 + bad_ag = ag with pytest.raises(ValueError): center_in_box(bad_ag)(ts) -@pytest.mark.parametrize('point', ( +@pytest.mark.parametrize('center_to', ( [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]])) + np.array([[0], [1], [2]]), + 'thisisnotacenter', + 1) ) -def test_center_in_box_bad_point(translate_universes, point): +def test_center_in_box_bad_center_to(translate_universes, center_to): ts = translate_universes[0].trajectory.ts ag = translate_universes[0].residues[0].atoms + bad_center_to = center_to # what if the box is in the wrong format? with pytest.raises(ValueError): - center_in_box(ag, point=point)(ts) + center_in_box(ag, center_to=bad_center_to)(ts) -def test_center_in_box_bad_pbc(translate_universes): +def test_center_in_box_bad_wrap(translate_universes): # this universe has a box size zero ts = translate_universes[0].trajectory.ts ag = translate_universes[0].residues[0].atoms - # is pbc passed to the center methods? + # is wrap passed to the center methods? # if yes it should raise an exception for boxes that are zero in size with pytest.raises(ValueError): center_in_box(ag, wrap=True)(ts) -def test_center_in_box_bad_center(translate_universes): +def test_center_in_box_bad_unwrap(translate_universes): # this universe has a box size zero ts = translate_universes[0].trajectory.ts ag = translate_universes[0].residues[0].atoms - # what if a wrong center type name is passed? - bad_center = " " + # is wrap passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size with pytest.raises(ValueError): - center_in_box(ag, center=bad_center)(ts) + center_in_box(ag, unwrap=True)(ts) + + +def test_center_in_box_bad_wrap_unwrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[1].trajectory.ts + ag = translate_universes[1].residues[0].atoms + # the two options are not compatible so it should throw an error + with pytest.raises(ValueError): + center_in_box(ag, wrap=True, unwrap=True)(ts) + + +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + np.array([0, 1]), + np.array([0, 1, 2, 3]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_center_in_box_bad_weights(translate_universes, weights): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if a wrong center type name is passed? + with pytest.raises(TypeError): + center_in_box(ag, weights=weights)(ts) def test_center_in_box_no_masses(translate_universes): @@ -129,10 +189,18 @@ def test_center_in_box_no_masses(translate_universes): ag = translate_universes[0].residues[0].atoms # if the universe has no masses and `mass` is passed as the center arg bad_center = "mass" - with pytest.raises(AttributeError): - center_in_box(ag, center=bad_center)(ts) + with pytest.raises(TypeError): + center_in_box(ag, center_of=bad_center)(ts) +def test_center_in_box_coords_dtype(translate_universes): + # is the dtype of the coordinates correct? + trans_u = translate_universes[1] + ag = trans_u.residues[0].atoms + trans = center_in_box(ag)(trans_u.trajectory.ts) + assert(trans.positions.dtype == np.float32) + + def test_center_in_box_coords_no_options(translate_universes): # what happens when we center the coordinates arround the center of geometry of a residue? ref_u, trans_u = translate_universes @@ -145,7 +213,7 @@ def test_center_in_box_coords_no_options(translate_universes): assert_array_almost_equal(trans.positions, ref.positions, decimal=6) -def test_center_in_box_coords_with_pbc(translate_universes): +def test_center_in_box_coords_with_wrap(translate_universes): # what happens when we center the coordinates arround the center of geometry of a residue? # using pbc into account for center of geometry calculation ref_u, trans_u = translate_universes @@ -159,15 +227,28 @@ def test_center_in_box_coords_with_pbc(translate_universes): assert_array_almost_equal(trans.positions, ref.positions, decimal=6) +def test_center_in_box_coords_with_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + ag = trans_u.atoms + box_center = np.float32([5, 5, 5]) + ref_center = ref_u.atoms.center(weights=None) + ref.positions += box_center - ref_center + trans = center_in_box(ag, unwrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + def test_center_in_box_coords_with_mass(translate_universes): # using masses for calculating the center of the atomgroup ref_u, trans_u = translate_universes ref = ref_u.trajectory.ts ag = trans_u.residues[24].atoms box_center = np.float32([186., 186.5, 187.]) - ref_center = ag.center_of_mass() + ref_center = ag.center(weights=ag.masses) ref.positions += box_center - ref_center - trans = center_in_box(ag, center="mass")(trans_u.trajectory.ts) + trans = center_in_box(ag, weights="mass")(trans_u.trajectory.ts) assert_array_almost_equal(trans.positions, ref.positions, decimal=6) @@ -180,11 +261,11 @@ def test_center_in_box_coords_with_box(translate_universes): box_center = np.float32(newpoint) ref_center = np.float32([6, 7, 8]) ref.positions += box_center - ref_center - trans = center_in_box(ag, point=newpoint)(trans_u.trajectory.ts) + trans = center_in_box(ag, center_to=newpoint)(trans_u.trajectory.ts) assert_array_almost_equal(trans.positions, ref.positions, decimal=6) -def test_center_in_box_coords_all_options(translate_universes): +def test_center_in_box_coords_all_options_wrap(translate_universes): # what happens when we center the coordinates arround the center of geometry of a residue? # using pbc into account for center of geometry calculation ref_u, trans_u = translate_universes @@ -192,13 +273,481 @@ def test_center_in_box_coords_all_options(translate_universes): ag = trans_u.residues[24].atoms newpoint = [1000, 1000, 1000] box_center = np.float32(newpoint) - ref_center = ag.center_of_mass(pbc=True) + ref_center = ag.center(weights=ag.masses, pbc=True) ref.positions += box_center - ref_center - trans = center_in_box(ag, center='mass', wrap=True, point=newpoint)(trans_u.trajectory.ts) + trans = center_in_box(ag, weights='mass', wrap=True, center_to=newpoint)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_box_coords_all_options_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + ag = trans_u.atoms + newpoint = [1000, 1000, 1000] + box_center = np.float32(newpoint) + ref_center = ref_u.atoms.center(weights=ref_u.atoms.masses) + ref.positions += box_center - ref_center + trans = center_in_box(ag, weights='mass', unwrap=True, center_to=newpoint)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.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_center_in_axis_bad_ag(translate_universes, ag): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + # what happens if something other than an AtomGroup is given? + bad_ag = ag + axis = 'x' + with pytest.raises(ValueError): + center_in_axis(bad_ag, axis)(ts) + + +@pytest.mark.parametrize('center_to', ( + [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]]), + 'thisisnotacenter', + 1) +) +def test_center_in_axis_bad_center_to(translate_universes, center_to): + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if the box is in the wrong format? + axis = 'x' + with pytest.raises(ValueError): + center_in_axis(ag, axis, center_to=center_to)(ts) + + +@pytest.mark.parametrize('axis', ( + 1, + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + "xyz", + "notanaxis") +) +def test_center_in_axis_bad_axis(translate_universes, axis): + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if the box is in the wrong format? + bad_axis = axis + with pytest.raises(ValueError): + center_in_axis(ag, bad_axis)(ts) + + +def test_center_in_axis_bad_wrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # is wrap passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size + axis = 'x' + with pytest.raises(ValueError): + center_in_axis(ag,axis, wrap=True)(ts) + + +def test_center_in_axis_bad_unwrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # is unwrap passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size + axis = 'x' + with pytest.raises(ValueError): + center_in_axis(ag,axis, unwrap=True)(ts) + + +def test_center_in_axis_bad_wrap_unwrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[1].trajectory.ts + ag = translate_universes[1].residues[0].atoms + # the two options are not compatible so it should throw an error + axis = 'x' + with pytest.raises(ValueError): + center_in_axis(ag,axis, wrap=True, unwrap=True)(ts) + + +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + np.array([0, 1]), + np.array([0, 1, 2, 3]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_center_in_axis_bad_weights(translate_universes, weights): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if a wrong center type name is passed? + axis = 'x' + with pytest.raises(TypeError): + center_in_axis(ag, axis, weights=weights)(ts) + + +def test_center_in_axis_no_masses(translate_universes): + # this universe has no masses + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # if the universe has no masses and `mass` is passed as the center arg + bad_weights = "mass" + axis = 'x' + with pytest.raises(TypeError): + center_in_axis(ag, axis, weights=bad_weights)(ts) + +def test_center_in_axis_coords_dtype(translate_universes): + # is the dtype of the coordinates correct? + trans_u = translate_universes[1] + axis = 'x' + ag = trans_u.residues[0].atoms + trans = center_in_axis(ag, axis)(trans_u.trajectory.ts) + assert(trans.positions.dtype == np.float32) + + +def test_center_in_axis_coords_no_options(translate_universes): + # what happens when we center the coordinates on an axis without any other options? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[0].atoms + box_center = np.float32([186., 186.5, 187.]) + center_to = box_center + ref_center = np.float32([6, 7, 8]) + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_axis_coords_with_wrap(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue + # taking pbc into account for center of geometry calculation? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + trans_u.dimensions = [363., 364., 365., 90., 90., 90.] + ag = trans_u.residues[24].atoms + box_center = np.float32([181.5, 182., 182.5]) + center_to = box_center + ref_center = ag.center_of_geometry(pbc=True) + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, wrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_axis_coords_with_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue + # taking pbc into account for center of geometry calculation? + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + ag = trans_u.atoms + box_center = np.float32([5, 5, 5]) + center_to = box_center + ref_center = ref_u.atoms.center(weights=None) + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, unwrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_axis_coords_with_mass(translate_universes): + # using masses for calculating the center of the atomgroup + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[24].atoms + ref_center = ag.center(weights=ag.masses) + box_center = np.float32([186., 186.5, 187.]) + center_to = box_center + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, weights="mass")(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_axis_coords_with_coord_center_to(translate_universes): + # what happens when we use a custom coordinate as center_to point? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[0].atoms + center_to = [1000, 1000, 1000] + ref_center = np.float32([6, 7, 8]) + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, center_to=center_to)(trans_u.trajectory.ts) assert_array_almost_equal(trans.positions, ref.positions, decimal=6) -def test_center_transformations_api(translate_universes): +def test_center_in_axis_coords_all_options_wrap(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[24].atoms + neworigin = [1000, 1000, 1000] + ref_center = ag.center(weights=ag.masses, pbc=True) + axis = 'x' + axis_point = np.float32([ref_center[0], neworigin[1], neworigin[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, weights='mass', wrap=True, center_to=neworigin)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_axis_coords_all_options_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + ag = trans_u.atoms + neworigin = [1000, 1000, 1000] + ref_center = ref_u.atoms.center(weights=ref_u.atoms.masses) + axis = 'x' + axis_point = np.float32([ref_center[0], neworigin[1], neworigin[2]]) + ref.positions += axis_point - ref_center + trans = center_in_axis(ag, axis, weights='mass', unwrap=True, center_to=neworigin)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.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_center_in_plane_bad_ag(translate_universes, ag): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + # what happens if something other than an AtomGroup is given? + bad_ag = ag + plane = 'yz' + with pytest.raises(ValueError): + center_in_plane(bad_ag, plane)(ts) + + +@pytest.mark.parametrize('center_to', ( + [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]]), + 'thisisnotacenter', + 1) +) +def test_center_in_plane_bad_center_to(translate_universes, center_to): + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if the box is in the wrong format? + plane = 'yz' + with pytest.raises(ValueError): + center_in_plane(ag, plane, center_to=center_to)(ts) + + +@pytest.mark.parametrize('plane', ( + 1, + [0, 1], + [0, 1, 2, 3, 4], + np.array([0, 1]), + "xyz", + "notaplane") +) +def test_center_in_plane_bad_plane(translate_universes, plane): + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if the box is in the wrong format? + bad_plane = plane + with pytest.raises(ValueError): + center_in_plane(ag, bad_plane)(ts) + + +def test_center_in_plane_bad_wrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # is wrap passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size + plane = 'yz' + with pytest.raises(ValueError): + center_in_plane(ag, plane, wrap=True)(ts) + + +def test_center_in_plane_bad_unwrap(translate_universes): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # is unwrap passed to the center methods? + # if yes it should raise an exception for boxes that are zero in size + plane = 'yz' + with pytest.raises(ValueError): + center_in_plane(ag, plane, unwrap=True)(ts) + + +@pytest.mark.parametrize('weights', ( + " ", + "totallynotmasses", + 123456789, + np.array([0, 1]), + np.array([0, 1, 2, 3]), + np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])) +) +def test_center_in_plane_bad_weights(translate_universes, weights): + # this universe has a box size zero + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # what if a wrong center type name is passed? + plane = 'yz' + with pytest.raises(TypeError): + center_in_plane(ag, plane, weights=weights)(ts) + + +def test_center_in_planes_no_masses(translate_universes): + # this universe has no masses + ts = translate_universes[0].trajectory.ts + ag = translate_universes[0].residues[0].atoms + # if the universe has no masses and `mass` is passed as the center arg + weights = "mass" + plane = 'yz' + with pytest.raises(TypeError): + center_in_plane(ag, plane, weights=weights)(ts) + +def test_center_in_plane_coords_dtype(translate_universes): + # is the dtype of the coordinates correct? + trans_u = translate_universes[1] + plane = 'yz' + ag = trans_u.residues[0].atoms + trans = center_in_plane(ag, plane)(trans_u.trajectory.ts) + assert(trans.positions.dtype == np.float32) + + +def test_center_in_plane_coords_no_options(translate_universes): + # what happens when we center the coordinates on a plane without special options? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[0].atoms + box_center = np.float32([186., 186.5, 187.]) + ref_center = np.float32([6, 7, 8]) + plane = 'yz' + plane_point = np.asarray([box_center[0], ref_center[1], ref_center[2]], np.float32) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_with_wrap(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + trans_u.dimensions = [363., 364., 365., 90., 90., 90.] + box_center = np.float32([181.5, 182., 182.5]) + ag = trans_u.residues[24].atoms + ref_center = ag.center(weights=None, pbc=True) + plane = 'yz' + plane_point = np.asarray([box_center[0], ref_center[1], ref_center[2]], np.float32) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, wrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_with_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + box_center = np.float32([5, 5, 5]) + ag = trans_u.atoms + ref_center = ref_u.atoms.center(weights=None) + plane = 'yz' + plane_point = np.asarray([box_center[0], ref_center[1], ref_center[2]], np.float32) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, unwrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_with_mass(translate_universes): + # using masses for calculating the center of the atomgroup + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[24].atoms + box_center = np.float32([186., 186.5, 187.]) + ref_center = ag.center(weights=ag.masses) + plane = 'yz' + plane_point = np.asarray([box_center[0], ref_center[1], ref_center[2]], np.float32) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, weights="mass")(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_with_coord_center_to(translate_universes): + # what happens when we use a custom coordinate as center_to point? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[0].atoms + center_to = [1000, 1000, 1000] + box_center = np.float32([186., 186.5, 187.]) + ref_center = np.float32([6, 7, 8]) + plane = 'yz' + plane_point = np.float32([center_to[0], ref_center[1], ref_center[2]]) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, center_to=center_to)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_all_options_wrap(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ag = trans_u.residues[24].atoms + trans_u.dimensions = [363., 364., 365., 90., 90., 90.] + box_center = np.float32([181.5, 182., 182.5]) + center_to = [1000, 1000, 1000] + ref_center = ag.center(weights=ag.masses, pbc=True) + plane = 'yz' + plane_point = np.float32([center_to[0], ref_center[1], ref_center[2]]) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, center_to=center_to, weights='mass', wrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_plane_coords_all_options_unwrap(translate_unwrap_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + # using pbc into account for center of geometry calculation + ref_u, trans_u = translate_unwrap_universes + ref = ref_u.trajectory.ts + ag = trans_u.atoms + box_center = np.float32([5, 5, 5]) + center_to = [1000, 1000, 1000] + ref_center = ref_u.atoms.center(weights=ref_u.atoms.masses) + plane = 'yz' + plane_point = np.float32([center_to[0], ref_center[1], ref_center[2]]) + ref.positions += plane_point - ref_center + trans = center_in_plane(ag, plane, center_to=center_to, weights='mass', unwrap=True)(trans_u.trajectory.ts) + assert_array_almost_equal(trans.positions, ref.positions, decimal=6) + + +def test_center_in_box_transformations_api(translate_universes): # test if the translate transformation works when using the # on-the-fly transformations API ref_u, trans_u = translate_universes @@ -209,3 +758,32 @@ def test_center_transformations_api(translate_universes): ag = trans_u.residues[0].atoms trans_u.trajectory.add_transformations(center_in_box(ag)) assert_array_almost_equal(trans_u.trajectory.ts.positions, ref.positions, decimal=6) + + +def test_center_in_axis_transformations_api(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ref_center = np.float32([6, 7, 8]) + box_center = np.float32([186., 186.5, 187.]) + center_to = box_center + axis = 'x' + axis_point = np.float32([ref_center[0], center_to[1], center_to[2]]) + ref.positions += axis_point - ref_center + ag = trans_u.residues[0].atoms + trans_u.trajectory.add_transformations(center_in_axis(ag, axis)) + assert_array_almost_equal(trans_u.trajectory.ts.positions, ref.positions, decimal=6) + + +def test_center_in_plane_transformations_api(translate_universes): + # what happens when we center the coordinates arround the center of geometry of a residue? + ref_u, trans_u = translate_universes + ref = ref_u.trajectory.ts + ref_center = np.float32([6, 7, 8]) + box_center = np.float32([186., 186.5, 187.]) + plane = 'yz' + plane_point = np.asarray([box_center[0], ref_center[1], ref_center[2]], np.float32) + ref.positions += plane_point - ref_center + ag = trans_u.residues[0].atoms + trans_u.trajectory.add_transformations(center_in_plane(ag, plane)) + assert_array_almost_equal(trans_u.trajectory.ts.positions, ref.positions, decimal=6)