Source code for MDAnalysis.lib.correlations

# -*- 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
#

"""Correlations utilities --- :mod:`MDAnalysis.lib.correlations`
=================================================================================


:Authors: Paul Smith & Mateusz Bieniek
:Year: 2020
:Copyright: GNU Public License v2

.. versionadded:: 1.0.0

This module is primarily for internal use by other analysis modules. It
provides functionality for calculating the time autocorrelation function
of a binary variable (i.e one that is either true or false at each
frame for a given atom/molecule/set of molecules). This module includes
functions for calculating both the time continuous autocorrelation and
the intermittent autocorrelation. The function :func:`autocorrelation`
calculates the continuous autocorrelation only. The data may be
pre-processed using the function :func:`intermittency` in order to
acount for intermittency before passing the results to
:func:`autocorrelation`.

This module is inspired by seemingly disparate analyses that rely on the same
underlying calculation, including the survival probability of water around
proteins [Araya-Secchi2014]_, hydrogen bond lifetimes [Gowers2015]_, [Araya-Secchi2014]_,
and the rate of cholesterol flip-flop in lipid bilayers [Gu2019]_.

.. seeAlso::

    Analysis tools that make use of modules:

        * :class:`MDAnalysis.analysis.waterdynamics.SurvivalProbability`
            Calculates the continuous or intermittent survival probability
            of an atom group in a region of interest.

        * :class:`MDAnalysis.analysis.hbonds.hbond_analysis`
            Calculates the continuous or intermittent hydrogen bond
            lifetime.

.. rubric:: References

.. [Gu2019] Gu, R.-X.; Baoukina, S.; Tieleman, D. P. (2019)
            Cholesterol Flip-Flop in Heterogeneous Membranes.
            J. Chem. Theory Comput. 15 (3), 2064–2070.
            https://doi.org/10.1021/acs.jctc.8b00933.

"""

from __future__ import absolute_import, division
import numpy as np
from copy import deepcopy


