next up previous contents
Next: Parameters Up: Mean-Square Displacement Previous: Mean-Square Displacement   Contents


Theory and implementation

Molecules in liquids and gases do not stay in the same place, but move constantly. This process is called diffusion and it happens quite naturally in fluids at equilibrium. During this process, the motion of an individual molecule does not follow a simple path [48]. As it travels, the molecule undergoes some collisions with other molecules which prevent it from following a straight line. If the path is examined in close detail, it will be seen to be a good approximation to a random walk. Mathematically, a random walk is a series of steps where each step is taken in a completely random direction from the one before. This kind of path was famously analysed by Albert Einstein in a study of Brownian motion. He showed that the MSD of a particle following a random walk is proportional to the time elapsed. This relationship can be written as
\begin{displaymath}
<r^2> = 6Dt + C
\end{displaymath} (4.14)

where $<r^2>$ is the MSD and t is the time. D and C are constants. The constant D defines the so-called diffusion coefficient. The figure 4.34 shows an example of a MSD analysis performed on a waterbox of 768 water molecules.
Figure 4.34: MSD calculated for a 100 ps MD simulation of 256 water molecules using NPT condition at 1 bar and 300 K.
\includegraphics[width=10cm]{Figures/msd_water.eps}

To get the diffusion coefficient out of this plot, the slope of the linear part of the plot should be calculated.

Defining,

\begin{displaymath}
\textbf{d}_{\alpha}(t,t_0) \doteq \textbf{R}_\alpha(t_0 + t) - \textbf{R}_\alpha(t_0).
\end{displaymath} (4.15)

the MSD of particle $\alpha$ can be defined as:
\begin{displaymath}
\Delta^2_\alpha(t) = \left\langle \textbf{d}^2_\alpha(t,t_0) \right\rangle_{t_0}
\end{displaymath} (4.16)

where $\textbf{R}_\alpha(t_0)$ and $\textbf{R}_\alpha(t_0 + t)$ are respectively the position of particle $\alpha$ at times $t_0$ and $t_0+t$. One can introduce a MSD with respect to a given axis n:
\begin{displaymath}
\Delta^2_\alpha(t,t_0;\textbf{n}) \doteq \left\langle
d^2_\alpha(t,\tau;\textbf{n})\right\rangle_{t_0}
\end{displaymath} (4.17)

with
\begin{displaymath}
d_{\alpha}(t,t_0;\textbf{n}) \doteq \textbf{n}\cdot \textbf{d}_{\alpha}(t,t_0).
\end{displaymath} (4.18)

The calculation of MSD is the standard way to obtain diffusion coefficients from MD simulations. Assuming Einstein-diffusion in the long time limit one has for isotropic systems
\begin{displaymath}
D_\alpha = \lim_{t\to\infty} \frac{1}{6t} \Delta^2_\alpha(t).
\end{displaymath} (4.19)

There exists also a well-known relation between the MSD and the velocity autocorrelation function. Writing $\textbf{d}_\alpha(t) = \int_{0}^{t}d\tau \textbf{v}_\alpha(\tau)$ in Eq. (4.16) one can show (see e.g. [50]) that

\begin{displaymath}
\Delta^2_\alpha(t) = 6 \int_{0}^{t}d\tau (t - \tau)
C_{vv ; \alpha\alpha}(\tau).
\end{displaymath} (4.20)

Using now the definition (4.19) of the diffusion coefficient one obtains the relation
\begin{displaymath}
D_\alpha = \int_{0}^{t}d\tau  C_{vv ; \alpha\alpha}(\tau).
\end{displaymath} (4.21)

With Eq. (4.41) this can also be written as
\begin{displaymath}
D_\alpha = \pi \tilde C_{vv ; \alpha\alpha}(0).
\end{displaymath} (4.22)

Computationally, the MSD is calculated using the Fast Correlation Algorithm (FCA) [51]. In this framework, in the discrete case, the mean-square displacement of a particle is given by

\begin{displaymath}
\Delta^2(m) = \frac{1}{N_t - m}\sum_{k=0}^{N_t-m-1}
[\textbf{r}(k+m) - \textbf{r}(k)]^2, \qquad m = 0\ldots N_t-1,
\end{displaymath} (4.23)

where $\textbf{r}(k)$ is the particle trajectory and $N_t$ is the the number of frames of the trajectory. We define now the auxiliary function
\begin{displaymath}
S(m) \doteq \sum_{k=0}^{N_t-m-1}[\textbf{r}(k+m) - \textbf{r}(k)]^2,
\qquad m = 0\ldots N_t-1,
\end{displaymath} (4.24)

which is splitted as follow:
$\displaystyle S(m)$ $\textstyle =$ $\displaystyle S_{AA+BB}(m) - 2 S_{AB}(m),$ (4.25)
$\displaystyle S_{AA+BB}(m)$ $\textstyle =$ $\displaystyle \sum_{k=0}^{N_t-m-1}[\textbf{r}^2(k+m) + \textbf{r}^2(k)],$ (4.26)
$\displaystyle S_{AB}(m)$ $\textstyle =$ $\displaystyle \sum_{k=0}^{N_t-m-1} \textbf{r}(k)\cdot\textbf{r}(k+m).$ (4.27)

The function $S_{AB}(m)$ can be computed using the FCA method described in Section A. For $S_{AA+BB}(m)$ the following recursion relation holds:
$\displaystyle S_{AA+BB}(m)$ $\textstyle =$ $\displaystyle S_{AA+BB}(m-1) - \textbf{r}^2(m-1) - \textbf{r}^2(N_t - m),$ (4.28)
$\displaystyle S_{AA+BB}(0)$ $\textstyle =$ $\displaystyle \sum_{k=0}^{N_t-1} \textbf{r}^2(k).$ (4.29)

This allows one to construct the following efficient scheme for the computation of the MSD:
  1. Compute $DSQ(k) = \textbf{r}^2(k),\qquad k = 0\ldots N_t-1;\qquad
DSQ(-1) = DSQ(N_t) = 0$.
  2. Compute $SUMSQ = 2\cdot\sum_{k=0}^{N_t-1} DSQ(k)$.
  3. Compute $S_{AB}(m)$ using the FFT method.
  4. Compute MSD(m) in the following loop:

    \begin{displaymath}
\begin{array}{ll}
&\begin{array}{ll}
SUMSQ &\leftarrow SU...
... \mbox{running from}\; 0 \;\mbox{to}\; N_t - 1\\
\end{array} \end{displaymath}

It should be noted that the efficiency of this algorithm is the same as for the FCA computation of time correlation functions since the number of operations in step (1), (2), and (4) grows linearly with $N_t$.
next up previous contents
Next: Parameters Up: Mean-Square Displacement Previous: Mean-Square Displacement   Contents
pellegrini eric 2009-10-06