Nongaussian parameter

package/MDAnalysis/analysis/msd.py

@@ -20,7 +20,6 @@
 # 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`
 ==============================================================
@@ -227,7 +226,18 @@
 References
 ----------
 
+<<<<<<< HEAD
+.. bibliography::
+    :filter: False
+    :style: MDA
+
+    Maginn2019
+    Yeh2004
+    Bulow2020
+    Fuchs1998
+=======
 .. footbibliography::
+>>>>>>> develop
 
 
 Classes
@@ -262,21 +272,34 @@
 class EinsteinMSD(AnalysisBase):
     r"""Class to calculate Mean Squared Displacement by the Einstein relation.
 
+    If `fft=False` so that the "windowed" algorithm is used, the second order
+    nongaussian parameter can also computed with ``nongaussian==True`` to 
+    capture heterogeneity as different time scales. See :cite:p:`Fuchs1998` 
+    for more information. 
+
+    .. math::
+
+       \alpha_{2}(r_{d}) = \frac{3}{5} \frac{ \bigg{\langle} \frac{1}{N} \sum_{i=1}^{N} |r_{d}
+        - r_{d}(t_0)|^4 \bigg{\rangle} } { \bigg{\langle} \frac{1}{N} \sum_{i=1}^{N} |r_{d}
+        - r_{d}(t_0)|^2 \bigg{\rangle}^2 } - 1
+
     Parameters
     ----------
     u : Universe or AtomGroup
         An MDAnalysis :class:`Universe` or :class:`AtomGroup`.
         Note that :class:`UpdatingAtomGroup` instances are not accepted.
-    select : str
+    select : str (optional)
         A selection string. Defaults to "all" in which case
         all atoms are selected.
-    msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'}
+    msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'} (optional)
         Desired dimensions to be included in the MSD. Defaults to 'xyz'.
-    fft : bool
+    fft : bool (optional)
         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``.
+    nongaussian : bool (optional)
+        If ``True`` the second order nongaussian parameter is calculated.
 
     Attributes
     ----------
@@ -286,6 +309,12 @@ class EinsteinMSD(AnalysisBase):
         The averaged MSD over all the particles with respect to lag-time.
     results.msds_by_particle : :class:`numpy.ndarray`
         The MSD of each individual particle with respect to lag-time.
+    results.nongaussian_parameter : :class:`numpy.ndarray`
+        Can only be ``True`` when `fft==False`. The averaged nongaussian 
+        parameter  over all the particles with respect to lag-time.
+    results.nongaussian_by_particle : :class:`numpy.ndarray`
+        Only available when `fft==False`. The nongaussian parameter 
+        of each individual particle with respect to lag-time.
     ag : :class:`AtomGroup`
         The :class:`AtomGroup` resulting from your selection
     n_frames : int
@@ -294,24 +323,36 @@ class EinsteinMSD(AnalysisBase):
         Number of particles MSD was calculated over.
 
 
+    .. versionchanged:: 2.4.0
+       Added optional second order non-gaussian parameter to output for 
+       "simple" algorithm.
     .. versionadded:: 2.0.0
     """
 
-    def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
+    def __init__(self,
+                 u,
+                 select='all',
+                 msd_type='xyz',
+                 fft=True,
+                 nongaussian=False,
+                 **kwargs):
         r"""
         Parameters
         ----------
         u : Universe or AtomGroup
             An MDAnalysis :class:`Universe` or :class:`AtomGroup`.
-        select : str
+        select : str (optional)
             A selection string. Defaults to "all" in which case
             all atoms are selected.
-        msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'}
-            Desired dimensions to be included in the MSD.
-        fft : bool
+        msd_type : {'xyz', 'xy', 'yz', 'xz', 'x', 'y', 'z'} (optional)
+            Desired dimensions to be included in the MSD. Defaults to 'xyz'.
+        fft : bool (optional)
             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`.
+        nongaussian : bool (optional)
+            If ``True`` the nongaussian parameter is calculated.
+            Can only be ``True`` when `fft==False`. The averaged nongaussian 
         """
         if isinstance(u, groups.UpdatingAtomGroup):
             raise TypeError("UpdatingAtomGroups are not valid for MSD "
@@ -324,6 +365,7 @@ def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
         self.msd_type = msd_type
         self._parse_msd_type()
         self.fft = fft
+        self.nongaussian = nongaussian
 
         # local
         self.ag = u.select_atoms(self.select)
@@ -333,12 +375,21 @@ def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
         # result
         self.results.msds_by_particle = None
         self.results.timeseries = None
+        if fft and nongaussian:
+            raise ValueError("The nongaussian parameter can only be computed"
+                             " when `fft=False`")
+        elif nongaussian:
+            self.results.nongaussian_by_particle = None
+            self.results.nongaussian_parameter = None
 
     def _prepare(self):
         # self.n_frames only available here
         # these need to be zeroed prior to each run() call
-        self.results.msds_by_particle = np.zeros((self.n_frames,
-                                                  self.n_particles))
+        self.results.msds_by_particle = np.zeros(
+            (self.n_frames, self.n_particles))
+        if self.nongaussian:
+            self.results.nongaussian_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.results.timeseries not set here
@@ -347,8 +398,15 @@ 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]}
+        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()
 
