diff --git a/.appveyor.yml b/.appveyor.yml index b3a24c9b0f7..dfb227c372d 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -10,7 +10,7 @@ cache: environment: global: CONDA_CHANNELS: conda-forge - CONDA_DEPENDENCIES: pip setuptools wheel cython mock six biopython networkx joblib matplotlib scipy vs2015_runtime pytest mmtf-python GridDataFormats hypothesis pytest-cov codecov chemfiles tqdm + CONDA_DEPENDENCIES: pip setuptools wheel cython mock six biopython networkx joblib matplotlib scipy vs2015_runtime pytest mmtf-python GridDataFormats hypothesis pytest-cov codecov chemfiles tqdm tidynamics>=1.0.0 PIP_DEPENDENCIES: gsd==1.9.3 duecredit parmed DEBUG: "False" MINGW_64: C:\mingw-w64\x86_64-6.3.0-posix-seh-rt_v5-rev1\mingw64\bin diff --git a/.travis.yml b/.travis.yml index 75816f20aad..ba98443157a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ env: - SETUP_CMD="${PYTEST_FLAGS}" - BUILD_CMD="pip install -e package/ && (cd testsuite/ && python setup.py build)" - CONDA_MIN_DEPENDENCIES="mmtf-python mock six biopython networkx cython matplotlib scipy griddataformats hypothesis gsd codecov" - - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles tqdm>=4.43.0" + - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles tqdm>=4.43.0 tidynamics>=1.0.0" - CONDA_CHANNELS='biobuilds conda-forge' - CONDA_CHANNEL_PRIORITY=True - PIP_DEPENDENCIES="duecredit parmed" diff --git a/package/CHANGELOG b/package/CHANGELOG index 419421b576c..27bc15e460e 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -13,7 +13,7 @@ The rules for this file: * release numbers follow "Semantic Versioning" http://semver.org ------------------------------------------------------------------------------ -??/??/?? richardjgowers, IAlibay +??/??/?? richardjgowers, IAlibay, hmacdope * 2.0.0 @@ -22,6 +22,7 @@ Fixes (Issues #2449, #2651) Enhancements + * Added computation of Mean Squared Displacements (#2438, PR #2619) Changes * Removes deprecated density_from_Universe, density_from_PDB, Bfactor2RMSF, diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py new file mode 100644 index 00000000000..96202647271 --- /dev/null +++ b/package/MDAnalysis/analysis/msd.py @@ -0,0 +1,384 @@ +# -*- 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 +# + +r""" +Mean Squared Displacement --- :mod:`MDAnalysis.analysis.msd` +============================================================== + +:Authors: Hugo MacDermott-Opeskin +:Year: 2020 +:Copyright: GNU Public License v2 + +This module implements the calculation of Mean Squared Displacements (MSDs) +by the Einstein relation. MSDs can be used to characterize the speed at +which particles move and has its roots in the study of Brownian motion. +For a full explanation of the theory behind MSDs and the subsequent +calculation of self-diffusivities the reader is directed to [Maginn2019]_. +MSDs can be computed from the following expression, known as the +**Einstein formula**: + +.. math:: + + MSD(r_{d}) = \bigg{\langle} \frac{1}{N} \sum_{i=1}^{N} |r_{d} + - r_{d}(t_0)|^2 \bigg{\rangle}_{t_{0}} + +where :math:`N` is the number of equivalent particles the MSD is calculated +over, :math:`r` are their coordinates and :math:`d` the desired dimensionality +of the MSD. Note that while the definition of the MSD is universal, there are +many practical considerations to computing the MSD that vary between +implementations. In this module, we compute a "windowed" MSD, where the MSD +is averaged over all possible lag-times :math:`\tau \le \tau_{max}`, +where :math:`\tau_{max}` is the length of the trajectory, thereby maximizing +the number of samples. + +The computation of the MSD in this way can be computationally intensive due to +its :math:`N^2` scaling with respect to :math:`\tau_{max}`. An algorithm to +compute the MSD with :math:`N log(N)` scaling based on a Fast Fourier +Transform is known and can be accessed by setting ``fft=True`` [Calandri2011]_ +[Buyl2018]_. The FFT-based approach requires that the +`tidynamics `_ package is +installed; otherwise the code will raise an :exc:`ImportError`. + +Please cite [Calandri2011]_ [Buyl2018]_ if you use this module in addition to +the normal MDAnalysis citations. + +Computing an MSD +---------------- +This example computes a 3D MSD for the movement of 100 particles undergoing a +random walk. Files provided as part of the MDAnalysis test suite are used +(in the variables :data:`~MDAnalysis.tests.datafiles.RANDOM_WALK` and +:data:`~MDAnalysis.tests.datafiles.RANDOM_WALK_TOPO`) + +First load all modules and test data + +.. code-block:: python + + import MDAnalysis as mda + import MDAnalysis.analysis.msd as msd + from MDAnalysis.tests.datafiles import RANDOM_WALK, RANDOM_WALK_TOPO + +Given a universe containing trajectory data we can extract the MSD +analysis by using the class :class:`EinsteinMSD` + +.. code-block:: python + + u = mda.Universe(RANDOM_WALK, RANDOM_WALK_TOPO) + MSD = msd.EinsteinMSD(u, 'all', msd_type='xyz', fft=True) + MSD.run() + +The MSD can then be accessed as + +.. code-block:: python + + msd = MSD.timeseries + +Visual inspection of the MSD is important, so let's take a look at it with a + simple plot. + +.. code-block:: python + + import matplotlib.pyplot as plt + nframes = MSD.n_frames + timestep = 1 # this needs to be the actual time between frames + lagtimes = np.arange(nframes)*timestep # make the lag-time axis + fig = plt.figure() + ax = plt.axes() + # plot the actual MSD + ax.plot(lagtimes, msd, lc="black", ls="-", label=r'3D random walk') + exact = lagtimes*6 + # plot the exact result + ax.plot(lagtimes, exact, lc="black", ls="--", label=r'$y=2 D\tau$') + plt.show() + +This gives us the plot of the MSD with respect to lag-time (:math:`\tau`). +We can see that the MSD is approximately linear with respect to :math:`\tau`. +This is a numerical example of a known theoretical result that the MSD of a +random walk is linear with respect to lag-time, with a slope of :math:`2d`. +In this expression :math:`d` is the dimensionality of the MSD. For our 3D MSD, +this is 3. For comparison we have plotted the line :math:`y=6\tau` to which an +ensemble of 3D random walks should converge. + +.. _figure-msd: + +.. figure:: /images/msd_demo_plot.png + :scale: 100 % + :alt: MSD plot + +Note that a segment of the MSD is required to be linear to accurately +determine self-diffusivity. This linear segment represents the so called +"middle" of the MSD plot, where ballistic trajectories at short time-lags are +excluded along with poorly averaged data at long time-lags. We can select the +"middle" of the MSD by indexing the MSD and the time-lags. Appropriately +linear segments of the MSD can be confirmed with a log-log plot as is often +reccomended [Maginn2019]_ where the "middle" segment can be identified as +having a slope of 1. + +.. code-block:: python + + plt.loglog(lagtimes, msd) + plt.show() + +Now that we have identified what segment of our MSD to analyse, let's compute +a self-diffusivity. + +Computing Self-Diffusivity +-------------------------------- +Self-diffusivity is closely related to the MSD. + +.. math:: + + D_d = \frac{1}{2d} \lim_{t \to \infty} \frac{d}{dt} MSD(r_{d}) + +From the MSD, self-diffusivities :math:`D` with the desired dimensionality +:math:`d` can be computed by fitting the MSD with respect to the lag-time to +a linear model. An example of this is shown below, using the MSD computed in +the example above. The segment between :math:`\tau = 20` and :math:`\tau = 60` +is used to demonstrate selection of a MSD segment. + +.. code-block:: python + + from scipy.stats import linregress + start_time = 20 + start_index = int(start_time/timestep) + end_time = 60 + linear_model = linregress(lagtimes[start_index:end_index], + msd[start_index:end_index]) + slope = linear_model.slope + error = linear_model.rvalue + # dim_fac is 3 as we computed a 3D msd with 'xyz' + D = slope * 1/(2*MSD.dim_fac) + +We have now computed a self-diffusivity! + + +Notes +_____ + +There are several factors that must be taken into account when setting up and +processing trajectories for computation of self-diffusivities. +These include specific instructions around simulation settings, using +unwrapped trajectories and maintaining a relatively small elapsed time between +saved frames. Additionally, corrections for finite size effects are sometimes +employed along with various means of estimating errors [Yeh2004]_ [Bulow2020]_. +The reader is directed to the following review, which describes many of the +common pitfalls [Maginn2019]_. There are other ways to compute +self-diffusivity, such as from a Green-Kubo integral. At this point in time, +these methods are beyond the scope of this module. + + +Note also that computation of MSDs is highly memory intensive. If this is +proving a problem, judicious use of the ``start``, ``stop``, ``step`` keywords to control which frames are incorporated may be required. + +References +---------- + +.. [Maginn2019] Maginn, E. J., Messerly, R. A., Carlson, D. J.; Roe, D. R., + Elliott, J. R. Best Practices for Computing Transport + Properties 1. Self-Diffusivity and Viscosity from Equilibrium + Molecular Dynamics [Article v1.0]. Living J. Comput. Mol. Sci. + 2019, 1 (1). + +.. [Yeh2004] Yeh, I. C.; Hummer, G. System-Size Dependence of Diffusion + Coefficients and Viscosities from Molecular Dynamics + Simulations with Periodic Boundary Conditions. + J. Phys. Chem. B 2004, 108 (40), 15873–15879. + +.. [Bulow2020] von Bülow, S.; Bullerjahn, J. T.; Hummer, G. Systematic + Errors in Diffusion Coefficients from Long-Time Molecular + Dynamics Simulations at Constant Pressure. 2020. + arXiv:2003.09205 [Cond-Mat, Physics:Physics]. + + + +Classes +------- + +.. autoclass:: EinsteinMSD + :members: + :inherited-members: + +""" + +import numpy as np +import logging +from ..due import due, Doi +from .base import AnalysisBase + +logger = logging.getLogger('MDAnalysis.analysis.msd') + +due.cite(Doi("10.21105/joss.00877"), + description="Mean Squared Displacements with tidynamics", + path="MDAnalysis.analysis.msd", + cite_module=True) +due.cite(Doi("10.1051/sfn/201112010"), + description="FCA fast correlation algorithm", + path="MDAnalysis.analysis.msd", + cite_module=True) +del Doi + + +class EinsteinMSD(AnalysisBase): + r"""Class to calculate Mean Squared Displacement by the Einstein relation. + + Parameters + ---------- + u : Universe + An MDAnalysis :class:`Universe`. + select : str + A selection string. Defaults to `None` in which case + all atoms are selected. + msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'} + Desired dimensions to be included in the MSD. Defaults to 'xyz'. + fft : bool + If ``True``, uses a fast FFT based algorithm for computation of + the MSD. Otherwise, use the simple "windowed" algorithm. + The tidynamics package is required for `fft=True`. + Defaults to ``True``. + + Attributes + ---------- + dim_fac : int + Dimensionality :math:`d` of the MSD. + timeseries : :class:`numpy.ndarray` + The averaged MSD over all the particles with respect to lag-time. + msd_per_particle : :class:`numpy.ndarray` + The MSD of each individual particle with respect to lag-time. + n_frames : int + Number of frames included in the analysis. + n_particles : int + Number of particles MSD was calculated over. + + """ + + def __init__(self, u, select, msd_type='xyz', fft=True, **kwargs): + r""" + Parameters + ---------- + u : Universe + An MDAnalysis :class:`Universe`. + select : str + A selection string. Defaults to `None` in which case + all atoms are selected. + msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'} + Desired dimensions to be included in the MSD. Defaults to 'xyz'. + fft : bool + If ``True``, uses a fast FFT based algorithm for computation of + the MSD. Otherwise, use the simple "windowed" algorithm. + The tidynamics package is required for `fft=True`. + Defaults to ``True``. + + """ + self.u = u + + super(EinsteinMSD, self).__init__(self.u.trajectory, **kwargs) + + # args + self.select = select + self.msd_type = msd_type + self._parse_msd_type() + self.fft = fft + + # local + self._atoms = self.u.select_atoms(self.select) + self.n_particles = len(self._atoms) + self._position_array = None + + # result + self.msds_by_particle = None + self.timeseries = None + + def _prepare(self): + # self.n_frames only available here + # these need to be zeroed prior to each run() call + self.msds_by_particle = np.zeros((self.n_frames, self.n_particles)) + self._position_array = np.zeros( + (self.n_frames, self.n_particles, self.dim_fac)) + # self.timeseries not set here + + def _parse_msd_type(self): + r""" Sets up the desired dimensionality of the MSD. + + """ + keys = {'x': [0], 'y': [1], 'z': [2], 'xy': [0, 1], + 'xz': [0, 2], 'yz': [1, 2], 'xyz': [0, 1, 2]} + + self.msd_type = self.msd_type.lower() + + try: + self._dim = keys[self.msd_type] + except KeyError: + raise ValueError( + 'invalid msd_type: {} specified, please specify one of xyz, ' + 'xy, xz, yz, x, y, z'.format(self.msd_type)) + + self.dim_fac = len(self._dim) + + def _single_frame(self): + r""" Constructs array of positions for MSD calculation. + + """ + # shape of position array set here, use span in last dimension + # from this point on + self._position_array[self._frame_index] = ( + self._atoms.positions[:, self._dim]) + + def _conclude(self): + if self.fft: + self._conclude_fft() + else: + self._conclude_simple() + + def _conclude_simple(self): + r""" Calculates the MSD via the simple "windowed" algorithm. + + """ + lagtimes = np.arange(1, self.n_frames) + positions = self._position_array.astype(np.float64) + for lag in lagtimes: + disp = positions[:-lag, :, :] - positions[lag:, :, :] + sqdist = np.square(disp).sum(axis=-1) + self.msds_by_particle[lag, :] = np.mean(sqdist, axis=0) + self.timeseries = self.msds_by_particle.mean(axis=1) + + def _conclude_fft(self): # with FFT, np.float64 bit prescision required. + r""" Calculates the MSD via the FCA fast correlation algorithm. + + """ + try: + import tidynamics + except ImportError: + raise ImportError("""ERROR --- tidynamics was not found! + + tidynamics is required to compute an FFT based MSD (default) + + try installing it using pip eg: + + pip install tidynamics + + or set fft=False""") + + positions = self._position_array.astype(np.float64) + for n in range(self.n_particles): + self.msds_by_particle[:, n] = tidynamics.msd( + positions[:, n, :]) + self.timeseries = self.msds_by_particle.mean(axis=1) diff --git a/package/doc/sphinx/source/documentation_pages/analysis/msd.rst b/package/doc/sphinx/source/documentation_pages/analysis/msd.rst new file mode 100644 index 00000000000..4fcc341139b --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/msd.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.analysis.msd diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index 19081f1b5dd..ac4c7c94f28 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -101,14 +101,24 @@ Polymers Structure ========= +Macromolecules +-------------- + .. toctree:: :maxdepth: 1 analysis/gnm analysis/helanal - analysis/rdf analysis/dihedrals +Liquids +------- + +.. toctree:: + :maxdepth: 1 + + analysis/rdf + analysis/msd Volumetric analysis =================== diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index 594f0719748..1a00a2d8b99 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -166,7 +166,19 @@ If you use :meth:`~MDAnalysis.analysis.pca.PCA.cumulative_overlap` or .. _`10.1016/j.str.2007.12.011`: https://dx.doi.org/10.1016/j.str.2007.12.011 +If you use the Mean Squared Displacement analysis code in +:mod:`MDAnalysis.analysis.msd` please cite [Calandri2011]_ and [Buyl2018]_. +.. [Calandri2011] Calandrini, V., Pellegrini, E., Calligari, P., Hinsen, K., Kneller, G. R. + NMoldyn-Interfacing Spectroscopic Experiments, Molecular Dynamics Simulations and Models for Time Correlation Functions. + *Collect. SFN*, **12**, 201–232 (2011). doi: `10.1051/sfn/201112010`_ + +.. _`10.1051/sfn/201112010`: https://doi.org/10.1051/sfn/201112010 + +.. [Buyl2018] Buyl, P. tidynamics: A tiny package to compute the dynamics of stochastic and molecular simulations. Journal of Open Source Software, + 3(28), 877 (2018). doi: `10.21105/joss.00877`_ + +.. _`10.21105/joss.00877`: https://doi.org/10.21105/joss.00877 .. _citations-using-duecredit: diff --git a/package/doc/sphinx/source/images/msd_demo_plot.png b/package/doc/sphinx/source/images/msd_demo_plot.png new file mode 100644 index 00000000000..f3ecc25c252 Binary files /dev/null and b/package/doc/sphinx/source/images/msd_demo_plot.png differ diff --git a/package/setup.py b/package/setup.py index 371be571dc7..b97524d0e85 100755 --- a/package/setup.py +++ b/package/setup.py @@ -559,6 +559,7 @@ def dynamic_author_list(): 'scipy>=1.0.0', 'matplotlib>=1.5.1', 'mock', + 'tidynamics>=1.0.0', 'tqdm>=4.43.0', ] if not os.name == 'nt': diff --git a/testsuite/MDAnalysisTests/analysis/test_msd.py b/testsuite/MDAnalysisTests/analysis/test_msd.py new file mode 100644 index 00000000000..7b6090021b4 --- /dev/null +++ b/testsuite/MDAnalysisTests/analysis/test_msd.py @@ -0,0 +1,215 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + + +import MDAnalysis as mda +from MDAnalysis.analysis.msd import EinsteinMSD as MSD + +from numpy.testing import (assert_almost_equal, assert_equal) +import numpy as np + +from MDAnalysisTests.datafiles import PSF, DCD, RANDOM_WALK, RANDOM_WALK_TOPO + +import pytest +import tidynamics + + +@pytest.fixture(scope='module') +def SELECTION(): + selection = 'backbone and name CA and resid 1-10' + return selection + + +@pytest.fixture(scope='module') +def NSTEP(): + nstep = 5000 + return nstep + + +@pytest.fixture(scope='module') +def u(): + return mda.Universe(PSF, DCD) + + +@pytest.fixture(scope='module') +def random_walk_u(): + # 100x100 + return mda.Universe(RANDOM_WALK_TOPO, RANDOM_WALK) + + +@pytest.fixture(scope='module') +def msd(u, SELECTION): + # non fft msd + m = MSD(u, SELECTION, msd_type='xyz', fft=False) + m.run() + return m + + +@pytest.fixture(scope='module') +def msd_fft(u, SELECTION): + # fft msd + m = MSD(u, SELECTION, msd_type='xyz', fft=True) + m.run() + return m + + +@pytest.fixture(scope='module') +def step_traj(NSTEP): # constant velocity + x = np.arange(NSTEP) + traj = np.vstack([x, x, x]).T + traj_reshape = traj.reshape([NSTEP, 1, 3]) + u = mda.Universe.empty(1) + u.load_new(traj_reshape) + return u + + +@pytest.fixture(scope='module') +def step_traj_arr(NSTEP): # constant velocity + x = np.arange(NSTEP) + traj = np.vstack([x, x, x]).T + return traj + + +def random_walk_3d(NSTEP): + np.random.seed(1) + steps = -1 + 2*np.random.randint(0, 2, size=(NSTEP, 3)) + traj = np.cumsum(steps, axis=0) + traj_reshape = traj.reshape([NSTEP, 1, 3]) + u = mda.Universe.empty(1) + u.load_new(traj_reshape) + return u, traj + + +def characteristic_poly(n, d): # polynomial that describes unit step traj MSD + x = np.arange(0, n) + y = d*x*x + return y + + +def test_selection_works(msd): + # test some basic size and shape things + assert_equal(msd.n_particles, 10) + + +def test_fft_vs_simple_default(msd, msd_fft): + # testing on the PSF, DCD trajectory + timeseries_simple = msd.timeseries + timeseries_fft = msd_fft.timeseries + assert_almost_equal(timeseries_simple, timeseries_fft, decimal=4) + + +def test_fft_vs_simple_default_per_particle(msd, msd_fft): + # check fft and simple give same result per particle + per_particle_simple = msd.msds_by_particle + per_particle_fft = msd_fft.msds_by_particle + assert_almost_equal(per_particle_simple, per_particle_fft, decimal=4) + + +@pytest.mark.parametrize("dim", ['xyz', 'xy', 'xz', 'yz', 'x', 'y', 'z']) +def test_fft_vs_simple_all_dims(u, SELECTION, dim): + # check fft and simple give same result for each dimensionality + m_simple = MSD(u, SELECTION, msd_type=dim, fft=False) + m_simple.run() + timeseries_simple = m_simple.timeseries + m_fft = MSD(u, SELECTION, msd_type=dim, fft=True) + m_fft.run() + timeseries_fft = m_fft.timeseries + assert_almost_equal(timeseries_simple, timeseries_fft, decimal=4) + + +@pytest.mark.parametrize("dim", ['xyz', 'xy', 'xz', 'yz', 'x', 'y', 'z']) +def test_fft_vs_simple_all_dims_per_particle(u, SELECTION, dim): + # check fft and simple give same result for each particle in each dimension + m_simple = MSD(u, SELECTION, msd_type=dim, fft=False) + m_simple.run() + per_particle_simple = m_simple.msds_by_particle + m_fft = MSD(u, SELECTION, msd_type=dim, fft=True) + m_fft.run() + per_particle_fft = m_fft.msds_by_particle + assert_almost_equal(per_particle_simple, per_particle_fft, decimal=4) + + +@pytest.mark.parametrize("dim, dim_factor", \ +[('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1), ('z', 1)]) +def test_simple_step_traj_all_dims(step_traj, NSTEP, dim, dim_factor): + # testing the "simple" algorithm on constant velocity trajectory + # should fit the polynomial y=dim_factor*x**2 + m_simple = MSD(step_traj, 'all', msd_type=dim, fft=False) + m_simple.run() + poly = characteristic_poly(NSTEP, dim_factor) + assert_almost_equal(m_simple.timeseries, poly, decimal=4) + +@pytest.mark.parametrize("dim, dim_factor", \ +[('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1), ('z', 1)]) +def test_simple_start_stop_step_all_dims(step_traj, NSTEP, dim, dim_factor): + # testing the "simple" algorithm on constant velocity trajectory + # test start stop step is working correctly + m_simple = MSD(step_traj, 'all', msd_type=dim, fft=False) + m_simple.run(start=10, stop=1000, step=10) + poly = characteristic_poly(NSTEP, dim_factor) + # polynomial must take offset start into account + assert_almost_equal(m_simple.timeseries, poly[0:990:10], decimal=4) + + +@pytest.mark.parametrize("dim, dim_factor", \ +[('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1), ('z', 1)]) +def test_fft_step_traj_all_dims(step_traj, NSTEP, dim, dim_factor): + # testing the fft algorithm on constant velocity trajectory + # this should fit the polynomial y=dim_factor*x**2 + # fft based tests require a slight decrease in expected prescision + # primarily due to roundoff in fft(ifft()) calls. + # relative accuracy expected to be around ~1e-12 + m_simple = MSD(step_traj, 'all', msd_type=dim, fft=True) + m_simple.run() + poly = characteristic_poly(NSTEP, dim_factor) + # this was relaxed from decimal=4 for numpy=1.13 test + assert_almost_equal(m_simple.timeseries, poly, decimal=3) + +@pytest.mark.parametrize("dim, dim_factor", \ +[('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1), ('z', 1)]) +def test_fft_start_stop_step_all_dims(step_traj, NSTEP, dim, dim_factor): + # testing the fft algorithm on constant velocity trajectory + # test start stop step is working correctly + m_simple = MSD(step_traj, 'all', msd_type=dim, fft=True) + m_simple.run(start=10, stop=1000, step=10) + poly = characteristic_poly(NSTEP, dim_factor) + # polynomial must take offset start into account + assert_almost_equal(m_simple.timeseries, poly[0:990:10], decimal=3) + + +def test_random_walk_u_simple(random_walk_u): + # regress against random_walk test data + msd_rw = MSD(random_walk_u, 'all', msd_type='xyz', fft=False) + msd_rw.run() + norm = np.linalg.norm(msd_rw.timeseries) + val = 3932.39927487146 + assert_almost_equal(norm, val, decimal=5) + + +def test_random_walk_u_fft(random_walk_u): + # regress against random_walk test data + msd_rw = MSD(random_walk_u, 'all', msd_type='xyz', fft=True) + msd_rw.run() + norm = np.linalg.norm(msd_rw.timeseries) + val = 3932.39927487146 + assert_almost_equal(norm, val, decimal=5) diff --git a/testsuite/setup.py b/testsuite/setup.py index 67d3b3d25ca..1b4d88b181a 100755 --- a/testsuite/setup.py +++ b/testsuite/setup.py @@ -185,6 +185,7 @@ def run(self): 'hypothesis', 'psutil>=4.0.2', 'mock>=2.0.0', # replace with unittest.mock in python 3 only version + 'tidynamics>=1.0.0' ], # had 'KeyError' as zipped egg (2MB savings are not worth the # trouble)