From e7f6b6f5185cca1d1eee9b1ea787ea47d3805ff6 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sat, 18 Feb 2023 21:11:36 -0800 Subject: [PATCH 01/54] First attempt at a nojump transformation --- .../MDAnalysis/transformations/__init__.py | 1 + package/MDAnalysis/transformations/nojump.py | 100 ++++++++++++++++++ 2 files changed, 101 insertions(+) create mode 100644 package/MDAnalysis/transformations/nojump.py diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index e9508cddb35..f359363157d 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -129,6 +129,7 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager +from .nojump import NoJump from .fit import fit_translation, fit_rot_trans from .wrap import wrap, unwrap from .boxdimensions import set_dimensions diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py new file mode 100644 index 00000000000..3b2e883ebfa --- /dev/null +++ b/package/MDAnalysis/transformations/nojump.py @@ -0,0 +1,100 @@ +# -*- 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. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + +"""\ +No Jump Trajectory Unwrapping --- :mod:`MDAnalysis.transformations.nojump` +======================================================================================= + +Unwraps the trajectory such that atoms never move more than half a periodic box length. +The algorithm used is based on Kulke & Vermaas 2022, 10.1021/acs.jctc.2c00327 + +.. autoclass:: NoJump + +""" +import numpy as np +import warnings + +from .base import TransformationBase + + +class NoJump(TransformationBase): + """ + Returns transformed coordinates for the given timestep so that an atom + does not move more than half the periodic box size, and will not jump + across the periodic boundary. The algorithm used is based on Kulke & + Vermaas 2022, 10.1021/acs.jctc.2c00327, equation B8 for non-orthogonal + systems. + + Example + ------- + + To calculate diffusion, one of the first steps is to compute a trajectory + where the atoms do not cross the periodic boundary. + + .. code-block:: python + + transformation = NoJump() + u.trajectory.add_transformations(transformation) + for ts in u.trajectory: + print(ts.positions) + + In this case, ``ts.positions`` will return the NoJump unwrapped trajectory. + To reverse the process, wrap the coordinates again. + + Returns + ------- + MDAnalysis.coordinates.timestep.Timestep + + """ + + def __init__(self, + max_threads=None, + #NoJump transforms are inherently unparallelizable, since + #it depends on the previous frame's unwrapped coordinates + parallelizable=False): + super().__init__(max_threads=max_threads, + parallelizable=parallelizable) + self.Lprevinverse = np.eye(3) + self.prevwrap = None + self.prev = None + + def _transform(self, ts): + if (self.prev == None): + self.prev = ts.positions.copy() + self.prevwrap = ts.positions.copy() + self.Lprevinverse = np.linalg.inv(ts.triclinic_dimensions) + return ts + #Remember that @ is a shorthand for matrix multiplication. + #np.matmul(a, b) is equivalent to a @ b + alpha = ts.triclinic_dimensions @ self.Lprevinverse + Linverse = np.linalg.inv(ts.triclinic_dimensions) + current = ts.positions.copy() + ts.positions = alpha @ self.prev + + current - alpha @ self.prevwrap + - ts.triclinic_dimensions @ np.round(Linverse @ current + - self.Lprevinverse @ prevwrap) + self.Lprevinverse = Linverse + self.prevwrap = current + self.prev = ts.positions.copy() + + return ts From 96f43c39fc3b0624e6b1167d299c05d4deac2cdb Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 19 Feb 2023 12:14:02 -0500 Subject: [PATCH 02/54] This works now, matching the output from VMD for orthogonal trajectories --- package/MDAnalysis/transformations/nojump.py | 22 ++++++++++++++------ 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 3b2e883ebfa..ce5adfd8e5d 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -67,7 +67,7 @@ class NoJump(TransformationBase): """ - def __init__(self, + def __init__(self, check_continuity=True, max_threads=None, #NoJump transforms are inherently unparallelizable, since #it depends on the previous frame's unwrapped coordinates @@ -77,24 +77,34 @@ def __init__(self, self.Lprevinverse = np.eye(3) self.prevwrap = None self.prev = None + self.old_frame = 0 + self.check_c = check_continuity def _transform(self, ts): - if (self.prev == None): + if (self.prev is None): self.prev = ts.positions.copy() self.prevwrap = ts.positions.copy() self.Lprevinverse = np.linalg.inv(ts.triclinic_dimensions) + self.old_frame = ts.frame return ts + if self.check_c and np.abs(self.old_frame - ts.frame) != 1: + warnings.warn("NoJump transform is only accurate when positions" + "do not move by more than half a box length." + "Currently jumping between frames with a step of more than 1." + "This might be fine, but depending on the trajectory stride," + "this might be inaccurate.", + Warning) #Remember that @ is a shorthand for matrix multiplication. #np.matmul(a, b) is equivalent to a @ b alpha = ts.triclinic_dimensions @ self.Lprevinverse Linverse = np.linalg.inv(ts.triclinic_dimensions) current = ts.positions.copy() - ts.positions = alpha @ self.prev - + current - alpha @ self.prevwrap - - ts.triclinic_dimensions @ np.round(Linverse @ current - - self.Lprevinverse @ prevwrap) + ts.positions = (alpha @ self.prev.T).T + ts.positions += current - (alpha @ self.prevwrap.T).T + ts.positions -= (ts.triclinic_dimensions @ np.round(Linverse @ current.T - self.Lprevinverse @ self.prevwrap.T)).T self.Lprevinverse = Linverse self.prevwrap = current self.prev = ts.positions.copy() + self.old_frame = ts.frame return ts From 14b724ecb34394c45897093b907adf9f040b547f Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 19 Feb 2023 09:50:26 -0800 Subject: [PATCH 03/54] Add self to AUTHORS, add to CHANGELOG, make style happier. --- package/AUTHORS | 1 + package/CHANGELOG | 2 + package/MDAnalysis/transformations/nojump.py | 46 ++++++++++---------- 3 files changed, 27 insertions(+), 22 deletions(-) diff --git a/package/AUTHORS b/package/AUTHORS index 5b347d08539..f713cb32a07 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -206,6 +206,7 @@ Chronological list of authors - Christian Pfaendner - Pratham Chauhan - Meet Brijwani + - Josh Vermaas External code diff --git a/package/CHANGELOG b/package/CHANGELOG index 0a34834cf4e..94f7a6bbd90 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -26,6 +26,8 @@ Fixes (Issue #3336) Enhancements + * Add a nojump transformation, which unwraps trajectories so that particle + paths are continuous. (Issue #3703, PR #4031) * Add distopia distance calculation library bindings as a selectable backend for `calc_bonds` in `MDA.lib.distances`. (Issue #3783, PR #3914) * AuxReaders are now pickle-able and copy-able (Issue #1785, PR #3887) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index ce5adfd8e5d..a80f8ef66ea 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -41,7 +41,7 @@ class NoJump(TransformationBase): """ Returns transformed coordinates for the given timestep so that an atom does not move more than half the periodic box size, and will not jump - across the periodic boundary. The algorithm used is based on Kulke & + across the periodic boundary. The algorithm used is based on Kulke & Vermaas 2022, 10.1021/acs.jctc.2c00327, equation B8 for non-orthogonal systems. @@ -54,12 +54,12 @@ class NoJump(TransformationBase): .. code-block:: python transformation = NoJump() - u.trajectory.add_transformations(transformation) + u.trajectory.add_transformations(transformation) for ts in u.trajectory: print(ts.positions) In this case, ``ts.positions`` will return the NoJump unwrapped trajectory. - To reverse the process, wrap the coordinates again. + To reverse the process, wrap the coordinates again. Returns ------- @@ -67,13 +67,15 @@ class NoJump(TransformationBase): """ - def __init__(self, check_continuity=True, - max_threads=None, - #NoJump transforms are inherently unparallelizable, since - #it depends on the previous frame's unwrapped coordinates - parallelizable=False): - super().__init__(max_threads=max_threads, - parallelizable=parallelizable) + def __init__( + self, + check_continuity=True, + max_threads=None, + # NoJump transforms are inherently unparallelizable, since + # it depends on the previous frame's unwrapped coordinates + parallelizable=False, + ): + super().__init__(max_threads=max_threads, parallelizable=parallelizable) self.Lprevinverse = np.eye(3) self.prevwrap = None self.prev = None @@ -81,27 +83,27 @@ def __init__(self, check_continuity=True, self.check_c = check_continuity def _transform(self, ts): - if (self.prev is None): + if self.prev is None: self.prev = ts.positions.copy() self.prevwrap = ts.positions.copy() self.Lprevinverse = np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame return ts if self.check_c and np.abs(self.old_frame - ts.frame) != 1: - warnings.warn("NoJump transform is only accurate when positions" - "do not move by more than half a box length." - "Currently jumping between frames with a step of more than 1." - "This might be fine, but depending on the trajectory stride," - "this might be inaccurate.", - Warning) - #Remember that @ is a shorthand for matrix multiplication. - #np.matmul(a, b) is equivalent to a @ b + warnings.warn( + "NoJump transform is only accurate when positions" + "do not move by more than half a box length." + "Currently jumping between frames with a step of more than 1." + "This might be fine, but depending on the trajectory stride," + "this might be inaccurate.", + Warning, + ) + # Remember that @ is a shorthand for matrix multiplication. + # np.matmul(a, b) is equivalent to a @ b alpha = ts.triclinic_dimensions @ self.Lprevinverse Linverse = np.linalg.inv(ts.triclinic_dimensions) current = ts.positions.copy() - ts.positions = (alpha @ self.prev.T).T - ts.positions += current - (alpha @ self.prevwrap.T).T - ts.positions -= (ts.triclinic_dimensions @ np.round(Linverse @ current.T - self.Lprevinverse @ self.prevwrap.T)).T + ts.positions = (alpha @ self.prev.T).T + current - (alpha @ self.prevwrap.T).T - (ts.triclinic_dimensions @ np.round(Linverse @ current.T - self.Lprevinverse @ self.prevwrap.T)).T self.Lprevinverse = Linverse self.prevwrap = current self.prev = ts.positions.copy() From 8a0c60d3d3c99d9e06fa5dbacc9d2e50ff08eb81 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 02:07:11 -0500 Subject: [PATCH 04/54] This implementation works for non-orthogonal boxes, and is simpler than I thought --- package/MDAnalysis/transformations/nojump.py | 29 ++++++++++---------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index ce5adfd8e5d..6ee573a1170 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -74,17 +74,13 @@ def __init__(self, check_continuity=True, parallelizable=False): super().__init__(max_threads=max_threads, parallelizable=parallelizable) - self.Lprevinverse = np.eye(3) - self.prevwrap = None self.prev = None self.old_frame = 0 self.check_c = check_continuity def _transform(self, ts): if (self.prev is None): - self.prev = ts.positions.copy() - self.prevwrap = ts.positions.copy() - self.Lprevinverse = np.linalg.inv(ts.triclinic_dimensions) + self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame return ts if self.check_c and np.abs(self.old_frame - ts.frame) != 1: @@ -96,15 +92,20 @@ def _transform(self, ts): Warning) #Remember that @ is a shorthand for matrix multiplication. #np.matmul(a, b) is equivalent to a @ b - alpha = ts.triclinic_dimensions @ self.Lprevinverse - Linverse = np.linalg.inv(ts.triclinic_dimensions) - current = ts.positions.copy() - ts.positions = (alpha @ self.prev.T).T - ts.positions += current - (alpha @ self.prevwrap.T).T - ts.positions -= (ts.triclinic_dimensions @ np.round(Linverse @ current.T - self.Lprevinverse @ self.prevwrap.T)).T - self.Lprevinverse = Linverse - self.prevwrap = current - self.prev = ts.positions.copy() + L = ts.triclinic_dimensions + Linverse = np.linalg.inv(L) + #Convert into reduced coordinate space + fcurrent = ts.positions @ Linverse + fprev = self.prev #Previous unwrapped coordinates in reduced box coordinates. + #Calculate the new positions in reduced coordinate space (Equation B6 from + #10.1021/acs.jctc.2c00327). As it turns out, the displacement term can + #be moved inside the round function in this coordinate space, as the + #difference between wrapped and unwrapped coordinates is an integer. + newpositions = fcurrent - np.round(fcurrent - fprev) + #Convert back into real space + ts.positions = newpositions @ L + #Set things we need to save. + self.prev = newpositions #Note that this is in reduced coordinate space. self.old_frame = ts.frame return ts From d16bbda43c942a4929f45f097bc5dc023eea151a Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 02:32:41 -0500 Subject: [PATCH 05/54] Added a test --- .../transformations/test_nojump.py | 47 +++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 testsuite/MDAnalysisTests/transformations/test_nojump.py diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py new file mode 100644 index 00000000000..65af88931cd --- /dev/null +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -0,0 +1,47 @@ +### + +import numpy as np +import pytest +from numpy.testing import assert_array_almost_equal + +import MDAnalysis as md +from MDAnalysis.transformations import NoJump +from MDAnalysisTests import datafiles + +@pytest.fixture() +def nojump_universes(): + ''' + Create the universe objects for the tests. + ''' + u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + transformation = NoJump() + u.trajectory.add_transformations(transformation) + return u + +@pytest.fixture() +def nojump_universes_nocheck(): + ''' + Create the universe objects for the tests. + Do not check for continunity + ''' + u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + transformation = NoJump(check_continuity=False) + u.trajectory.add_transformations(transformation) + return u + +def test_nojump_fwd(nojump_universes): + ''' + Test if the nojump transform is returning the correct + values when iterating forwards over the sample trajectory. + ''' + ref_matrix_fwd1 = np.asarray([9.159101 -6.2827725 5.156695]) + ref_matrix_fwd2 = np.asarray([5.455143 -9.155944 6.155884]) + size = (nojump_universes.trajectory.ts.positions.shape[0], + nojump_universes.trajectory.ts.positions.shape[1], + len(nojump_universes.trajectory)) + parr = np.empty(size) + for ts in nojump_universes.trajectory: + parr[...,ts.frame] = ts.positions.copy() + #Atoms 166 and 362 happen to move alot. + assert_array_almost_equal(ref_matrix_fwd1, parr[166,:,-1], decimal=5) + assert_array_almost_equal(ref_matrix_fwd2, parr[362,:,-1], decimal=5) From e05ecfd771a44d5cacc4f89627f28ec5caa45ba9 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 02:36:48 -0500 Subject: [PATCH 06/54] Updates for Linter --- package/MDAnalysis/transformations/nojump.py | 22 ++++++++++--------- .../transformations/test_nojump.py | 4 ++-- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 66aac443c9d..b1db4c67bda 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -59,7 +59,7 @@ class NoJump(TransformationBase): print(ts.positions) In this case, ``ts.positions`` will return the NoJump unwrapped trajectory. - To reverse the process, wrap the coordinates again. + To reverse the process, wrap the coordinates again. Returns ------- @@ -81,19 +81,21 @@ def __init__( self.check_c = check_continuity def _transform(self, ts): - if (self.prev is None): + if self.prev is None: self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame return ts if self.check_c and np.abs(self.old_frame - ts.frame) != 1: - warnings.warn("NoJump transform is only accurate when positions" - "do not move by more than half a box length." - "Currently jumping between frames with a step of more than 1." - "This might be fine, but depending on the trajectory stride," - "this might be inaccurate.", - Warning) - #Remember that @ is a shorthand for matrix multiplication. - #np.matmul(a, b) is equivalent to a @ b + warnings.warn( + "NoJump transform is only accurate when positions" + "do not move by more than half a box length." + "Currently jumping between frames with a step of more than 1." + "This might be fine, but depending on the trajectory stride," + "this might be inaccurate.", + Warning, + ) + # Remember that @ is a shorthand for matrix multiplication. + # np.matmul(a, b) is equivalent to a @ b L = ts.triclinic_dimensions Linverse = np.linalg.inv(L) #Convert into reduced coordinate space diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 65af88931cd..20cfbab52f5 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -34,8 +34,8 @@ def test_nojump_fwd(nojump_universes): Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. ''' - ref_matrix_fwd1 = np.asarray([9.159101 -6.2827725 5.156695]) - ref_matrix_fwd2 = np.asarray([5.455143 -9.155944 6.155884]) + ref_matrix_fwd1 = np.asarray([9.159101, -6.2827725, 5.156695]) + ref_matrix_fwd2 = np.asarray([5.455143 , -9.155944 , 6.155884]) size = (nojump_universes.trajectory.ts.positions.shape[0], nojump_universes.trajectory.ts.positions.shape[1], len(nojump_universes.trajectory)) From 701cb10f0e5c3c1808118263b7d7d78aae15a2d2 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 02:46:42 -0500 Subject: [PATCH 07/54] How about now linter? --- package/MDAnalysis/transformations/nojump.py | 18 +++++----- .../transformations/test_nojump.py | 34 ++++++++++--------- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index b1db4c67bda..a710aa57dc5 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -98,18 +98,18 @@ def _transform(self, ts): # np.matmul(a, b) is equivalent to a @ b L = ts.triclinic_dimensions Linverse = np.linalg.inv(L) - #Convert into reduced coordinate space + # Convert into reduced coordinate space fcurrent = ts.positions @ Linverse - fprev = self.prev #Previous unwrapped coordinates in reduced box coordinates. - #Calculate the new positions in reduced coordinate space (Equation B6 from - #10.1021/acs.jctc.2c00327). As it turns out, the displacement term can - #be moved inside the round function in this coordinate space, as the - #difference between wrapped and unwrapped coordinates is an integer. + fprev = self.prev # Previous unwrapped coordinates in reduced box coordinates. + # Calculate the new positions in reduced coordinate space (Equation B6 from + # 10.1021/acs.jctc.2c00327). As it turns out, the displacement term can + # be moved inside the round function in this coordinate space, as the + # difference between wrapped and unwrapped coordinates is an integer. newpositions = fcurrent - np.round(fcurrent - fprev) - #Convert back into real space + # Convert back into real space ts.positions = newpositions @ L - #Set things we need to save. - self.prev = newpositions #Note that this is in reduced coordinate space. + # Set things we need to save for the next frame. + self.prev = newpositions # Note that this is in reduced coordinate space. self.old_frame = ts.frame return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 20cfbab52f5..49357b40530 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -6,42 +6,44 @@ import MDAnalysis as md from MDAnalysis.transformations import NoJump -from MDAnalysisTests import datafiles +from MDAnalysisTests import datafiles as data @pytest.fixture() def nojump_universes(): - ''' + """ Create the universe objects for the tests. - ''' - u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + """ + u = md.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) transformation = NoJump() u.trajectory.add_transformations(transformation) return u @pytest.fixture() def nojump_universes_nocheck(): - ''' + """ Create the universe objects for the tests. Do not check for continunity - ''' - u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + """ + u = md.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) transformation = NoJump(check_continuity=False) u.trajectory.add_transformations(transformation) return u def test_nojump_fwd(nojump_universes): - ''' + """ Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. - ''' - ref_matrix_fwd1 = np.asarray([9.159101, -6.2827725, 5.156695]) - ref_matrix_fwd2 = np.asarray([5.455143 , -9.155944 , 6.155884]) - size = (nojump_universes.trajectory.ts.positions.shape[0], - nojump_universes.trajectory.ts.positions.shape[1], - len(nojump_universes.trajectory)) - parr = np.empty(size) + """ + ref_matrix_fwd1 = np.asarray([9.159101, -6.2827725, 5.156695]) + ref_matrix_fwd2 = np.asarray([5.455143, -9.155944, 6.155884]) + size = ( + nojump_universes.trajectory.ts.positions.shape[0], + nojump_universes.trajectory.ts.positions.shape[1], + len(nojump_universes.trajectory), + ) + parr = np.empty(size) for ts in nojump_universes.trajectory: parr[...,ts.frame] = ts.positions.copy() - #Atoms 166 and 362 happen to move alot. + # Atoms 166 and 362 happen to move alot. assert_array_almost_equal(ref_matrix_fwd1, parr[166,:,-1], decimal=5) assert_array_almost_equal(ref_matrix_fwd2, parr[362,:,-1], decimal=5) From 5a4bd77c8daf47155d84a91d153942f1da0e3cca Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 02:50:20 -0500 Subject: [PATCH 08/54] Whitespace fixes --- package/MDAnalysis/transformations/nojump.py | 2 +- testsuite/MDAnalysisTests/transformations/test_nojump.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index a710aa57dc5..e247f6f7487 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -109,7 +109,7 @@ def _transform(self, ts): # Convert back into real space ts.positions = newpositions @ L # Set things we need to save for the next frame. - self.prev = newpositions # Note that this is in reduced coordinate space. + self.prev = newpositions # Note that this is in reduced coordinate space. self.old_frame = ts.frame return ts diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 49357b40530..9b0e17721d8 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -41,9 +41,9 @@ def test_nojump_fwd(nojump_universes): nojump_universes.trajectory.ts.positions.shape[1], len(nojump_universes.trajectory), ) - parr = np.empty(size) + parr = np.empty(size) for ts in nojump_universes.trajectory: - parr[...,ts.frame] = ts.positions.copy() + parr[..., ts.frame] = ts.positions.copy() # Atoms 166 and 362 happen to move alot. assert_array_almost_equal(ref_matrix_fwd1, parr[166,:,-1], decimal=5) - assert_array_almost_equal(ref_matrix_fwd2, parr[362,:,-1], decimal=5) + assert_array_almost_equal(ref_matrix_fwd2, parr[362,:,-1], decimal=5) From 88f43952a0bc69e375a88730d826c8fb08e1652d Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 03:03:00 -0500 Subject: [PATCH 09/54] Copy the right value --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 9b0e17721d8..e3fa160b7fb 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -34,7 +34,7 @@ def test_nojump_fwd(nojump_universes): Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. """ - ref_matrix_fwd1 = np.asarray([9.159101, -6.2827725, 5.156695]) + ref_matrix_fwd1 = np.asarray([-3.42261, 11.28495, -1.37211]) ref_matrix_fwd2 = np.asarray([5.455143, -9.155944, 6.155884]) size = ( nojump_universes.trajectory.ts.positions.shape[0], From 081c1273faf8db65b708e7e9dc65df13f295c998 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 03:06:01 -0500 Subject: [PATCH 10/54] Nope, made a mistake with the other one too --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index e3fa160b7fb..f369335c433 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -35,7 +35,7 @@ def test_nojump_fwd(nojump_universes): values when iterating forwards over the sample trajectory. """ ref_matrix_fwd1 = np.asarray([-3.42261, 11.28495, -1.37211]) - ref_matrix_fwd2 = np.asarray([5.455143, -9.155944, 6.155884]) + ref_matrix_fwd2 = np.asarray([3.674243, -8.725193, -0.07884017]) size = ( nojump_universes.trajectory.ts.positions.shape[0], nojump_universes.trajectory.ts.positions.shape[1], From dc060527cfd764afa5a4efcc6cdbe011f5a490e0 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 10:24:48 -0800 Subject: [PATCH 11/54] Add references, make darker happy. --- package/MDAnalysis/transformations/nojump.py | 20 ++++++++++++++----- package/doc/sphinx/source/references.bib | 11 ++++++++++ .../transformations/test_nojump.py | 18 ++++++++++++----- 3 files changed, 39 insertions(+), 10 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index e247f6f7487..5fe9b395b6f 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -22,11 +22,11 @@ # """\ -No Jump Trajectory Unwrapping --- :mod:`MDAnalysis.transformations.nojump` +No Jump Trajectory Unwrapping --- :mod:`MDAnalysis.transformations.nojump` ======================================================================================= Unwraps the trajectory such that atoms never move more than half a periodic box length. -The algorithm used is based on Kulke & Vermaas 2022, 10.1021/acs.jctc.2c00327 +The algorithm used is based on :cite:p:`Kulke2022`. .. autoclass:: NoJump @@ -36,14 +36,16 @@ from .base import TransformationBase +@due.dcite(Doi("10.1021/acs.jctc.2c00327"), + description="Works through the orthogonal case, " + "and proposes the non-orthogonal appraoch.", path=__name__) class NoJump(TransformationBase): """ Returns transformed coordinates for the given timestep so that an atom does not move more than half the periodic box size, and will not jump - across the periodic boundary. The algorithm used is based on Kulke & - Vermaas 2022, 10.1021/acs.jctc.2c00327, equation B8 for non-orthogonal - systems. + across the periodic boundary. The algorithm used is based on :cite:p:`Kulke2022`, + equation B6 for non-orthogonal systems. Example ------- @@ -65,6 +67,14 @@ class NoJump(TransformationBase): ------- MDAnalysis.coordinates.timestep.Timestep + References + ---------- + .. bibliography:: + :filter: False + :style: MDA + + Kulke2022 + """ def __init__( diff --git a/package/doc/sphinx/source/references.bib b/package/doc/sphinx/source/references.bib index 49ce246bdc8..68c95e9133b 100644 --- a/package/doc/sphinx/source/references.bib +++ b/package/doc/sphinx/source/references.bib @@ -759,3 +759,14 @@ @incollection{Waveren2005 booktitle = {}, id = {transformations} } + +@article{Kulke2022, + title = {Reversible Unwrapping Algorithm for Constant-Pressure Molecular Dynamics Simulations}, + author = {Kulke, Martin and Vermaas, Josh V.}, + year = {2022}, + journal = {Journal of Chemical Theory and Computation}, + volume = {18}, + number = {10}, + pages = {6161--6171}, + doi = {10.1021/acs.jctc.2c00327} +} diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index f369335c433..30746c77015 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -4,36 +4,44 @@ import pytest from numpy.testing import assert_array_almost_equal -import MDAnalysis as md +import MDAnalysis as mda from MDAnalysis.transformations import NoJump from MDAnalysisTests import datafiles as data + @pytest.fixture() def nojump_universes(): """ Create the universe objects for the tests. """ - u = md.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) transformation = NoJump() u.trajectory.add_transformations(transformation) return u + @pytest.fixture() def nojump_universes_nocheck(): """ Create the universe objects for the tests. Do not check for continunity """ - u = md.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) transformation = NoJump(check_continuity=False) u.trajectory.add_transformations(transformation) return u + def test_nojump_fwd(nojump_universes): """ Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. """ + # These were determined based on determining what the TRICLINIC system + # would return on a validated version that worked on triclinic and + # orthogonal systems for a small water trajectory. + # These are specific to atoms 166 and 362, which moved the most in the + # trajectory. ref_matrix_fwd1 = np.asarray([-3.42261, 11.28495, -1.37211]) ref_matrix_fwd2 = np.asarray([3.674243, -8.725193, -0.07884017]) size = ( @@ -45,5 +53,5 @@ def test_nojump_fwd(nojump_universes): for ts in nojump_universes.trajectory: parr[..., ts.frame] = ts.positions.copy() # Atoms 166 and 362 happen to move alot. - assert_array_almost_equal(ref_matrix_fwd1, parr[166,:,-1], decimal=5) - assert_array_almost_equal(ref_matrix_fwd2, parr[362,:,-1], decimal=5) + assert_array_almost_equal(ref_matrix_fwd1, parr[166, :, -1], decimal=5) + assert_array_almost_equal(ref_matrix_fwd2, parr[362, :, -1], decimal=5) From 65c95c0b725ee34aff18fedbb8e04f4a952990f7 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 10:31:27 -0800 Subject: [PATCH 12/54] Spellcheck, actually define due and Doi --- package/MDAnalysis/transformations/nojump.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 5fe9b395b6f..95b8183e815 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -35,10 +35,8 @@ import warnings from .base import TransformationBase +from ..due import due, Doi -@due.dcite(Doi("10.1021/acs.jctc.2c00327"), - description="Works through the orthogonal case, " - "and proposes the non-orthogonal appraoch.", path=__name__) class NoJump(TransformationBase): """ @@ -77,6 +75,10 @@ class NoJump(TransformationBase): """ + @due.dcite(Doi("10.1021/acs.jctc.2c00327"), + description="Works through the orthogonal case, " + "and proposes the non-orthogonal approach.", path=__name__) + def __init__( self, check_continuity=True, From 5acec8a3aac0982c08ec29d1a5fd23dba9afa295 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 21:16:11 -0800 Subject: [PATCH 13/54] Update test_nojump.py with a much nicer synthetic example --- .../transformations/test_nojump.py | 65 +++++++++---------- 1 file changed, 31 insertions(+), 34 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 30746c77015..7a666067aa1 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -1,8 +1,6 @@ -### - import numpy as np import pytest -from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_allclose import MDAnalysis as mda from MDAnalysis.transformations import NoJump @@ -10,48 +8,47 @@ @pytest.fixture() -def nojump_universes(): +def nojump_universe(): """ Create the universe objects for the tests. """ - u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) - transformation = NoJump() - u.trajectory.add_transformations(transformation) + u = mda.Universe.empty(N, trajectory=True) + coordinates = np.empty((100, # number of frames + u.atoms.n_atoms, + 3)) + coordinates[::3,0] = 0 * np.ones(3) + coordinates[1::3,0] = 0.3 * np.ones(3) + coordinates[2::3,0] = 0.6 * np.ones(3) + u.load_new(coordinates, order='fac') return u -@pytest.fixture() -def nojump_universes_nocheck(): +def test_nojump_orthogonal_fwd(nojump_universe): """ - Create the universe objects for the tests. - Do not check for continunity + Test if the nojump transform is returning the correct + values when iterating forwards over the sample trajectory. """ - u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) - transformation = NoJump(check_continuity=False) - u.trajectory.add_transformations(transformation) - return u + u = nojump_universe + dim = np.asarray([1, 1, 1, 90, 90, 90], np.float32) + workflow = [mda.transformations.boxdimensions.set_dimensions(dim), NoJump()] + u.add_transformations(*workflow) + transformed_coordinates = u.trajectory.timeseries()[0] + # Step is 1 unit every 3 steps. After 99 steps from the origin, + # we'll end up at 33. + assert_allclose(transformed_coordinates[99], 33 * np.ones(3)) -def test_nojump_fwd(nojump_universes): +def test_nojump_nonorthogonal_fwd(nojump_universes_nonorthogonal): """ Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. """ - # These were determined based on determining what the TRICLINIC system - # would return on a validated version that worked on triclinic and - # orthogonal systems for a small water trajectory. - # These are specific to atoms 166 and 362, which moved the most in the - # trajectory. - ref_matrix_fwd1 = np.asarray([-3.42261, 11.28495, -1.37211]) - ref_matrix_fwd2 = np.asarray([3.674243, -8.725193, -0.07884017]) - size = ( - nojump_universes.trajectory.ts.positions.shape[0], - nojump_universes.trajectory.ts.positions.shape[1], - len(nojump_universes.trajectory), - ) - parr = np.empty(size) - for ts in nojump_universes.trajectory: - parr[..., ts.frame] = ts.positions.copy() - # Atoms 166 and 362 happen to move alot. - assert_array_almost_equal(ref_matrix_fwd1, parr[166, :, -1], decimal=5) - assert_array_almost_equal(ref_matrix_fwd2, parr[362, :, -1], decimal=5) + u = nojump_universe + dim = np.asarray([1, 1, 1, 90, 60, 90], np.float32) + workflow = [mda.transformations.boxdimensions.set_dimensions(dim), NoJump()] + u.add_transformations(*workflow) + transformed_coordinates = u.trajectory.timeseries()[0] + # Step is 1 unit every 3 steps. After 99 steps from the origin, + # we'll end up at 33. However, since the unit cell is non-orthogonal, + # we'll end up at a distorted place. + assert_allclose(transformed_coordinates[99], 33 * np.array([0.5, 1, np.sqrt(3)/2])) From 5b2f067b0e924653f962410cb8a94709381853cd Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 21:50:56 -0800 Subject: [PATCH 14/54] Add to the docs --- package/MDAnalysis/transformations/nojump.py | 4 ++-- .../source/documentation_pages/trajectory_transformations.rst | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 95b8183e815..3f9c84e39f9 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -76,7 +76,7 @@ class NoJump(TransformationBase): """ @due.dcite(Doi("10.1021/acs.jctc.2c00327"), - description="Works through the orthogonal case, " + description="Works through the orthogonal case for unwrapping, " "and proposes the non-orthogonal approach.", path=__name__) def __init__( @@ -104,7 +104,7 @@ def _transform(self, ts): "Currently jumping between frames with a step of more than 1." "This might be fine, but depending on the trajectory stride," "this might be inaccurate.", - Warning, + UserWarning, ) # Remember that @ is a shorthand for matrix multiplication. # np.matmul(a, b) is equivalent to a @ b diff --git a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst index aafc31976f8..425762bb594 100644 --- a/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst +++ b/package/doc/sphinx/source/documentation_pages/trajectory_transformations.rst @@ -270,5 +270,6 @@ Currently implemented transformations ./transformations/positionaveraging ./transformations/fit ./transformations/wrap + ./transformations/nojump ./transformations/boxdimensions From f0d57ce17c41cc51af09c93d3c71786a0d222a91 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 22:04:45 -0800 Subject: [PATCH 15/54] Add a check that we don't go backwards --- package/MDAnalysis/transformations/nojump.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 3f9c84e39f9..1f264f6bee5 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -90,6 +90,7 @@ def __init__( super().__init__(max_threads=max_threads, parallelizable=parallelizable) self.prev = None self.old_frame = 0 + self.older_frame = 'A' self.check_c = check_continuity def _transform(self, ts): @@ -97,6 +98,16 @@ def _transform(self, ts): self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame return ts + if self.check_c and self.older_frame != 'A' and (self.old_frame - self.older_frame) != (ts.frame - self.old_frame): + warnings.warn( + "NoJump detected that the interval between frames is unequal." + "We are resetting the reference frame to the current timestep.", + UserWarning, + ) + self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) + self.old_frame = ts.frame + self.older_frame = 'A' + return ts if self.check_c and np.abs(self.old_frame - ts.frame) != 1: warnings.warn( "NoJump transform is only accurate when positions" @@ -122,6 +133,7 @@ def _transform(self, ts): ts.positions = newpositions @ L # Set things we need to save for the next frame. self.prev = newpositions # Note that this is in reduced coordinate space. + self.older_frame = self.old_frame self.old_frame = ts.frame return ts From 146cc07725a491ee3ab058c274e013b6b6739a91 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 22:12:48 -0800 Subject: [PATCH 16/54] Automake the nojump documentation --- .../sphinx/source/documentation_pages/transformations/nojump.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 package/doc/sphinx/source/documentation_pages/transformations/nojump.rst diff --git a/package/doc/sphinx/source/documentation_pages/transformations/nojump.rst b/package/doc/sphinx/source/documentation_pages/transformations/nojump.rst new file mode 100644 index 00000000000..d1543ca4cb7 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/transformations/nojump.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.transformations.nojump From c8263f1ebc0a8ec9bfe5d0dc708080aec05c4a62 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 22:37:47 -0800 Subject: [PATCH 17/54] Test every third value --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 7a666067aa1..2cf6d77559d 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -35,10 +35,10 @@ def test_nojump_orthogonal_fwd(nojump_universe): transformed_coordinates = u.trajectory.timeseries()[0] # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. - assert_allclose(transformed_coordinates[99], 33 * np.ones(3)) + assert_allclose(transformed_coordinates[::3], np.outer(np.arange(33.5),np.ones(3))) -def test_nojump_nonorthogonal_fwd(nojump_universes_nonorthogonal): +def test_nojump_nonorthogonal_fwd(nojump_universe): """ Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. @@ -51,4 +51,4 @@ def test_nojump_nonorthogonal_fwd(nojump_universes_nonorthogonal): # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. However, since the unit cell is non-orthogonal, # we'll end up at a distorted place. - assert_allclose(transformed_coordinates[99], 33 * np.array([0.5, 1, np.sqrt(3)/2])) + assert_allclose(transformed_coordinates[::3], np.outer(np.arange(33.5),np.array([0.5, 1, np.sqrt(3)/2]))) From 84527c8b8b0228a7247d79f72d30ce3a9052c519 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 22:43:13 -0800 Subject: [PATCH 18/54] Nojump --- package/MDAnalysis/transformations/nojump.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 1f264f6bee5..4d73272f263 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -26,6 +26,9 @@ ======================================================================================= Unwraps the trajectory such that atoms never move more than half a periodic box length. +The consequence of this is that particles will diffuse across periodic boundaries when +needed. This unwrapping method is suitable as a preprocessing step to calculate +molecular diffusion. The algorithm used is based on :cite:p:`Kulke2022`. .. autoclass:: NoJump @@ -41,8 +44,8 @@ class NoJump(TransformationBase): """ Returns transformed coordinates for the given timestep so that an atom - does not move more than half the periodic box size, and will not jump - across the periodic boundary. The algorithm used is based on :cite:p:`Kulke2022`, + does not move more than half the periodic box size, and will move + across periodic boundary edges. The algorithm used is based on :cite:p:`Kulke2022`, equation B6 for non-orthogonal systems. Example From 635d15a9d2b31efcb83cb8668ce5390ac44f39fd Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 22:50:10 -0800 Subject: [PATCH 19/54] Test values in more places --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 2cf6d77559d..8ad94639547 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -16,9 +16,9 @@ def nojump_universe(): coordinates = np.empty((100, # number of frames u.atoms.n_atoms, 3)) - coordinates[::3,0] = 0 * np.ones(3) - coordinates[1::3,0] = 0.3 * np.ones(3) - coordinates[2::3,0] = 0.6 * np.ones(3) + coordinates[::3,0] = 0 * np.ones(3) / 3 + coordinates[1::3,0] = 1 * np.ones(3) / 3 + coordinates[2::3,0] = 2 * np.ones(3) / 3 u.load_new(coordinates, order='fac') return u @@ -35,7 +35,7 @@ def test_nojump_orthogonal_fwd(nojump_universe): transformed_coordinates = u.trajectory.timeseries()[0] # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. - assert_allclose(transformed_coordinates[::3], np.outer(np.arange(33.5),np.ones(3))) + assert_allclose(transformed_coordinates, np.outer(np.linspace(0,33,100),np.ones(3))) def test_nojump_nonorthogonal_fwd(nojump_universe): From 662a31609fde5260631451d7b913f21e402c0221 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 23:44:00 -0800 Subject: [PATCH 20/54] Actually add constant velocity test --- .../transformations/test_nojump.py | 32 ++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 8ad94639547..a3e269a8dfc 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -12,7 +12,7 @@ def nojump_universe(): """ Create the universe objects for the tests. """ - u = mda.Universe.empty(N, trajectory=True) + u = mda.Universe.empty(1, trajectory=True) coordinates = np.empty((100, # number of frames u.atoms.n_atoms, 3)) @@ -23,6 +23,26 @@ def nojump_universe(): return u +@pytest.fixture() +def nojump_constantvel_universe(): + """ + Create the universe objects for the tests. + """ + Natom = 1 + Nframe = 100 + coordinates = np.empty((Nframe, # number of frames + u.atoms.n_atoms, + 3)) + coordinates[:,0,0] = np.linspace(0,45,Nframe) + coordinates[:,0,1] = np.linspace(-19,15,Nframe)[::-1] + coordinates[:,0,2] = np.linspace(-10,10,Nframe) + reference = mda.Universe.empty(Natom, trajectory=True) + reference.load_new(coordinates, order='fac') + towrap = mda.Universe.empty(Natom, trajectory=True) + towrap.load_new(coordinates, order='fac') + return reference, towrap + + def test_nojump_orthogonal_fwd(nojump_universe): """ Test if the nojump transform is returning the correct @@ -52,3 +72,13 @@ def test_nojump_nonorthogonal_fwd(nojump_universe): # we'll end up at 33. However, since the unit cell is non-orthogonal, # we'll end up at a distorted place. assert_allclose(transformed_coordinates[::3], np.outer(np.arange(33.5),np.array([0.5, 1, np.sqrt(3)/2]))) + +def test_nojump_constantvel(nojump_constantvel_universe): + """ + Test if the nojump transform is returning the correct + values when iterating forwards over the sample trajectory. + """ + ref, towrap = nojump_constantvel_universe + towrap.add_transformations(workflow = [mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump()]) + assert_allclose(towrap.trajectory.ts.positions, ref.trajectory.ts.positions) + From 1ee49479839c5f885fbfdf68a91d05aa0c2fd37b Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 20 Feb 2023 23:52:39 -0800 Subject: [PATCH 21/54] Actually provide a dimension --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index a3e269a8dfc..07ce5438cbd 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -79,6 +79,7 @@ def test_nojump_constantvel(nojump_constantvel_universe): values when iterating forwards over the sample trajectory. """ ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) towrap.add_transformations(workflow = [mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump()]) assert_allclose(towrap.trajectory.ts.positions, ref.trajectory.ts.positions) From e5a10e63875d9a6b87e12b34d1c022942682a2dd Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 00:04:23 -0800 Subject: [PATCH 22/54] Trajectories can be transformed, universes cannot be --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 07ce5438cbd..5f64fb19ddc 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -66,7 +66,7 @@ def test_nojump_nonorthogonal_fwd(nojump_universe): u = nojump_universe dim = np.asarray([1, 1, 1, 90, 60, 90], np.float32) workflow = [mda.transformations.boxdimensions.set_dimensions(dim), NoJump()] - u.add_transformations(*workflow) + u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. However, since the unit cell is non-orthogonal, @@ -80,6 +80,6 @@ def test_nojump_constantvel(nojump_constantvel_universe): """ ref, towrap = nojump_constantvel_universe dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - towrap.add_transformations(workflow = [mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump()]) - assert_allclose(towrap.trajectory.ts.positions, ref.trajectory.ts.positions) + towrap.trajectory.add_transformations(workflow = [mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump()]) + assert_allclose(towrap.trajectory.timeseries(), ref.trajectory.timeseries()) From 112cb7d94089168afd35ba6032a2609c93e11b21 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 06:57:02 -0800 Subject: [PATCH 23/54] Hopefully the last issue? --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 5f64fb19ddc..708b2c1f380 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -51,7 +51,7 @@ def test_nojump_orthogonal_fwd(nojump_universe): u = nojump_universe dim = np.asarray([1, 1, 1, 90, 90, 90], np.float32) workflow = [mda.transformations.boxdimensions.set_dimensions(dim), NoJump()] - u.add_transformations(*workflow) + u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. From 43f3d0cc47ed1c9b7f51a3c7bbcc2509636496ca Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 07:03:39 -0800 Subject: [PATCH 24/54] Fix variables --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 708b2c1f380..9a8de7b352d 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -31,7 +31,7 @@ def nojump_constantvel_universe(): Natom = 1 Nframe = 100 coordinates = np.empty((Nframe, # number of frames - u.atoms.n_atoms, + Natom, 3)) coordinates[:,0,0] = np.linspace(0,45,Nframe) coordinates[:,0,1] = np.linspace(-19,15,Nframe)[::-1] From 232be2ae5b1b69bb7eacead188a38daa6ef9ac52 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 07:13:37 -0800 Subject: [PATCH 25/54] add wrap to the test --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 9a8de7b352d..e42062f6aab 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -3,8 +3,7 @@ from numpy.testing import assert_allclose import MDAnalysis as mda -from MDAnalysis.transformations import NoJump -from MDAnalysisTests import datafiles as data +from MDAnalysis.transformations import NoJump, wrap @pytest.fixture() From f3e18a27bae236519f41ec09287a5c5812421c65 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 10:21:36 -0500 Subject: [PATCH 26/54] Ok, I think I was supposed to use Black to do the code reformatting to make the linter happy --- package/MDAnalysis/transformations/nojump.py | 22 +++++---- .../transformations/test_nojump.py | 45 +++++++++++-------- 2 files changed, 40 insertions(+), 27 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 4d73272f263..3a578b4a713 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -45,7 +45,7 @@ class NoJump(TransformationBase): """ Returns transformed coordinates for the given timestep so that an atom does not move more than half the periodic box size, and will move - across periodic boundary edges. The algorithm used is based on :cite:p:`Kulke2022`, + across periodic boundary edges. The algorithm used is based on :cite:p:`Kulke2022`, equation B6 for non-orthogonal systems. Example @@ -78,10 +78,12 @@ class NoJump(TransformationBase): """ - @due.dcite(Doi("10.1021/acs.jctc.2c00327"), - description="Works through the orthogonal case for unwrapping, " - "and proposes the non-orthogonal approach.", path=__name__) - + @due.dcite( + Doi("10.1021/acs.jctc.2c00327"), + description="Works through the orthogonal case for unwrapping, " + "and proposes the non-orthogonal approach.", + path=__name__, + ) def __init__( self, check_continuity=True, @@ -93,7 +95,7 @@ def __init__( super().__init__(max_threads=max_threads, parallelizable=parallelizable) self.prev = None self.old_frame = 0 - self.older_frame = 'A' + self.older_frame = "A" self.check_c = check_continuity def _transform(self, ts): @@ -101,7 +103,11 @@ def _transform(self, ts): self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame return ts - if self.check_c and self.older_frame != 'A' and (self.old_frame - self.older_frame) != (ts.frame - self.old_frame): + if ( + self.check_c + and self.older_frame != "A" + and (self.old_frame - self.older_frame) != (ts.frame - self.old_frame) + ): warnings.warn( "NoJump detected that the interval between frames is unequal." "We are resetting the reference frame to the current timestep.", @@ -109,7 +115,7 @@ def _transform(self, ts): ) self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) self.old_frame = ts.frame - self.older_frame = 'A' + self.older_frame = "A" return ts if self.check_c and np.abs(self.old_frame - ts.frame) != 1: warnings.warn( diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index e42062f6aab..ffe284fd3fc 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -12,13 +12,11 @@ def nojump_universe(): Create the universe objects for the tests. """ u = mda.Universe.empty(1, trajectory=True) - coordinates = np.empty((100, # number of frames - u.atoms.n_atoms, - 3)) - coordinates[::3,0] = 0 * np.ones(3) / 3 - coordinates[1::3,0] = 1 * np.ones(3) / 3 - coordinates[2::3,0] = 2 * np.ones(3) / 3 - u.load_new(coordinates, order='fac') + coordinates = np.empty((100, u.atoms.n_atoms, 3)) # number of frames + coordinates[::3, 0] = 0 * np.ones(3) / 3 + coordinates[1::3, 0] = 1 * np.ones(3) / 3 + coordinates[2::3, 0] = 2 * np.ones(3) / 3 + u.load_new(coordinates, order="fac") return u @@ -29,16 +27,14 @@ def nojump_constantvel_universe(): """ Natom = 1 Nframe = 100 - coordinates = np.empty((Nframe, # number of frames - Natom, - 3)) - coordinates[:,0,0] = np.linspace(0,45,Nframe) - coordinates[:,0,1] = np.linspace(-19,15,Nframe)[::-1] - coordinates[:,0,2] = np.linspace(-10,10,Nframe) + coordinates = np.empty((Nframe, Natom, 3)) # number of frames + coordinates[:, 0, 0] = np.linspace(0, 45, Nframe) + coordinates[:, 0, 1] = np.linspace(-19, 15, Nframe)[::-1] + coordinates[:, 0, 2] = np.linspace(-10, 10, Nframe) reference = mda.Universe.empty(Natom, trajectory=True) - reference.load_new(coordinates, order='fac') + reference.load_new(coordinates, order="fac") towrap = mda.Universe.empty(Natom, trajectory=True) - towrap.load_new(coordinates, order='fac') + towrap.load_new(coordinates, order="fac") return reference, towrap @@ -54,7 +50,9 @@ def test_nojump_orthogonal_fwd(nojump_universe): transformed_coordinates = u.trajectory.timeseries()[0] # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. - assert_allclose(transformed_coordinates, np.outer(np.linspace(0,33,100),np.ones(3))) + assert_allclose( + transformed_coordinates, np.outer(np.linspace(0, 33, 100), np.ones(3)) + ) def test_nojump_nonorthogonal_fwd(nojump_universe): @@ -70,7 +68,11 @@ def test_nojump_nonorthogonal_fwd(nojump_universe): # Step is 1 unit every 3 steps. After 99 steps from the origin, # we'll end up at 33. However, since the unit cell is non-orthogonal, # we'll end up at a distorted place. - assert_allclose(transformed_coordinates[::3], np.outer(np.arange(33.5),np.array([0.5, 1, np.sqrt(3)/2]))) + assert_allclose( + transformed_coordinates[::3], + np.outer(np.arange(33.5), np.array([0.5, 1, np.sqrt(3) / 2])), + ) + def test_nojump_constantvel(nojump_constantvel_universe): """ @@ -79,6 +81,11 @@ def test_nojump_constantvel(nojump_constantvel_universe): """ ref, towrap = nojump_constantvel_universe dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - towrap.trajectory.add_transformations(workflow = [mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump()]) + towrap.trajectory.add_transformations( + workflow=[ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + ) assert_allclose(towrap.trajectory.timeseries(), ref.trajectory.timeseries()) - From 79927ca89b9df58fdb3b6d78857af3a6da8eff3a Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 07:33:21 -0800 Subject: [PATCH 27/54] Black did something dumb I think --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index ffe284fd3fc..d7979780299 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -81,11 +81,10 @@ def test_nojump_constantvel(nojump_constantvel_universe): """ ref, towrap = nojump_constantvel_universe dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - towrap.trajectory.add_transformations( - workflow=[ + workflow=[ mda.transformations.boxdimensions.set_dimensions(dim), wrap(towrap.atoms), NoJump(), ] - ) + towrap.trajectory.add_transformations(*workflow) assert_allclose(towrap.trajectory.timeseries(), ref.trajectory.timeseries()) From 08477ee428706ffcf0b4935a8571ab84a18a4492 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 07:58:32 -0800 Subject: [PATCH 28/54] I forgot that wrapping in negative space can lead to problems --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index d7979780299..4e3f34123ed 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -28,9 +28,9 @@ def nojump_constantvel_universe(): Natom = 1 Nframe = 100 coordinates = np.empty((Nframe, Natom, 3)) # number of frames - coordinates[:, 0, 0] = np.linspace(0, 45, Nframe) - coordinates[:, 0, 1] = np.linspace(-19, 15, Nframe)[::-1] - coordinates[:, 0, 2] = np.linspace(-10, 10, Nframe) + coordinates[:, 0, 0] = np.linspace(0,45,Nframe) + coordinates[:, 0, 1] = np.linspace(0,15,Nframe) + coordinates[:, 0, 2] = np.linspace(0,10,Nframe) reference = mda.Universe.empty(Natom, trajectory=True) reference.load_new(coordinates, order="fac") towrap = mda.Universe.empty(Natom, trajectory=True) From 64a4b1dee9a476097d6618466c3646cd86b03993 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 08:17:05 -0800 Subject: [PATCH 29/54] Bump up the tolerance a bit. --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 4e3f34123ed..40db2b287d1 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -87,4 +87,9 @@ def test_nojump_constantvel(nojump_constantvel_universe): NoJump(), ] towrap.trajectory.add_transformations(*workflow) - assert_allclose(towrap.trajectory.timeseries(), ref.trajectory.timeseries()) + assert_allclose( + towrap.trajectory.timeseries(), + ref.trajectory.timeseries(), + rtol=5e-07, + atol=5e-06, + ) From f1b3588bea9d7937a79a43dcb210b4a11a913139 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 21 Feb 2023 09:32:28 -0800 Subject: [PATCH 30/54] Testing lint fix --- .../transformations/test_nojump.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 40db2b287d1..ee495700d8c 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -28,9 +28,9 @@ def nojump_constantvel_universe(): Natom = 1 Nframe = 100 coordinates = np.empty((Nframe, Natom, 3)) # number of frames - coordinates[:, 0, 0] = np.linspace(0,45,Nframe) - coordinates[:, 0, 1] = np.linspace(0,15,Nframe) - coordinates[:, 0, 2] = np.linspace(0,10,Nframe) + coordinates[:, 0, 0] = np.linspace(0, 45, Nframe) + coordinates[:, 0, 1] = np.linspace(0, 15, Nframe) + coordinates[:, 0, 2] = np.linspace(0, 10, Nframe) reference = mda.Universe.empty(Natom, trajectory=True) reference.load_new(coordinates, order="fac") towrap = mda.Universe.empty(Natom, trajectory=True) @@ -81,11 +81,11 @@ def test_nojump_constantvel(nojump_constantvel_universe): """ ref, towrap = nojump_constantvel_universe dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow=[ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] towrap.trajectory.add_transformations(*workflow) assert_allclose( towrap.trajectory.timeseries(), From 0eccb0a9a753040436bd289aedf9f1b314abd6d4 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Wed, 22 Feb 2023 21:20:28 -0800 Subject: [PATCH 31/54] Added some more tests to cover the branches. --- .../transformations/test_nojump.py | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index ee495700d8c..52a08cafde7 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -93,3 +93,55 @@ def test_nojump_constantvel(nojump_constantvel_universe): rtol=5e-07, atol=5e-06, ) + + +def test_nojump_constantvel_skip(nojump_constantvel_universe): + """ + Test if the nojump transform is returning the correct + values when iterating forwards over the sample trajectory, + skipping by 2. + """ + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[::2]], + [ts for ts in towrap.trajectory[::2]], + ): + assert_allclose( + r.trajectory.ts.positions, + c.trajectory.ts.positions, + rtol=5e-07, + atol=5e-06, + ) + + +def test_nojump_constantvel_jumparound(nojump_constantvel_universe): + """ + Test if the nojump transform is returning the correct + values when iterating forwards over the sample trajectory, + skipping by 2. + """ + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[0,1,3,5,4]], + [ts for ts in towrap.trajectory[0,1,3,5,4]], + ): + assert_allclose( + r.trajectory.ts.positions, + c.trajectory.ts.positions, + rtol=5e-07, + atol=5e-06, + ) From 42c61121e2a14d9db5e6cabd2e0e60393453cbfa Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Wed, 22 Feb 2023 21:52:12 -0800 Subject: [PATCH 32/54] Fix the tests again --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 52a08cafde7..3714eff955f 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -114,8 +114,8 @@ def test_nojump_constantvel_skip(nojump_constantvel_universe): [ts for ts in towrap.trajectory[::2]], ): assert_allclose( - r.trajectory.ts.positions, - c.trajectory.ts.positions, + r.positions, + c.positions, rtol=5e-07, atol=5e-06, ) @@ -140,8 +140,8 @@ def test_nojump_constantvel_jumparound(nojump_constantvel_universe): [ts for ts in towrap.trajectory[0,1,3,5,4]], ): assert_allclose( - r.trajectory.ts.positions, - c.trajectory.ts.positions, + r.positions, + c.positions, rtol=5e-07, atol=5e-06, ) From f4904c681dee1e030d8bda6b74191d9722d136a9 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Thu, 23 Feb 2023 07:02:47 -0800 Subject: [PATCH 33/54] I'm definitely not used to how MDAnalysis does things. --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 3714eff955f..99882deef6d 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -136,8 +136,8 @@ def test_nojump_constantvel_jumparound(nojump_constantvel_universe): ] towrap.trajectory.add_transformations(*workflow) for r, c in zip( - [ts for ts in ref.trajectory[0,1,3,5,4]], - [ts for ts in towrap.trajectory[0,1,3,5,4]], + [ts for ts in ref.trajectory[[0,1,3,5,4]]], + [ts for ts in towrap.trajectory[[0,1,3,5,4]]], ): assert_allclose( r.positions, From a3b1a2122a7858a608dd2206b9686d2fbfa811d7 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Thu, 23 Feb 2023 07:31:38 -0800 Subject: [PATCH 34/54] Apparently I'm not covering a branch. --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 99882deef6d..b3f2639385e 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -136,8 +136,8 @@ def test_nojump_constantvel_jumparound(nojump_constantvel_universe): ] towrap.trajectory.add_transformations(*workflow) for r, c in zip( - [ts for ts in ref.trajectory[[0,1,3,5,4]]], - [ts for ts in towrap.trajectory[[0,1,3,5,4]]], + [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], + [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], ): assert_allclose( r.positions, From 3c2a90b8b52e673f1a38004d8106eb525fb9c6f0 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Thu, 23 Feb 2023 22:12:50 -0800 Subject: [PATCH 35/54] Make the tests better, also give the user a better error message if there is an illegal box. --- package/MDAnalysis/transformations/nojump.py | 7 +- .../transformations/test_nojump.py | 82 +++++++++++-------- 2 files changed, 52 insertions(+), 37 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 3a578b4a713..5a1c4a3e8ca 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -39,6 +39,7 @@ from .base import TransformationBase from ..due import due, Doi +from ..exceptions import NoDataError class NoJump(TransformationBase): @@ -129,7 +130,11 @@ def _transform(self, ts): # Remember that @ is a shorthand for matrix multiplication. # np.matmul(a, b) is equivalent to a @ b L = ts.triclinic_dimensions - Linverse = np.linalg.inv(L) + try: + Linverse = np.linalg.inv(L) + except np.linalg.LinAlgError: + msg = "Periodic box dimensions are not invertible at step %d." % ts.frame + raise NoDataError(msg) # Convert into reduced coordinate space fcurrent = ts.positions @ Linverse fprev = self.prev # Previous unwrapped coordinates in reduced box coordinates. diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index b3f2639385e..c8c41d649a1 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -101,24 +101,25 @@ def test_nojump_constantvel_skip(nojump_constantvel_universe): values when iterating forwards over the sample trajectory, skipping by 2. """ - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[::2]], - [ts for ts in towrap.trajectory[::2]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) + with pytest.warns(UserWarning): + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[::2]], + [ts for ts in towrap.trajectory[::2]], + ): + assert_allclose( + r.positions, + c.positions, + rtol=5e-07, + atol=5e-06, + ) def test_nojump_constantvel_jumparound(nojump_constantvel_universe): @@ -127,21 +128,30 @@ def test_nojump_constantvel_jumparound(nojump_constantvel_universe): values when iterating forwards over the sample trajectory, skipping by 2. """ - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], - [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) + with pytest.warns(UserWarning): + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], + [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], + ): + assert_allclose( + r.positions, + c.positions, + rtol=5e-07, + atol=5e-06, + ) + + +def test_missing_dimensions(nojump_universe): + with pytest.raises(mda.exceptions.NoDataError): + u = nojump_universe + workflow = [NoJump()] + u.trajectory.add_transformations(*workflow) + transformed_coordinates = u.trajectory.timeseries()[0] From 22953e04b138f11112109cb73e046d7ed899c568 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 06:10:16 -0800 Subject: [PATCH 36/54] I CANNOT trigger the warnings. --- .../transformations/test_nojump.py | 74 +++++++++---------- 1 file changed, 36 insertions(+), 38 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index c8c41d649a1..6e63a32a622 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -101,25 +101,24 @@ def test_nojump_constantvel_skip(nojump_constantvel_universe): values when iterating forwards over the sample trajectory, skipping by 2. """ - with pytest.warns(UserWarning): - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[::2]], - [ts for ts in towrap.trajectory[::2]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[::2]], + [ts for ts in towrap.trajectory[::2]], + ): + assert_allclose( + r.positions, + c.positions, + rtol=5e-07, + atol=5e-06, + ) def test_nojump_constantvel_jumparound(nojump_constantvel_universe): @@ -128,25 +127,24 @@ def test_nojump_constantvel_jumparound(nojump_constantvel_universe): values when iterating forwards over the sample trajectory, skipping by 2. """ - with pytest.warns(UserWarning): - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], - [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) + ref, towrap = nojump_constantvel_universe + dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) + workflow = [ + mda.transformations.boxdimensions.set_dimensions(dim), + wrap(towrap.atoms), + NoJump(), + ] + towrap.trajectory.add_transformations(*workflow) + for r, c in zip( + [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], + [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], + ): + assert_allclose( + r.positions, + c.positions, + rtol=5e-07, + atol=5e-06, + ) def test_missing_dimensions(nojump_universe): From faf8db91ac38ed0c9c07120f5928727c5ca86fbe Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 07:54:52 -0800 Subject: [PATCH 37/54] I am now testing the warnings via a loaded file path --- .../transformations/test_nojump.py | 72 +++++++------------ 1 file changed, 26 insertions(+), 46 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 6e63a32a622..d8109e263fa 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -6,6 +6,17 @@ from MDAnalysis.transformations import NoJump, wrap +@pytest.fixture() +def nojump_universes_fromfile(): + ''' + Create the universe objects for the tests. + ''' + u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + transformation = NoJump() + u.trajectory.add_transformations(transformation) + return u + + @pytest.fixture() def nojump_universe(): """ @@ -95,56 +106,25 @@ def test_nojump_constantvel(nojump_constantvel_universe): ) -def test_nojump_constantvel_skip(nojump_constantvel_universe): +def test_nojump_constantvel_skip(nojump_universes_fromfile): """ - Test if the nojump transform is returning the correct - values when iterating forwards over the sample trajectory, - skipping by 2. + Test if the nojump transform warning is emitted. """ - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[::2]], - [ts for ts in towrap.trajectory[::2]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) - - -def test_nojump_constantvel_jumparound(nojump_constantvel_universe): + with pytest.warns(UserWarning): + u = nojump_universes_fromfile + u.trajectory[0] + u.trajectory[99] #Exercises the warning. + + +def test_nojump_constantvel_jumparound(nojump_universes_fromfile): """ - Test if the nojump transform is returning the correct - values when iterating forwards over the sample trajectory, - skipping by 2. + Test if the nojump transform is emitting a warning. """ - ref, towrap = nojump_constantvel_universe - dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) - workflow = [ - mda.transformations.boxdimensions.set_dimensions(dim), - wrap(towrap.atoms), - NoJump(), - ] - towrap.trajectory.add_transformations(*workflow) - for r, c in zip( - [ts for ts in ref.trajectory[[0,1,2,3,5,4]]], - [ts for ts in towrap.trajectory[[0,1,2,3,5,4]]], - ): - assert_allclose( - r.positions, - c.positions, - rtol=5e-07, - atol=5e-06, - ) + with pytest.warns(UserWarning): + u = nojump_universes_fromfile + u.trajectory[0] + u.trajectory[1] + u.trajectory[99] def test_missing_dimensions(nojump_universe): From b8afc9699feaa7c4546fe9312b99f2375a0d7f50 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 13:46:46 -0800 Subject: [PATCH 38/54] MDAnalysis->mda: --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index d8109e263fa..a002780bd79 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -11,7 +11,7 @@ def nojump_universes_fromfile(): ''' Create the universe objects for the tests. ''' - u = MDAnalysis.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) + u = mda.Universe(data.PSF_TRICLINIC, data.DCD_TRICLINIC) transformation = NoJump() u.trajectory.add_transformations(transformation) return u From 903cb290d34e1fc58909d2322a846c304922b358 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 14:08:23 -0800 Subject: [PATCH 39/54] Ok, misplaced an import --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index a002780bd79..d168d3de050 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -4,6 +4,7 @@ import MDAnalysis as mda from MDAnalysis.transformations import NoJump, wrap +from MDAnalysisTests import datafiles as data @pytest.fixture() From 535132e5992ad8aa8652a61ffc21a6bc0b00273b Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 14:18:02 -0800 Subject: [PATCH 40/54] No longer out of bounds --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index d168d3de050..aa0b19a9611 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -125,7 +125,7 @@ def test_nojump_constantvel_jumparound(nojump_universes_fromfile): u = nojump_universes_fromfile u.trajectory[0] u.trajectory[1] - u.trajectory[99] + u.trajectory[9] def test_missing_dimensions(nojump_universe): From e0e2b1878383ec59446d409f25439e6300c06650 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 14:20:38 -0800 Subject: [PATCH 41/54] Catch busted Linverse during initialization. --- package/MDAnalysis/transformations/nojump.py | 7 ++++++- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 5a1c4a3e8ca..5fb92fb3298 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -101,7 +101,12 @@ def __init__( def _transform(self, ts): if self.prev is None: - self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) + try: + Linverse = np.linalg.inv(ts.triclinic_dimensions) + except np.linalg.LinAlgError: + msg = "Periodic box dimensions are not invertible at step %d." % ts.frame + raise NoDataError(msg) + self.prev = ts.positions @ Linverse self.old_frame = ts.frame return ts if ( diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index aa0b19a9611..65bf3fd680e 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -114,7 +114,7 @@ def test_nojump_constantvel_skip(nojump_universes_fromfile): with pytest.warns(UserWarning): u = nojump_universes_fromfile u.trajectory[0] - u.trajectory[99] #Exercises the warning. + u.trajectory[9] #Exercises the warning. def test_nojump_constantvel_jumparound(nojump_universes_fromfile): From eec1758ca4eb90189148027ca4f8131f3a8818a6 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 14:29:31 -0800 Subject: [PATCH 42/54] Add an initial dimension, which hopefully gets ignored for subsequent frames --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 65bf3fd680e..f485df5ad78 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -128,9 +128,18 @@ def test_nojump_constantvel_jumparound(nojump_universes_fromfile): u.trajectory[9] +def test_missing_dimensions_init(nojump_universe): + with pytest.raises(mda.exceptions.NoDataError): + u = nojump_universe + workflow = [NoJump()] + u.trajectory.add_transformations(*workflow) + transformed_coordinates = u.trajectory.timeseries()[0] + + def test_missing_dimensions(nojump_universe): with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe + u.dimensions = [73, 73, 73, 90, 90, 90] workflow = [NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] From 0e6431283614934b52e006c9086ec6e0973244e1 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Fri, 24 Feb 2023 14:32:54 -0800 Subject: [PATCH 43/54] Add text description for what the intended difference between the two missing dimension tests is. --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index f485df5ad78..e378ce1dc22 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -129,6 +129,10 @@ def test_nojump_constantvel_jumparound(nojump_universes_fromfile): def test_missing_dimensions_init(nojump_universe): + """ + Test if the nojump transform raises a NoDataError if there is no + initial dimension for the periodic unit cell. + """ with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe workflow = [NoJump()] @@ -137,6 +141,10 @@ def test_missing_dimensions_init(nojump_universe): def test_missing_dimensions(nojump_universe): + """ + Test if the nojump transform raises a NoDataError if there is no + dimension for the periodic unit cell in a subsequent timestep. + """ with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe u.dimensions = [73, 73, 73, 90, 90, 90] From be901c3f2f4196776e61ca0f6d29dada6fce266d Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 26 Feb 2023 14:47:47 -0500 Subject: [PATCH 44/54] Separating no dimensions from invalid dimensions tests. --- package/MDAnalysis/transformations/nojump.py | 24 +++++++-------- .../transformations/test_nojump.py | 30 +++++++++++++++++-- 2 files changed, 37 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 5fb92fb3298..337653aacfd 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -100,12 +100,16 @@ def __init__( self.check_c = check_continuity def _transform(self, ts): + L = ts.triclinic_dimensions + if L is None: + msg = "Periodic box dimensions not provided at step %d." % ts.frame + raise NoDataError(msg) + try: + Linverse = np.linalg.inv(L) + except np.linalg.LinAlgError: + msg = "Periodic box dimensions are not invertible at step %d." % ts.frame + raise NoDataError(msg) if self.prev is None: - try: - Linverse = np.linalg.inv(ts.triclinic_dimensions) - except np.linalg.LinAlgError: - msg = "Periodic box dimensions are not invertible at step %d." % ts.frame - raise NoDataError(msg) self.prev = ts.positions @ Linverse self.old_frame = ts.frame return ts @@ -119,7 +123,7 @@ def _transform(self, ts): "We are resetting the reference frame to the current timestep.", UserWarning, ) - self.prev = ts.positions @ np.linalg.inv(ts.triclinic_dimensions) + self.prev = ts.positions @ Linverse self.old_frame = ts.frame self.older_frame = "A" return ts @@ -132,14 +136,6 @@ def _transform(self, ts): "this might be inaccurate.", UserWarning, ) - # Remember that @ is a shorthand for matrix multiplication. - # np.matmul(a, b) is equivalent to a @ b - L = ts.triclinic_dimensions - try: - Linverse = np.linalg.inv(L) - except np.linalg.LinAlgError: - msg = "Periodic box dimensions are not invertible at step %d." % ts.frame - raise NoDataError(msg) # Convert into reduced coordinate space fcurrent = ts.positions @ Linverse fprev = self.prev # Previous unwrapped coordinates in reduced box coordinates. diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index f485df5ad78..6464f0468ee 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -73,17 +73,32 @@ def test_nojump_nonorthogonal_fwd(nojump_universe): values when iterating forwards over the sample trajectory. """ u = nojump_universe + # Set a non-orthogonal box dimension. The box below works out to be this cell: + # [[1. 0. 0. ] + # [0. 1. 0. ] + # [0.5 0. 0.8660254]] dim = np.asarray([1, 1, 1, 90, 60, 90], np.float32) workflow = [mda.transformations.boxdimensions.set_dimensions(dim), NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] - # Step is 1 unit every 3 steps. After 99 steps from the origin, - # we'll end up at 33. However, since the unit cell is non-orthogonal, - # we'll end up at a distorted place. + # After the transformation, you should end up in a repeating pattern, since you are + # working in a hexagonal unit cell system. Since you jump every third timestep across + # a periodic boundary, the shift in each axis is saved. As a consequence, the correct + # jump every third step is just the original position + the size of the periodic cells. assert_allclose( transformed_coordinates[::3], np.outer(np.arange(33.5), np.array([0.5, 1, np.sqrt(3) / 2])), ) + assert_allclose( + transformed_coordinates[1::3], + np.outer(np.arange(32.5), np.array([0.5, 1, np.sqrt(3) / 2])) + 1 * np.ones(3) / 3, + rtol=1.2e-7 + ) + assert_allclose( + transformed_coordinates[2::3], + np.outer(np.arange(32.5), np.array([0.5, 1, np.sqrt(3) / 2])) + 2 * np.ones(3) / 3, + rtol=1.2e-7 + ) def test_nojump_constantvel(nojump_constantvel_universe): @@ -143,3 +158,12 @@ def test_missing_dimensions(nojump_universe): workflow = [NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] + + +def test_notinvertible(nojump_universe): + with pytest.raises(mda.exceptions.NoDataError): + u = nojump_universe + u.dimensions = [0, 0, 0, 90, 90, 90] + workflow = [NoJump()] + u.trajectory.add_transformations(*workflow) + transformed_coordinates = u.trajectory.timeseries()[0] From 54b5651865845d79f48e373d5bf9d31dd6f2d368 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 26 Feb 2023 14:49:25 -0500 Subject: [PATCH 45/54] Add some text explaining the notinvertible transform --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index 2e6a67b3a78..fd3b64314d8 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -169,6 +169,10 @@ def test_missing_dimensions(nojump_universe): def test_notinvertible(nojump_universe): + """ + Test if the nojump transform raises a NoDataError if the dimensions + are invalid for the periodic unit cell. + """ with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe u.dimensions = [0, 0, 0, 90, 90, 90] From a484becfeae084ae7f995e44959f1c1662e99111 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 26 Feb 2023 15:50:14 -0500 Subject: [PATCH 46/54] Apply zero dimensions --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index fd3b64314d8..f0924b39432 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -176,6 +176,6 @@ def test_notinvertible(nojump_universe): with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe u.dimensions = [0, 0, 0, 90, 90, 90] - workflow = [NoJump()] + workflow = [mda.transformations.boxdimensions.set_dimensions(dim),NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] From 3ffdc27d0c9a8bebdae95a70abaa19543c17406e Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sat, 4 Mar 2023 21:46:35 -0500 Subject: [PATCH 47/54] Whoops. Actually make something non-invertible --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index f0924b39432..e1615b744b2 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -175,7 +175,7 @@ def test_notinvertible(nojump_universe): """ with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe - u.dimensions = [0, 0, 0, 90, 90, 90] + dim = [0, 0, 0, 90, 90, 90] workflow = [mda.transformations.boxdimensions.set_dimensions(dim),NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] From 47b0e39f121e02308d36493e383fae62dbad7536 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Wed, 15 Mar 2023 20:07:08 -0400 Subject: [PATCH 48/54] Update __init__.py Darker complains, so I guess the import has to go. --- package/MDAnalysis/transformations/__init__.py | 1 - 1 file changed, 1 deletion(-) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index f359363157d..e9508cddb35 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -129,7 +129,6 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager -from .nojump import NoJump from .fit import fit_translation, fit_rot_trans from .wrap import wrap, unwrap from .boxdimensions import set_dimensions From 89d9eb77f704583a3002529c3ff6d058c6a117c0 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Wed, 15 Mar 2023 20:16:31 -0400 Subject: [PATCH 49/54] Revert "Update __init__.py" This reverts commit 47b0e39f121e02308d36493e383fae62dbad7536. --- package/MDAnalysis/transformations/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/package/MDAnalysis/transformations/__init__.py b/package/MDAnalysis/transformations/__init__.py index e9508cddb35..f359363157d 100644 --- a/package/MDAnalysis/transformations/__init__.py +++ b/package/MDAnalysis/transformations/__init__.py @@ -129,6 +129,7 @@ def wrapped(ts): from .translate import translate, center_in_box from .rotate import rotateby from .positionaveraging import PositionAverager +from .nojump import NoJump from .fit import fit_translation, fit_rot_trans from .wrap import wrap, unwrap from .boxdimensions import set_dimensions From 997b221246e038b428f5c5b8d3bac185228caffa Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Sun, 26 Mar 2023 06:07:17 -0400 Subject: [PATCH 50/54] Need to pass at least one non-zero dimension to trigger a LinAlgError --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index e1615b744b2..bd6332c03d5 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -175,7 +175,7 @@ def test_notinvertible(nojump_universe): """ with pytest.raises(mda.exceptions.NoDataError): u = nojump_universe - dim = [0, 0, 0, 90, 90, 90] + dim = [1, 0, 0, 90, 90, 90] workflow = [mda.transformations.boxdimensions.set_dimensions(dim),NoJump()] u.trajectory.add_transformations(*workflow) transformed_coordinates = u.trajectory.timeseries()[0] From 70d23ba8af4d0e07b2d35006abfe31e7811caacb Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Mon, 27 Mar 2023 09:58:09 -0400 Subject: [PATCH 51/54] Used mdanalysis.copy inside a test --- testsuite/MDAnalysisTests/transformations/test_nojump.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/testsuite/MDAnalysisTests/transformations/test_nojump.py b/testsuite/MDAnalysisTests/transformations/test_nojump.py index bd6332c03d5..4352c9c4bdc 100644 --- a/testsuite/MDAnalysisTests/transformations/test_nojump.py +++ b/testsuite/MDAnalysisTests/transformations/test_nojump.py @@ -45,9 +45,7 @@ def nojump_constantvel_universe(): coordinates[:, 0, 2] = np.linspace(0, 10, Nframe) reference = mda.Universe.empty(Natom, trajectory=True) reference.load_new(coordinates, order="fac") - towrap = mda.Universe.empty(Natom, trajectory=True) - towrap.load_new(coordinates, order="fac") - return reference, towrap + return reference def test_nojump_orthogonal_fwd(nojump_universe): @@ -106,7 +104,9 @@ def test_nojump_constantvel(nojump_constantvel_universe): Test if the nojump transform is returning the correct values when iterating forwards over the sample trajectory. """ - ref, towrap = nojump_constantvel_universe + ref = nojump_constantvel_universe + towrap = ref.copy() # This copy of the universe will be wrapped, then unwrapped, + # and should be equal to ref. dim = np.asarray([5, 5, 5, 54, 60, 90], np.float32) workflow = [ mda.transformations.boxdimensions.set_dimensions(dim), From 7ae4e29ca83e73b7a2cc1b911e57b1420efdc7a2 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 28 Mar 2023 14:58:57 -0400 Subject: [PATCH 52/54] Update package/MDAnalysis/transformations/nojump.py Co-authored-by: Oliver Beckstein --- package/MDAnalysis/transformations/nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index 337653aacfd..c37081846ec 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -102,7 +102,7 @@ def __init__( def _transform(self, ts): L = ts.triclinic_dimensions if L is None: - msg = "Periodic box dimensions not provided at step %d." % ts.frame + msg = f"Periodic box dimensions not provided at step {ts.frame}" raise NoDataError(msg) try: Linverse = np.linalg.inv(L) From 49f5c013541b09e568ba02dbd87768188cffe93c Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 28 Mar 2023 14:59:29 -0400 Subject: [PATCH 53/54] Update package/MDAnalysis/transformations/nojump.py Co-authored-by: Oliver Beckstein --- package/MDAnalysis/transformations/nojump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index c37081846ec..ef9f2383a80 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -107,7 +107,7 @@ def _transform(self, ts): try: Linverse = np.linalg.inv(L) except np.linalg.LinAlgError: - msg = "Periodic box dimensions are not invertible at step %d." % ts.frame + msg = f"Periodic box dimensions are not invertible at step {ts.frame}" raise NoDataError(msg) if self.prev is None: self.prev = ts.positions @ Linverse From 95f08a45faeb09593f3fffa46924e18366814364 Mon Sep 17 00:00:00 2001 From: Josh Vermaas Date: Tue, 28 Mar 2023 16:01:25 -0400 Subject: [PATCH 54/54] Update nojump.py Added a bit more to the docstring that highlights what can go wrong when applying the transform to a trajectory that is missing box dimensions. --- package/MDAnalysis/transformations/nojump.py | 26 +++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/package/MDAnalysis/transformations/nojump.py b/package/MDAnalysis/transformations/nojump.py index ef9f2383a80..adfa24ea7b7 100644 --- a/package/MDAnalysis/transformations/nojump.py +++ b/package/MDAnalysis/transformations/nojump.py @@ -28,7 +28,8 @@ Unwraps the trajectory such that atoms never move more than half a periodic box length. The consequence of this is that particles will diffuse across periodic boundaries when needed. This unwrapping method is suitable as a preprocessing step to calculate -molecular diffusion. +molecular diffusion, or more commonly to keep multi-domain proteins whole during trajectory +analysis. The algorithm used is based on :cite:p:`Kulke2022`. .. autoclass:: NoJump @@ -45,15 +46,26 @@ class NoJump(TransformationBase): """ Returns transformed coordinates for the given timestep so that an atom - does not move more than half the periodic box size, and will move + does not move more than half the periodic box size between two timesteps, and will move across periodic boundary edges. The algorithm used is based on :cite:p:`Kulke2022`, - equation B6 for non-orthogonal systems. + equation B6 for non-orthogonal systems, so it is general to most applications where + molecule trajectories should not "jump" from one side of a periodic box to another. + + Note that this transformation depends on a periodic box dimension being set for every + frame in the trajectory, and that this box dimension can be transformed to an orthonormal + unit cell. If not, an error is emitted. Since it is typical to transform all frames + sequentially when unwrapping trajectories, warnings are emitted if non-sequential timesteps + are transformed using NoJump. Example ------- - To calculate diffusion, one of the first steps is to compute a trajectory - where the atoms do not cross the periodic boundary. + Suppose you had a molecular trajectory with a protein dimer that moved across a periodic + boundary. This will normally appear to make the dimer split apart. This transformation + uses the molecular motions between frames in a wrapped trajectory to create an unwrapped + trajectory where molecules never move more than half a periodic box dimension between subsequent + frames. Thus, the normal use case is to apply the transformation to every frame of a trajectory, + and to do so sequentially. .. code-block:: python @@ -62,7 +74,9 @@ class NoJump(TransformationBase): for ts in u.trajectory: print(ts.positions) - In this case, ``ts.positions`` will return the NoJump unwrapped trajectory. + In this case, ``ts.positions`` will return the NoJump unwrapped trajectory, which would keep + protein dimers together if they are together in the inital timestep. The unwrapped trajectory may + also be useful to compute diffusion coefficients via the Einstein relation. To reverse the process, wrap the coordinates again. Returns