msd

mdcraft.lib.correlation.msd(r_i: np.ndarray[float_t], r_j: np.ndarray[float_t] | None = None, /, axis: int | None = None, *, average: bool = True, fft: bool = True) → np.ndarray[float_t]

Evaluates the mean squared displacement (MSD) or the cross mean squared displacement (CMSD) of positions \(\mathbf{r}_i(t)\) and \(\mathbf{r}_j(t)\).

Using the Einstein relation, the MSD for a set of positions \(\mathbf{r}_i(t)\) can be computed using

\[\mathrm{MSD}_i(\tau)=\langle[\mathbf{r}_i(t_0+\tau) -\mathbf{r}_i(t_0)]^2\rangle =\dfrac{1}{N_\tau}\sum_{k=1}^{N_\tau} [\mathbf{r}_i(t_k+\tau)-\mathbf{r}_i(t_k)]^2\]

where \(\tau\) is the time lag, \(t_j\) is an arbitrary reference time, and \(N_\tau\) is the number of possible reference times.

Similarly, the CMSD for two sets of positions \(\mathbf{r}_i(t)\) and \(\mathbf{r}_j(t)\) can be computed using

\[\begin{split}\mathrm{CMSD}_{ij}(\tau)&=\langle [\mathbf{r}_i(t_0+\tau)-\mathbf{r}_i(t_0)]\cdot [\mathbf{r}_j(t_0+\tau)-\mathbf{r}_j(t_0)]\rangle\\ &=\dfrac{1}{N_\tau}\sum_{k=1}^{N_\tau} [\mathbf{r}_i(t_k+\tau)-\mathbf{r}_i(t_k)]\cdot [\mathbf{r}_j(t_k+\tau)-\mathbf{r}_j(t_k)]\end{split}\]

To minimize statistical noise, the MSD/CMSD is calculated for and averaged over all possible reference times \(t_k\).

Alternatively, the MSD/CMSD can be efficiently computed using the fast convolution algorithm (FCA) [1] [2], which leverages the Wiener–Khinchin theorem. FCA uses fast Fourier transforms (FFT) and has a time complexity of \(\mathcal{O}(N\log{N})\).

Using FCA, the MSD for a set of positions \(\mathbf{r}_i(t)\) is computed using

\[\begin{split}\mathrm{MSD}_{i,m}&=\frac{1}{N_t-m} \sum_{k=0}^{N_t-m-1} [\mathbf{r}_{i,k+m}-\mathbf{r}_{i,k}]^2\\ &=\frac{1}{N_t-m}\sum_{k=0}^{N_t-m-1} \left[\mathbf{r}_{i,k+m}^2+\mathbf{r}_{i,k}^2\right] -\frac{2}{N_t-m}\sum_{k=0}^{N_t-m-1} \mathbf{r}_{i,k}\cdot\mathbf{r}_{i,k+m}\\ &=\mathrm{S}_{ii,m}-2\mathrm{R}_{ii,m}\end{split}\]

where \(i\) is the species index, \(m\) is the index corresponding to time lag \(\tau\), \(\mathrm{R}_{ii,m}\) is the autocorrelation of \(\mathbf{r}_i(t)\), and \(S_m\) is evaluated using the recursive relation

\[\begin{split}\begin{gather*} D_{ii,m}=\mathbf{r}_{i,m}^2\\ Q_{ii,-1}=2\sum_{k=0}^{N_t-1}D_{ii,k}\\ Q_{ii,m}=Q_{ii,m-1}-D_{ii,m-1}-D_{ii,N_t-m}\\ S_{ii,m}=\frac{Q_{ii,m}}{N_t-m} \end{gather*}\end{split}\]

Similarly, the CMSD for two sets of positions \(\mathbf{r}_i(t)\) and \(\mathbf{r}_j(t)\) is computed using

\[\mathrm{CMSD}_{ij,m}=S_{ij,m}-2\mathrm{R}_{ij,m}\]