@@ -377,7 +435,8 @@ def _conclude(self):
             self._conclude_simple()
 
     def _conclude_simple(self):
-        r""" Calculates the MSD via the simple "windowed" algorithm.
+        r""" Calculates the MSD via the simple "windowed" algorithm 
+        with nongaussian parameter.
 
         """
         lagtimes = np.arange(1, self.n_frames)
@@ -386,7 +445,18 @@ def _conclude_simple(self):
             disp = positions[:-lag, :, :] - positions[lag:, :, :]
             sqdist = np.square(disp).sum(axis=-1)
             self.results.msds_by_particle[lag, :] = np.mean(sqdist, axis=0)
+            if self.nongaussian:
+                self.results.nongaussian_by_particle[
+                    lag, :] = 3.0 / 5.0 * np.mean(
+                        np.square(sqdist), axis=0) / np.square(
+                            self.results.msds_by_particle[lag, :]) - 1.0
+
         self.results.timeseries = self.results.msds_by_particle.mean(axis=1)
+        if self.nongaussian:
+            self.results.nongaussian_parameter = (3.0 / 5.0 * np.mean(
+                (self.results.nongaussian_by_particle + 1.0) * 5.0 / 3.0 *
+                np.square(self.results.msds_by_particle),
+                axis=1) / np.square(self.results.timeseries) - 1.0)
 
     def _conclude_fft(self):  # with FFT, np.float64 bit prescision required.
         r""" Calculates the MSD via the FCA fast correlation algorithm.

package/doc/sphinx/source/references.bib

