Nongaussian parameter

package/CHANGELOG

@@ -31,6 +31,7 @@ Fixes
     (e.g. bonds, angles) (PR #3779).
 
 Enhancements
+  * Nongaussian parameter available from msd.run(fft=False), (#3895)
   * LAMMPSDump Reader optionally unwraps trajectories with image flags upon 
     loading (#3843
   * LAMMPSDump Reader now imports velocities and forces (#3843)

package/MDAnalysis/analysis/msd.py

@@ -234,6 +234,7 @@
     Maginn2019
     Yeh2004
     Bulow2020
+    Fuchs1998
 
 
 Classes
@@ -267,6 +268,10 @@
 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
+    nonbonded parameter is also computed to capture heterogiety as different
+    time scales. See :cite:p:`Fuchs1998` for more information. 
+
     Parameters
     ----------
     u : Universe or AtomGroup
@@ -291,6 +296,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 : :class:`numpy.ndarray`
+        Only available when `fft==False`. The averaged nongaussian 
+        parameter  over all the particles with respect to lag-time.
+    results.msds_by_nongaussian : :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
@@ -299,6 +310,9 @@ class EinsteinMSD(AnalysisBase):
         Number of particles MSD was calculated over.
 
 
+    .. versionchanged:: 2.4.0
+       Added second order non-gaussian parameter to output for "simple"
+       algorithm.
     .. versionadded:: 2.0.0
     """
 
@@ -338,12 +352,18 @@ def __init__(self, u, select='all', msd_type='xyz', fft=True, **kwargs):
         # result
         self.results.msds_by_particle = None
         self.results.timeseries = None
+        if not fft:
+            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))
+        if not self.fft:
+            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
@@ -382,16 +402,28 @@ 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)
         positions = self._position_array.astype(np.float64)
         for lag in lagtimes:
             disp = positions[:-lag, :, :] - positions[lag:, :, :]
             sqdist = np.square(disp).sum(axis=-1)
-            self.results.msds_by_particle[lag, :] = np.mean(sqdist, axis=0)
+            tmp = np.mean(sqdist, axis=0)
+            self.results.msds_by_particle[lag, :] = tmp
+            self.results.nongaussian_by_particle[lag, :] = 3.0/5.0*np.mean(
+                np.square(sqdist), axis=0)/np.square(tmp) - 1.0
+
         self.results.timeseries = self.results.msds_by_particle.mean(axis=1)
+        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}
+}
+
 @article{MichaudAgrawal2011,
     author = {Michaud-Agrawal, Naveen and Denning, Elizabeth J. and Woolf, Thomas B. and Beckstein, Oliver},
     title = {MDAnalysis: A toolkit for the analysis of molecular dynamics simulations},
@@ -734,4 +746,4 @@ @incollection{Waveren2005
     url = {http://www.intel.com/cd/ids/developer/asmo-na/eng/293748.htm},
     booktitle = {},
     id = {transformations}
-}
+}

testsuite/MDAnalysisTests/analysis/test_msd.py

@@ -123,7 +123,12 @@ def test_simple_step_traj_all_dims(self, step_traj, NSTEP, dim,
         m_simple = MSD(step_traj, 'all', msd_type=dim, fft=False)
         m_simple.run()
         poly = characteristic_poly(NSTEP, dim_factor)
+        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),
@@ -147,6 +152,15 @@ def test_random_walk_u_simple(self, random_walk_u):
         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 
+        # postiive values, but converge to zero at long times.
+        assert_almost_equal(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]), 
+                            decimal=5)
 
 
 @pytest.mark.skipif(import_not_available("tidynamics"),