where \(\mathrm{R}_{ij,m}\) is the cross-correlation of \(\mathbf{r}_i(t)\) and \(\mathbf{r}_j(t)\), and \(S_{ij,m}\) is evaluated using the recursive relation

\[\begin{split}\begin{gather*} D_{ij,m}=\mathbf{r}_{i,m}\cdot\mathbf{r}_{j,m}\\ Q_{ij,-1}=2\sum_{k=0}^{N_t-1}D_{ij,k}\\ Q_{ij,m}=Q_{ij,m-1}-D_{ij,m-1}-D_{ij,N_t-m}\\ S_{ij,m}=\frac{Q_{ij,m}}{N_t-m} \end{gather*}\end{split}\]

Note

r_i and r_j should be summed over all entities before being passed to this function if it is being used to evaluate the cross terms in the Onsager transport coefficients [3]

\[L_{ij}=\frac{1}{6k_\mathrm{B}T}\lim_{t\rightarrow\infty} \frac{d}{d\tau}\left\langle\sum_{\alpha=1}^{N_i} [\mathbf{r}_{i,\alpha}(t_0+\tau) -\mathbf{r}_{i,\alpha}(t_0)]\cdot \sum_{\beta=1}^{N_j}[\mathbf{r}_{j,\beta}(t_0+\tau) -\mathbf{r}_{j,\beta}(t_0)]\right\rangle\]

Parameters:

r_inumpy.ndarray, positional-only

Time evolution of individual or summed \(d\)-dimensional positions \(\mathbf{r}_i(t)\) for \(N\) entities over \(N_\mathrm{b}\) blocks of \(N_t\) times each.

Shape: \((N_t,d)\), \((N_t,N,d)\), \((N_\mathrm{b},N_t,d)\), or \((N_\mathrm{b},N_t,N,d)\).

Reference unit: \(\mathrm{nm}\).

r_jnumpy.ndarray, positional-only, optional

Time evolution of individual or summed \(d\)-dimensional positions \(\mathbf{r}_j(t)\) for another \(N\) entities over \(N_\mathrm{b}\) blocks of \(N_t\) times each.

Shape: Same as r_i.

Reference unit: \(\mathrm{nm}\).

axisint, optional

Axis along which time evolves. If not specified, the axis is determined automatically using the shape of r_i.

averagebool, keyword-only, default: True

Specifies whether to average the MSD/CMSD over all entities. Only available if r_i and r_j contain information for multiple entities.

fftbool, keyword-only, default: True

Specifies whether to use fast Fourier transforms (FFT) to evaluate the MSD/CMSD.

Returns:

dispnumpy.ndarray

\(\mathrm{MSD}_i(\tau)\) or \(\mathrm{CMSD}_{ij}(\tau)\).

Shape: Shape of r_i, except the last axis is no longer present. If average=True, the axis indexing the \(N\) entities is also no longer present.

Reference unit: \(\mathrm{nm}^2\).

References

[1]

Kneller, G. R.; Keiner, V.; Kneller, M.; Schiller, M. NMOLDYN: A Program Package for a Neutron Scattering Oriented Analysis of Molecular Dynamics Simulations. Computer Physics Communications 1995, 91 (1–3), 191–214. https://doi.org/10.1016/0010-4655(95)00048-K.

[2]

Calandrini, V.; Pellegrini, E.; Calligari, P.; Hinsen, K.; Kneller, G. R. NMoldyn - Interfacing Spectroscopic Experiments, Molecular Dynamics Simulations and Models for Time Correlation Functions. JDN 2011, 12, 201–232. https://doi.org/10.1051/sfn/201112010.

[3]

Fong, K. D.; Self, J.; McCloskey, B. D.; Persson, K. A. Onsager Transport Coefficients and Transference Numbers in Polyelectrolyte Solutions and Polymerized Ionic Liquids. Macromolecules 2020, 53 (21), 9503–9512. https://doi.org/10.1021/acs.macromol.0c02001.