[docs]def autocorrelation(list_of_sets, tau_max, window_step=1): r"""Implementation of a discrete autocorrelation function. The autocorrelation of a property :math:`x` from a time :math:`t=t_0` to :math:`t=t_0 + \tau` is given by: .. math:: C(\tau) = \langle \frac{ x(t_0)x(t_0 +\tau) }{ x(t_0)x(t_0) } \rangle where :math:`x` may represent any property of a particle, such as velocity or potential energy. This function is an implementation of a special case of the time autocorrelation function in which the property under consideration can be encoded with indicator variables, :math:`0` and :math:`1`, to represent the binary state of said property. This special case is often referred to as the survival probability (:math:`S(\tau)`). As an example, in calculating the survival probability of water molecules within 5 Å of a protein, each water molecule will either be within this cutoff range (:math:`1`) or not (:math:`0`). The total number of water molecules within the cutoff at time :math:`t_0` will be given by :math:`N(t_0)`. Other cases include the Hydrogen Bond Lifetime as well as the translocation rate of cholesterol across a bilayer. The survival probability of a property of a set of particles is given by: .. math:: S(\tau) = \langle \frac{ N(t_0, t_0 + \tau )} { N(t_0) }\rangle where :math:`N(t0)` is the number of particles at time :math:`t_0` for which the feature is observed, and :math:`N(t0, t_0 + \tau)` is the number of particles for which this feature is present at every frame from :math:`t_0` to :math:`N(t0, t_0 + \tau)`. The angular brackets represent an average over all time origins, :math:`t_0`. See [Araya-Secchi2014]_ for a description survival probability. Parameters ---------- list_of_sets : list List of sets. Each set corresponds to data from a single frame. Each element in a set may be, for example, an atom id or a tuple of atoms ids. In the case of calculating the survival probability of water around a protein, these atom ids in a given set will be those of the atoms which are within a cutoff distance of the protein at a given frame. tau_max : int The last tau (lag time, inclusive) for which to calculate the autocorrelation. e.g if tau_max = 20, the survival probability will be calculated over 20 frames. window_step : int, optional The step size for t0 to perform autocorrelation. Ideally, window_step will be larger than tau_max to ensure independence of each window for which the calculation is performed. Default is 1. Returns -------- tau_timeseries : list of int the values of tau for which the autocorrelation was calculated timeseries : list of int the autocorelation values for each of the tau values timeseries_data : list of list of int the raw data from which the autocorrelation is computed, i.e :math:`S(\tau)` at each window. This allows the time dependant evolution of :math:`S(\tau)` to be investigated. .. versionadded:: 0.19.2 """ # check types if (type(list_of_sets) != list and len(list_of_sets) != 0) or type(list_of_sets[0]) != set: raise TypeError("list_of_sets must be a one-dimensional list of sets") # pragma: no cover # Check dimensions of parameters if len(list_of_sets) < tau_max: raise ValueError("tau_max cannot be greater than the length of list_of_sets") # pragma: no cover tau_timeseries = list(range(1, tau_max + 1)) timeseries_data = [[] for _ in range(tau_max)] # calculate autocorrelation for t in range(0, len(list_of_sets), window_step): Nt = len(list_of_sets[t]) if Nt == 0: continue for tau in tau_timeseries: if t + tau >= len(list_of_sets): break # continuous: IDs that survive from t to t + tau and at every frame in between Ntau = len(set.intersection(*list_of_sets[t:t + tau + 1])) timeseries_data[tau - 1].append(Ntau / float(Nt)) timeseries = [np.mean(x) for x in timeseries_data] # at time 0 the value has to be one tau_timeseries.insert(0, 0) timeseries.insert(0, 1) return tau_timeseries, timeseries, timeseries_data
[docs]def correct_intermittency(list_of_sets, intermittency): r"""Preprocess data to allow intermittent behaviour prior to calling :func:`autocorrelation`. Survival probabilty may be calculated either with a strict continuous requirement or a less strict intermittency. If calculating the survival probability water around a protein for example, in the former case the water must be within a cutoff distance of the protein at every frame from :math:`t_0` to :math:`t_0 + \tau` in order for it to be considered present at :math:`t_0 + \tau`. In the intermittent case, the water molecule is allowed to leave the region of interest for up to a specified consecutive number of frames whilst still being considered present at :math:`t_0 + \tau`. This function pre-processes data, such as the atom ids of water molecules within a cutoff distance of a protein at each frame, in order to allow for intermittent behaviour, with a single pass over the data. For example, if an atom is absent for a number of frames equal or smaller than the parameter :attr:`intermittency`, then this absence will be removed and thus the atom is considered to have not left. e.g 7,A,A,7 with `intermittency=2` will be replaced by 7,7,7,7, where A=absence. The returned data can be used as input to the function :func:`autocorrelation` in order to calculate the survival probability with a given intermittency. See [Gowers2015]_ for a description of intermittency in the calculation of hydrogen bond lifetimes. # TODO - is intermittency consitent with list of sets of sets? (hydrogen bonds) Parameters ---------- list_of_sets: list In the simple case of e.g survival probability, a list of sets of atom ids present at each frame, where a single set contains atom ids at a given frame, e.g [{0, 1}, {0}, {0}, {0, 1}] intermittency : int The maximum gap allowed. The default `intermittency=0` means that if the datapoint is missing at any frame, no changes are made to the data. With the value of `intermittency=2`, all datapoints missing for up to two consecutive frames will be instead be considered present. Returns ------- list_of_sets: list returns a new list with the IDs with added IDs which disappeared for <= :attr:`intermittency`. e.g If [{0, 1}, {0}, {0}, {0, 1}] is a list of sets of atom ids present at each frame and `intermittency=2`, both atoms will be considered present throughout and thus the returned list of sets will be [{0, 1}, {0, 1}, {0, 1}, {0, 1}]. """ if intermittency == 0: return list_of_sets list_of_sets = deepcopy(list_of_sets) for i, ids in enumerate(list_of_sets): # initially update each frame as seen 0 ago (now) seen_frames_ago = {i: 0 for i in ids} for j in range(1, intermittency + 2): for atomid in seen_frames_ago.keys(): # no more frames if i + j >= len(list_of_sets): continue # if the atom is absent now if not atomid in list_of_sets[i + j]: # increase its absence counter seen_frames_ago[atomid] += 1 continue # the atom is found if seen_frames_ago[atomid] == 0: # the atom was present in the last frame continue # it was absent more times than allowed if seen_frames_ago[atomid] > intermittency: continue # the atom was absent but returned (within <= intermittency_value) # add it to the frames where it was absent. # Introduce the corrections. for k in range(seen_frames_ago[atomid], 0, -1): list_of_sets[i + j - k].add(atomid) seen_frames_ago[atomid] = 0 return list_of_sets