From 000f1e765c77b1cf68db999055fe1c19b067e3a2 Mon Sep 17 00:00:00 2001 From: davidercruz Date: Sat, 21 Jul 2018 02:43:38 +0100 Subject: [PATCH 01/19] fixed issue with rotation from non-origin point calculation --- package/MDAnalysis/transformations/rotate.py | 5 +++- .../transformations/test_rotate.py | 30 ++++++++++++------- 2 files changed, 24 insertions(+), 11 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index a9712d2a44a..3de862dcf2d 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -119,8 +119,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] + translation = matrix[:3, 3] ts.positions= np.dot(ts.positions, rotation) + ts.positions += translation return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index f22e7771770..9422b4787ca 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -70,8 +70,10 @@ def test_rotateby_custom_position(rotate_universes): vector = [1,0,0] pos = [0,0,0] angle = 90 - matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] - ref.positions = np.dot(ref.positions, matrix) + matrix = rotation_matrix(np.deg2rad(angle), vector, pos) + rotation = matrix[:3, :3] + translation = matrix[:3, 3] + ref.positions = np.dot(ref.positions, rotation) + translation transformed = rotateby(angle, vector, point=pos)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) @@ -85,8 +87,10 @@ 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) + rotation = matrix[:3, :3] + translation = matrix[:3, 3] + ref.positions = np.dot(ref.positions, rotation) + translation selection = trans_u.residues[0].atoms transformed = rotateby(angle, vector, ag=selection, center='geometry')(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) @@ -102,8 +106,10 @@ 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) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + rotation = matrix[:3, :3] + translation = matrix[:3, 3] + ref.positions = np.dot(ref.positions, rotation) + translation transformed = rotateby(angle, vector, ag=selection, center='mass')(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) @@ -118,8 +124,10 @@ 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) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + rotation = matrix[:3, :3] + translation = matrix[:3, 3] + ref.positions = np.dot(ref.positions, rotation) + translation transformed = rotateby(angle, vector, ag=selection, center='geometry', wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) @@ -134,8 +142,10 @@ 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) + matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) + rotation = matrix[:3, :3] + translation = matrix[:3, 3] + ref.positions = np.dot(ref.positions, rotation) + translation transformed = rotateby(angle, vector, ag=selection, center='mass', wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) From 4f23420cd12741795081bf7693be94aca1b918a3 Mon Sep 17 00:00:00 2001 From: davidercruz Date: Sun, 22 Jul 2018 00:29:20 +0100 Subject: [PATCH 02/19] rotaateby now behaves as rotate, updated tests to use transform(matrix) --- package/MDAnalysis/transformations/rotate.py | 14 ++++---- .../transformations/test_rotate.py | 32 +++++++------------ 2 files changed, 18 insertions(+), 28 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 3de862dcf2d..12b1069e442 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -36,7 +36,7 @@ from ..lib.transformations import rotation_matrix -def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=None): +def rotateby(angle, direction, point=None, center_of="geometry", wrap=False, ag=None): ''' 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 @@ -71,7 +71,7 @@ def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=Non 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 + center_of: str, optional used to choose the method of centering on the given atom group. Can be 'geometry' or 'mass' wrap: bool, optional @@ -100,14 +100,14 @@ def rotateby(angle, direction, point=None, center="geometry", wrap=False, ag=Non raise ValueError('{} is not a valid point'.format(point)) elif ag: try: - if center == 'geometry': + if center_of == 'geometry': center_method = partial(ag.center_of_geometry, pbc=pbc_arg) - elif center == 'mass': + elif center_of == 'mass': center_method = partial(ag.center_of_mass, pbc=pbc_arg) else: - raise ValueError('{} is not a valid argument for center'.format(center)) + raise ValueError('{} is not a valid argument for center'.format(center_of)) except AttributeError: - if center == 'mass': + if center_of == 'mass': raise AttributeError('{} is not an AtomGroup object with masses'.format(ag)) else: raise ValueError('{} is not an AtomGroup object'.format(ag)) @@ -120,7 +120,7 @@ def wrapped(ts): else: position = point matrix = rotation_matrix(angle, direction, position) - rotation = matrix[:3, :3] + rotation = matrix[:3, :3].T translation = matrix[:3, 3] ts.positions= np.dot(ts.positions, rotation) ts.positions += translation diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 9422b4787ca..51619281ff0 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -71,9 +71,7 @@ def test_rotateby_custom_position(rotate_universes): pos = [0,0,0] angle = 90 matrix = rotation_matrix(np.deg2rad(angle), vector, pos) - rotation = matrix[:3, :3] - translation = matrix[:3, 3] - ref.positions = np.dot(ref.positions, rotation) + translation + ref_u.atoms.transform(matrix) transformed = rotateby(angle, vector, point=pos)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) @@ -88,11 +86,9 @@ def test_rotateby_atomgroup_cog_nopbc(rotate_universes): vector = [1,0,0] angle = 90 matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) - rotation = matrix[:3, :3] - translation = matrix[:3, 3] - ref.positions = np.dot(ref.positions, rotation) + translation + 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, center_of='geometry')(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) def test_rotateby_atomgroup_com_nopbc(rotate_universes): @@ -107,10 +103,8 @@ def test_rotateby_atomgroup_com_nopbc(rotate_universes): selection = trans_u.residues[0].atoms center_pos = selection.center_of_mass() matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) - rotation = matrix[:3, :3] - translation = matrix[:3, 3] - ref.positions = np.dot(ref.positions, rotation) + translation - transformed = rotateby(angle, vector, ag=selection, center='mass')(trans) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, center_of='mass')(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) def test_rotateby_atomgroup_cog_pbc(rotate_universes): @@ -125,10 +119,8 @@ def test_rotateby_atomgroup_cog_pbc(rotate_universes): selection = trans_u.residues[0].atoms center_pos = selection.center_of_geometry(pbc=True) matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) - rotation = matrix[:3, :3] - translation = matrix[:3, 3] - ref.positions = np.dot(ref.positions, rotation) + translation - transformed = rotateby(angle, vector, ag=selection, center='geometry', wrap=True)(trans) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, center_of='geometry', wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) def test_rotateby_atomgroup_com_pbc(rotate_universes): @@ -143,10 +135,8 @@ def test_rotateby_atomgroup_com_pbc(rotate_universes): selection = trans_u.residues[0].atoms center_pos = selection.center_of_mass(pbc=True) matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) - rotation = matrix[:3, :3] - translation = matrix[:3, 3] - ref.positions = np.dot(ref.positions, rotation) + translation - transformed = rotateby(angle, vector, ag=selection, center='mass', wrap=True)(trans) + ref_u.atoms.transform(matrix) + transformed = rotateby(angle, vector, ag=selection, center_of='mass', wrap=True)(trans) assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) def test_rotateby_bad_ag(rotate_universes): @@ -190,7 +180,7 @@ def test_rotateby_bad_center(rotate_universes): vector = [0, 0, 1] bad_center = " " with pytest.raises(ValueError): - rotateby(angle, vector, ag = ag, center=bad_center)(ts) + rotateby(angle, vector, ag = ag, center_of=bad_center)(ts) def test_rotateby_no_masses(rotate_universes): # this universe as a box size zero @@ -201,7 +191,7 @@ def test_rotateby_no_masses(rotate_universes): vector = [0, 0, 1] bad_center = "mass" with pytest.raises(AttributeError): - rotateby(angle, vector, ag = ag, center=bad_center)(ts) + rotateby(angle, vector, ag = ag, center_of=bad_center)(ts) def test_rotateby_no_args(rotate_universes): # this universe as a box size zero From 6cffae9ab867940d56a480c80e1b8b799decb8c3 Mon Sep 17 00:00:00 2001 From: davidercruz Date: Sun, 22 Jul 2018 01:43:50 +0100 Subject: [PATCH 03/19] rotateby now accepts weights --- package/MDAnalysis/transformations/rotate.py | 34 +++++++++---------- .../transformations/test_rotate.py | 16 ++++----- 2 files changed, 25 insertions(+), 25 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 12b1069e442..53a87f38e50 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -35,8 +35,9 @@ from functools import partial from ..lib.transformations import rotation_matrix +from ..lib.util import get_weights -def rotateby(angle, direction, point=None, center_of="geometry", wrap=False, ag=None): +def rotateby(angle, direction, point=None, weights=None, wrap=False, ag=None): ''' 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 @@ -69,11 +70,14 @@ def rotateby(angle, direction, point=None, center_of="geometry", wrap=False, ag= vector that will define the direction of a custom axis of rotation from the provided point. 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_of: str, optional - used to choose the method of centering on the given atom group. Can be 'geometry' - or 'mass' + use the eighted center of an AtomGroup as the point from where the rotation axis + will be defined + 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 @@ -92,7 +96,6 @@ def rotateby(angle, direction, point=None, center_of="geometry", wrap=False, ag= after rotating the trajectory. ''' - pbc_arg = wrap angle = np.deg2rad(angle) if point: point = np.asarray(point, np.float32) @@ -100,17 +103,14 @@ def rotateby(angle, direction, point=None, center_of="geometry", wrap=False, ag= raise ValueError('{} is not a valid point'.format(point)) elif ag: try: - if center_of == 'geometry': - center_method = partial(ag.center_of_geometry, pbc=pbc_arg) - elif center_of == 'mass': - center_method = partial(ag.center_of_mass, pbc=pbc_arg) - else: - raise ValueError('{} is not a valid argument for center'.format(center_of)) + weights = get_weights(ag.atoms, weights=weights) + except (ValueError, TypeError): + raise ValueError("weights must be {'mass', None} or an iterable of the " + "same size as the atomgroup.") + try: + center_method = partial(ag.atoms.center, weights, pbc=wrap) except AttributeError: - if center_of == '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: raise ValueError('A point or an AtomGroup must be specified') diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 51619281ff0..2ed878a7614 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -88,7 +88,7 @@ def test_rotateby_atomgroup_cog_nopbc(rotate_universes): 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_of='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): @@ -104,7 +104,7 @@ def test_rotateby_atomgroup_com_nopbc(rotate_universes): center_pos = selection.center_of_mass() matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) ref_u.atoms.transform(matrix) - transformed = rotateby(angle, vector, ag=selection, center_of='mass')(trans) + 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): @@ -120,7 +120,7 @@ def test_rotateby_atomgroup_cog_pbc(rotate_universes): center_pos = selection.center_of_geometry(pbc=True) matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) ref_u.atoms.transform(matrix) - transformed = rotateby(angle, vector, ag=selection, center_of='geometry', wrap=True)(trans) + 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): @@ -136,7 +136,7 @@ def test_rotateby_atomgroup_com_pbc(rotate_universes): center_pos = selection.center_of_mass(pbc=True) matrix = rotation_matrix(np.deg2rad(angle), vector, center_pos) ref_u.atoms.transform(matrix) - transformed = rotateby(angle, vector, ag=selection, center_of='mass', wrap=True)(trans) + 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): @@ -147,7 +147,7 @@ def test_rotateby_bad_ag(rotate_universes): angle = 90 vector = [0, 0, 1] bad_ag = 1 - with pytest.raises(ValueError): + with pytest.raises(AttributeError): rotateby(angle, vector, ag = bad_ag)(ts) def test_rotateby_bad_position(rotate_universes): @@ -180,7 +180,7 @@ def test_rotateby_bad_center(rotate_universes): vector = [0, 0, 1] bad_center = " " with pytest.raises(ValueError): - rotateby(angle, vector, ag = ag, center_of=bad_center)(ts) + rotateby(angle, vector, ag = ag, weights=bad_center)(ts) def test_rotateby_no_masses(rotate_universes): # this universe as a box size zero @@ -190,8 +190,8 @@ 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_of=bad_center)(ts) + with pytest.raises(ValueError): + rotateby(angle, vector, ag = ag, weights=bad_center)(ts) def test_rotateby_no_args(rotate_universes): # this universe as a box size zero From 9ba25f1f07e173089b2edc449ef559c634d14f59 Mon Sep 17 00:00:00 2001 From: davidercruz Date: Sun, 22 Jul 2018 04:19:36 +0100 Subject: [PATCH 04/19] fixed typos --- .../transformations/test_rotate.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 2ed878a7614..425f5794384 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -171,16 +171,26 @@ def test_rotateby_bad_pbc(rotate_universes): 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 = " " + bad_weights = " " with pytest.raises(ValueError): - rotateby(angle, vector, ag = ag, weights=bad_center)(ts) + rotateby(angle, vector, ag = ag, weights=bad_weights)(ts) def test_rotateby_no_masses(rotate_universes): # this universe as a box size zero From 9e5db309609cd46642d06f647102393eb0ecaede Mon Sep 17 00:00:00 2001 From: davidercruz Date: Sun, 22 Jul 2018 12:51:31 +0100 Subject: [PATCH 05/19] code cleanup --- package/MDAnalysis/transformations/rotate.py | 14 ++++++++------ .../MDAnalysisTests/transformations/test_rotate.py | 6 +++--- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 53a87f38e50..4bd791281a9 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -103,14 +103,16 @@ def rotateby(angle, direction, point=None, weights=None, wrap=False, ag=None): raise ValueError('{} is not a valid point'.format(point)) elif ag: try: - weights = get_weights(ag.atoms, weights=weights) - except (ValueError, TypeError): - raise ValueError("weights must be {'mass', None} or an iterable of the " - "same size as the atomgroup.") - try: - center_method = partial(ag.atoms.center, weights, pbc=wrap) + 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.") + center_method = partial(atoms.center, weights, pbc=wrap) else: raise ValueError('A point or an AtomGroup must be specified') diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 425f5794384..5205f3a6998 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -147,7 +147,7 @@ def test_rotateby_bad_ag(rotate_universes): angle = 90 vector = [0, 0, 1] bad_ag = 1 - with pytest.raises(AttributeError): + with pytest.raises(ValueError): rotateby(angle, vector, ag = bad_ag)(ts) def test_rotateby_bad_position(rotate_universes): @@ -189,7 +189,7 @@ def test_rotateby_bad_weights(rotate_universes, weights): angle = 90 vector = [0, 0, 1] bad_weights = " " - with pytest.raises(ValueError): + with pytest.raises(TypeError): rotateby(angle, vector, ag = ag, weights=bad_weights)(ts) def test_rotateby_no_masses(rotate_universes): @@ -200,7 +200,7 @@ def test_rotateby_no_masses(rotate_universes): angle = 90 vector = [0, 0, 1] bad_center = "mass" - with pytest.raises(ValueError): + with pytest.raises(TypeError): rotateby(angle, vector, ag = ag, weights=bad_center)(ts) def test_rotateby_no_args(rotate_universes): From f8d49fbfbc5bd4080792169be20d624abff3175f Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Sun, 22 Jul 2018 22:33:17 +0100 Subject: [PATCH 06/19] added checks for direction ; added tests for direction --- package/MDAnalysis/transformations/rotate.py | 24 +++++-- .../transformations/test_rotate.py | 64 +++++++++++++++++-- 2 files changed, 76 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index 4bd791281a9..b16292d7ca0 100644 --- a/package/MDAnalysis/transformations/rotate.py +++ b/package/MDAnalysis/transformations/rotate.py @@ -37,20 +37,26 @@ from ..lib.transformations import rotation_matrix from ..lib.util import get_weights -def rotateby(angle, direction, point=None, weights=None, 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: + 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: @@ -60,7 +66,7 @@ def rotateby(angle, direction, point=None, weights=None, wrap=False, ag=None): 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 ---------- @@ -97,10 +103,18 @@ def rotateby(angle, direction, point=None, weights=None, wrap=False, ag=None): ''' angle = np.deg2rad(angle) + 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: 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: atoms = ag.atoms diff --git a/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 5205f3a6998..823c8be4310 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 @@ -74,7 +75,8 @@ def test_rotateby_custom_position(rotate_universes): ref_u.atoms.transform(matrix) transformed = rotateby(angle, vector, point=pos)(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? @@ -91,6 +93,7 @@ def test_rotateby_atomgroup_cog_nopbc(rotate_universes): 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? @@ -106,6 +109,7 @@ def test_rotateby_atomgroup_com_nopbc(rotate_universes): 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 @@ -123,6 +127,7 @@ def test_rotateby_atomgroup_cog_pbc(rotate_universes): 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? @@ -138,8 +143,19 @@ def test_rotateby_atomgroup_com_pbc(rotate_universes): 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,24 +166,56 @@ 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) @@ -191,6 +239,7 @@ def test_rotateby_bad_weights(rotate_universes, weights): 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 @@ -203,6 +252,7 @@ def test_rotateby_no_masses(rotate_universes): 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 ts = rotate_universes[0].trajectory.ts From e60b268f279b93210f19456b0a5738a1c8aae4e4 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Sat, 28 Jul 2018 17:02:39 +0100 Subject: [PATCH 07/19] Improved docs; test point and vector in 2 shapes --- package/MDAnalysis/transformations/rotate.py | 29 +++++++++-------- .../trajectory_transformations.rst | 1 + .../transformations/rotate.rst | 1 + .../transformations/test_rotate.py | 32 ++++++++++++++++--- 4 files changed, 45 insertions(+), 18 deletions(-) create mode 100644 package/doc/sphinx/source/documentation_pages/transformations/rotate.rst diff --git a/package/MDAnalysis/transformations/rotate.py b/package/MDAnalysis/transformations/rotate.py index b16292d7ca0..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 @@ -41,8 +43,8 @@ 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. + of geometry or the center of mass of a user defined AtomGroup, or an array defining + custom coordinates. Examples -------- @@ -74,10 +76,14 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): 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 the eighted center of an AtomGroup as the point from where the rotation axis - will be defined + 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. @@ -88,13 +94,10 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): 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 ------- @@ -110,7 +113,7 @@ def rotateby(angle, direction, point=None, ag=None, weights=None, wrap=False): direction = direction.reshape(3, ) except ValueError: raise ValueError('{} is not a valid direction'.format(direction)) - if point: + 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)) diff --git a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst index 71599ff44f9..6ed38cdcf62 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -96,3 +96,4 @@ e.g. giving a workflow as a keyword argument when defining the universe: .. toctree:: ./transformations/translate + ./transformations/rotate \ No newline at end of file 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/testsuite/MDAnalysisTests/transformations/test_rotate.py b/testsuite/MDAnalysisTests/transformations/test_rotate.py index 823c8be4310..39384712e51 100644 --- a/testsuite/MDAnalysisTests/transformations/test_rotate.py +++ b/testsuite/MDAnalysisTests/transformations/test_rotate.py @@ -61,19 +61,41 @@ def test_rotation_matrix(): matrix = rotation_matrix(np.deg2rad(angle), vector, pos)[:3, :3] assert_array_almost_equal(matrix, ref_matrix, decimal=6) - -def test_rotateby_custom_position(rotate_universes): +@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 = [0,0,0] + 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=pos)(trans) + transformed = rotateby(angle, vector, point=point)(trans) + assert_array_almost_equal(transformed.positions, ref.positions, decimal=6) + + +@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 + point = [0, 0, 0] + angle = 90 + 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) From 31fa88e41ac01ce495a80d8c11e4d065c2de91ba Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Wed, 1 Aug 2018 01:10:54 +0100 Subject: [PATCH 08/19] added wrap and unwrap --- package/MDAnalysis/transformations/wrap.py | 131 +++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 package/MDAnalysis/transformations/wrap.py diff --git a/package/MDAnalysis/transformations/wrap.py b/package/MDAnalysis/transformations/wrap.py new file mode 100644 index 00000000000..2e99da9cf00 --- /dev/null +++ b/package/MDAnalysis/transformations/wrap.py @@ -0,0 +1,131 @@ +# -*- 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 translation --- :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 this group back into the unit cell. + + +-----------+ +-----------+ + | | | | + | 3 | 6 | 6 3 | + | ! | ! | ! ! | + | 1-2-|-5-8 -> |-5-8 1-2-| + | ! | ! | ! ! | + | 4 | 7 | 7 4 | + | | | | + +-----------+ +-----------+ + + Usage + ----- + + .. 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 + | | | | + +-----------+ +-----------+ + + Parameters + ---------- + atomgroup : AtomGroup + The :class:`MDAnalysis.core.groups.AtomGroup` to work with. + The positions of this are modified in place. + """ + + 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 From 9c7b61e74041d1ddaf16e9983ce46326ca05b1b7 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 3 Aug 2018 18:20:37 +0100 Subject: [PATCH 09/19] added tests ; updated docs --- .../MDAnalysis/transformations/__init__.py | 1 + package/MDAnalysis/transformations/wrap.py | 33 +++-- .../trajectory_transformations.rst | 1 + .../transformations/wrap.rst | 1 + .../transformations/test_wrap.py | 129 ++++++++++++++++++ 5 files changed, 156 insertions(+), 9 deletions(-) create mode 100644 package/doc/sphinx/source/documentation_pages/transformations/wrap.rst create mode 100644 testsuite/MDAnalysisTests/transformations/test_wrap.py 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/wrap.py b/package/MDAnalysis/transformations/wrap.py index 2e99da9cf00..101591899f2 100644 --- a/package/MDAnalysis/transformations/wrap.py +++ b/package/MDAnalysis/transformations/wrap.py @@ -21,8 +21,8 @@ # """\ -Wrap/unwrap translation --- :mod:`MDAnalysis.transformations.wrap` -================================================================== +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 @@ -31,14 +31,14 @@ .. autofunction:: wrap .. autofunction:: unwrap - + """ from ..lib._cutil import make_whole def wrap(ag, compound='atoms'): """ - Shift the contents of this group back into the unit cell. + Shift the contents of a given AtomGroup back into the unit cell. :: +-----------+ +-----------+ | | | | @@ -50,8 +50,8 @@ def wrap(ag, compound='atoms'): | | | | +-----------+ +-----------+ - Usage - ----- + Example + ------- .. code-block:: python @@ -65,7 +65,7 @@ def wrap(ag, compound='atoms'): 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. + The group which will be kept together through the shifting process. Notes ----- @@ -98,7 +98,7 @@ def unwrap(ag): 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. + such as calculating the center of mass of the molecule. :: +-----------+ +-----------+ | | | | @@ -110,11 +110,25 @@ def unwrap(ag): | | | | +-----------+ +-----------+ + 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. + The positions of this are modified in place. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + """ try: @@ -129,3 +143,4 @@ def wrapped(ts): 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..0ef663d8dd6 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -96,3 +96,4 @@ e.g. giving a workflow as a keyword argument when defining the universe: .. toctree:: ./transformations/translate + ./transformations/wrap \ 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/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) From 917b60ddccf4e883547cafceb124e8683dbd106c Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 3 Aug 2018 18:25:33 +0100 Subject: [PATCH 10/19] updated CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 05e9b078d1c..d4effd645f7 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -20,6 +20,7 @@ The rules for this file: Enhancements + * Added wrap/unwrap transformations * Added augment functionality to create relevant images of particles in the vicinity of central cell to handle periodic boundary conditions (PR #1977) From c00a04a454da92acc6bf866a0b357f7987aefd20 Mon Sep 17 00:00:00 2001 From: Richard Gowers Date: Sat, 4 Aug 2018 08:15:22 -0500 Subject: [PATCH 11/19] removed use of pytest.approx (#2019) * removed use of pytest.approx * removed mock decorator, playing badly with pytest fixtures --- testsuite/MDAnalysisTests/topology/test_mmtf.py | 9 +++++---- testsuite/MDAnalysisTests/topology/test_pqr.py | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) 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) From 2c9d6c50dfeb75b07f8aa87e299a717536b51931 Mon Sep 17 00:00:00 2001 From: Jonathan Barnoud Date: Sun, 5 Aug 2018 18:50:09 +0200 Subject: [PATCH 12/19] Use stricter type for pkdtree indices (#2025) Using consistently np.int64 for the indices used in pkdtree avoids type errors when using unique_int_1d. This should reduce the number of failing tests on windows. --- package/MDAnalysis/lib/pkdtree.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From fe792a1e0aceb6ccd27ff49aa9f00bcdf607eeb4 Mon Sep 17 00:00:00 2001 From: Tyler Reddy Date: Sun, 5 Aug 2018 10:23:18 -0700 Subject: [PATCH 13/19] BUG: Force 64-bit storage for xdr file offsets * Fixes Issue #2021 with Windows xdr large file offset handling * on the Windows platform xdr file offsets are now forced to use a 64-bit storage type so that, for example, file seeks beyond 4 GB are correctly reported --- package/MDAnalysis/lib/formats/src/xdrfile.c | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) 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; From 6b74bac909141477a6576b4560755708db50c632 Mon Sep 17 00:00:00 2001 From: ayushsuhane <34154224+ayushsuhane@users.noreply.github.com> Date: Sun, 5 Aug 2018 12:45:11 -0700 Subject: [PATCH 14/19] Using Self Capped distance function in guess_bonds (#2006) * added an option to check for pairs within coordinate array using bruteforce_capped_self, pkdtree_capped_self + tests * Modified guess bonds to use capped function * seperated the self_capped_function, added tests, modified the guess_bonds function * modified CHANGELOG, and minor modifications * modified test for increased coverage, minor modifications to the self_capped_distance --- package/CHANGELOG | 6 + package/MDAnalysis/lib/distances.py | 304 ++++++++++++++++-- package/MDAnalysis/topology/guessers.py | 27 +- .../MDAnalysisTests/lib/test_distances.py | 64 +++- .../MDAnalysisTests/topology/test_guessers.py | 11 + 5 files changed, 358 insertions(+), 54 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 05e9b078d1c..08234a6cb6e 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -20,6 +20,12 @@ The rules for this file: Enhancements + + * 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) 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/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/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))) From f92ec5da70805897bd6f874204e822fe76cb8722 Mon Sep 17 00:00:00 2001 From: ayushsuhane <34154224+ayushsuhane@users.noreply.github.com> Date: Sun, 5 Aug 2018 12:49:29 -0700 Subject: [PATCH 15/19] Around keyword can use KDTree for periodic boundary conditions (#2022) * kdtree now works with pbc, modified _apply_kdtree method in pointselections * CHANGELOG --- package/CHANGELOG | 3 ++- package/MDAnalysis/core/selection.py | 16 +++++++++------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 08234a6cb6e..93524527cd4 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -20,7 +20,8 @@ The rules for this file: Enhancements - + * 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) 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) From bf03cd4571dab9ab16d95348dd6805de9ffd24f3 Mon Sep 17 00:00:00 2001 From: Micaela Matta <31009265+micaela-matta@users.noreply.github.com> Date: Mon, 6 Aug 2018 08:25:20 -0500 Subject: [PATCH 16/19] fixed Issue #2001 - support for Tinker periodic boundary box (#2011) * fixes Issue #2001 - support for Tinker periodic boundary box --- package/.pylintrc | 1 - package/CHANGELOG | 3 ++- package/MDAnalysis/coordinates/TXYZ.py | 21 ++++++++++++---- package/MDAnalysis/topology/TXYZParser.py | 20 ++++++++++++++-- .../MDAnalysisTests/coordinates/test_txyz.py | 24 +++++++++++++++++-- .../data/coordinates/new_hexane.arc | 24 +++++++++++++++++++ testsuite/MDAnalysisTests/datafiles.py | 3 ++- 7 files changed, 84 insertions(+), 12 deletions(-) create mode 100644 testsuite/MDAnalysisTests/data/coordinates/new_hexane.arc 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 93524527cd4..1fb8fee83a1 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -14,7 +14,7 @@ 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 @@ -75,6 +75,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/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/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') From 7c6a5275199c2da4e57f0a35eb9d67380a584ae3 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Wed, 1 Aug 2018 01:10:54 +0100 Subject: [PATCH 17/19] added wrap and unwrap --- package/MDAnalysis/transformations/wrap.py | 131 +++++++++++++++++++++ 1 file changed, 131 insertions(+) create mode 100644 package/MDAnalysis/transformations/wrap.py diff --git a/package/MDAnalysis/transformations/wrap.py b/package/MDAnalysis/transformations/wrap.py new file mode 100644 index 00000000000..2e99da9cf00 --- /dev/null +++ b/package/MDAnalysis/transformations/wrap.py @@ -0,0 +1,131 @@ +# -*- 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 translation --- :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 this group back into the unit cell. + + +-----------+ +-----------+ + | | | | + | 3 | 6 | 6 3 | + | ! | ! | ! ! | + | 1-2-|-5-8 -> |-5-8 1-2-| + | ! | ! | ! ! | + | 4 | 7 | 7 4 | + | | | | + +-----------+ +-----------+ + + Usage + ----- + + .. 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 + | | | | + +-----------+ +-----------+ + + Parameters + ---------- + atomgroup : AtomGroup + The :class:`MDAnalysis.core.groups.AtomGroup` to work with. + The positions of this are modified in place. + """ + + 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 From c9627b432d26217343c376cfde40f24e9d5def97 Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 3 Aug 2018 18:20:37 +0100 Subject: [PATCH 18/19] added tests ; updated docs --- .../MDAnalysis/transformations/__init__.py | 1 + package/MDAnalysis/transformations/wrap.py | 33 +++-- .../trajectory_transformations.rst | 3 +- .../transformations/wrap.rst | 1 + .../transformations/test_wrap.py | 129 ++++++++++++++++++ 5 files changed, 157 insertions(+), 10 deletions(-) create mode 100644 package/doc/sphinx/source/documentation_pages/transformations/wrap.rst create mode 100644 testsuite/MDAnalysisTests/transformations/test_wrap.py 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/wrap.py b/package/MDAnalysis/transformations/wrap.py index 2e99da9cf00..101591899f2 100644 --- a/package/MDAnalysis/transformations/wrap.py +++ b/package/MDAnalysis/transformations/wrap.py @@ -21,8 +21,8 @@ # """\ -Wrap/unwrap translation --- :mod:`MDAnalysis.transformations.wrap` -================================================================== +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 @@ -31,14 +31,14 @@ .. autofunction:: wrap .. autofunction:: unwrap - + """ from ..lib._cutil import make_whole def wrap(ag, compound='atoms'): """ - Shift the contents of this group back into the unit cell. + Shift the contents of a given AtomGroup back into the unit cell. :: +-----------+ +-----------+ | | | | @@ -50,8 +50,8 @@ def wrap(ag, compound='atoms'): | | | | +-----------+ +-----------+ - Usage - ----- + Example + ------- .. code-block:: python @@ -65,7 +65,7 @@ def wrap(ag, compound='atoms'): 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. + The group which will be kept together through the shifting process. Notes ----- @@ -98,7 +98,7 @@ def unwrap(ag): 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. + such as calculating the center of mass of the molecule. :: +-----------+ +-----------+ | | | | @@ -110,11 +110,25 @@ def unwrap(ag): | | | | +-----------+ +-----------+ + 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. + The positions of this are modified in place. + + Returns + ------- + MDAnalysis.coordinates.base.Timestep + """ try: @@ -129,3 +143,4 @@ def wrapped(ts): 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 6ed38cdcf62..ea4539c86ba 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -96,4 +96,5 @@ e.g. giving a workflow as a keyword argument when defining the universe: .. toctree:: ./transformations/translate - ./transformations/rotate \ No newline at end of file + ./transformations/rotate + ./transformations/wrap 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/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) From 5af4274215c109ec4fe8dd5329f8c2363e4925ac Mon Sep 17 00:00:00 2001 From: Davide Cruz Date: Fri, 3 Aug 2018 18:25:33 +0100 Subject: [PATCH 19/19] updated CHANGELOG --- package/CHANGELOG | 1 + 1 file changed, 1 insertion(+) diff --git a/package/CHANGELOG b/package/CHANGELOG index 1fb8fee83a1..c491f708128 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -20,6 +20,7 @@ The rules for this file: 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