@@ -1,3 +1,15 @@
+@article{Fuchs1998,
+        author = {Fuchs, M. and Gotze, W. and Mayr, M. R.},
+        doi = {10.1103/PhysRevE.58.3384},
+        journal = {Physical Review E},
+        Number = {3},
+        Pages = {3384--3399},
+        Publisher = {American Physical Society},
+        Title = {Asymptotic laws for tagged-particle motion in glassy systems},
+        Volume = {58},
+        year = {1998}
+}
+
 @book{Gray1984,
 	address = {Oxford ; New York},
 	series = {The {International} series of monographs on chemistry},

testsuite/MDAnalysisTests/analysis/test_msd.py

@@ -21,11 +21,10 @@
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 #
 
-
 from MDAnalysis.analysis.msd import EinsteinMSD as MSD
 import MDAnalysis as mda
 
-from numpy.testing import (assert_almost_equal, assert_equal)
+from numpy.testing import (assert_almost_equal, assert_equal, assert_allclose)
 import numpy as np
 
 from MDAnalysisTests.datafiles import PSF, DCD, RANDOM_WALK, RANDOM_WALK_TOPO
@@ -86,7 +85,7 @@ def test_notidynamics(u, SELECTION):
 def characteristic_poly(n, d):
     # polynomial that describes unit step traj MSD
     x = np.arange(0, n)
-    y = d*x*x
+    y = d * x * x
     return y
 
 
@@ -112,23 +111,32 @@ def test_msdtype_error(self, u, SELECTION, msdtype):
         with pytest.raises(ValueError, match=errmsg):
             m = MSD(u, SELECTION, msd_type=msdtype)
 
-    @pytest.mark.parametrize("dim, dim_factor", [
-        ('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1),
-        ('z', 1)
-    ])
+    @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(self, 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 = MSD(step_traj,
+                       'all',
+                       msd_type=dim,
+                       fft=False,
+                       nongaussian=True)
         m_simple.run()
         poly = characteristic_poly(NSTEP, dim_factor)
+        # In the simple alorithm with an exact solution, mean(r(d)^4) == MSD^2
+        # So, alpha = 3/5 - 1 = -0.4
+        alpha = -0.4 * np.ones(NSTEP)
+        alpha[0] = np.nan
         assert_almost_equal(m_simple.results.timeseries, poly, decimal=4)
+        assert_almost_equal(m_simple.results.nongaussian_parameter,
+                            alpha,
+                            decimal=5)
 
-    @pytest.mark.parametrize("dim, dim_factor", [
-        ('xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1),
-        ('z', 1)
-    ])
+    @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(self, step_traj, NSTEP, dim,
                                              dim_factor):
         # testing the "simple" algorithm on constant velocity trajectory
@@ -137,16 +145,38 @@ def test_simple_start_stop_step_all_dims(self, step_traj, NSTEP, dim,
         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.results.timeseries, poly[0:990:10],
+        assert_almost_equal(m_simple.results.timeseries,
+                            poly[0:990:10],
                             decimal=4)
 
     def test_random_walk_u_simple(self, random_walk_u):
         # regress against random_walk test data
-        msd_rw = MSD(random_walk_u, 'all', msd_type='xyz', fft=False)
+        msd_rw = MSD(random_walk_u,
+                     'all',
+                     msd_type='xyz',
+                     fft=False,
+                     nongaussian=True)
         msd_rw.run()
         norm = np.linalg.norm(msd_rw.results.timeseries)
         val = 3932.39927487146
         assert_almost_equal(norm, val, decimal=5)
+        # nongaussian parameter may have negative values (above -2.5) and
+        # positive values, but converge to zero at long times.
+        assert_allclose(msd_rw.results.nongaussian_parameter[:20],
+                        np.array([
+                            np.nan, 0.00049, 0.00048, 0.00536, -0.00299,
+                            0.00014, -0.0046, -0.00612, -0.00307, 0.00033,
+                            0.00584, 0.0131, 0.01516, 0.01666, 0.01924,
+                            0.02052, 0.01961, 0.01743, 0.01355, 0.00997
+                        ]),
+                        atol=1e-5)
+
+    def test_no_fft_if_nongaussian(self, random_walk_u):
+        # regress against random_walk test data
+        errmsg = ("The nongaussian parameter can only be computed"
+                  " when `fft=False`")
+        with pytest.raises(ValueError, match=errmsg):
+            msd_rw = MSD(random_walk_u, 'all', fft=True, nongaussian=True)
 
 
 @pytest.mark.skipif(import_not_available("tidynamics"),
@@ -195,10 +225,9 @@ def test_fft_vs_simple_all_dims_per_particle(self, u, SELECTION, dim):
         per_particle_fft = m_fft.results.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)
-    ])
+    @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(self, step_traj, NSTEP, dim, dim_factor):
         # testing the fft algorithm on constant velocity trajectory
         # this should fit the polynomial y=dim_factor*x**2
@@ -211,10 +240,9 @@ def test_fft_step_traj_all_dims(self, step_traj, NSTEP, dim, dim_factor):
         # this was relaxed from decimal=4 for numpy=1.13 test
         assert_almost_equal(m_simple.results.timeseries, poly, decimal=3)
 
-    @pytest.mark.parametrize("dim, dim_factor", [(
-        'xyz', 3), ('xy', 2), ('xz', 2), ('yz', 2), ('x', 1), ('y', 1),
-        ('z', 1)
-    ])
+    @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(self, step_traj, NSTEP, dim,
                                           dim_factor):
         # testing the fft algorithm on constant velocity trajectory
@@ -223,7 +251,8 @@ def test_fft_start_stop_step_all_dims(self, step_traj, NSTEP, dim,
         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.results.timeseries, poly[0:990:10],
+        assert_almost_equal(m_simple.results.timeseries,
+                            poly[0:990:10],
                             decimal=3)
 
     def test_random_walk_u_fft(self, random_walk_u):