Algorithm to calculate diffusion coefficient

Question

According to the Einstein relationship, the diffusion coefficient $D$ is

$$\lim _{t\rightarrow \infty} \frac{\langle \left(\mathbf{r}(t)-\mathbf{r}(0) \right) ^2\rangle}{6t} = D$$

I have run a MD simulation with $N$ particles, and I have a file which has the location of each particle at every time step.

If I want to calculate $D$, I would have to do the following:

$$\langle (\mathbf{r}(t)-\mathbf{r}(0))^2\rangle = \frac{1}{N} \sum _{i=1}^N (\mathbf{r}_i(t)-\mathbf{r}_i(0))^2$$

Divide the above by $6t$, then take the limit as $t$ goes to infinity - which in my case is the final time step of my simulation.

My questions:

1. Is this the right way of calculating MSD?
2. If it is the right way of doing this, and my system reaches equilibrium and doesn't diverge, as $t\rightarrow \infty$, shouldn't my diffusion coefficient always be $0$, because of that infinity in the denominator?

Answer 1

Your approach, while mathematically correct, is not a good one for determining $D$. You're basically using information about the endpoints of the trajectories, and discarding all the information in between. A standard method is to evaluate $r^2/6t$ for each possible timestep $t$, i.e. all displacements over one frame, two frames, etc. The slope of $r^2$ vs. $t$ gives $D$. This does, however, suffer from the fact that these differences are correlated -- i.e. the two-frame displacements involve the same motions as the one-frame displacements.

A better approach is to use a correlation-based estimator, as described in this excellent paper:

Optimal estimation of diffusion coefficients from single-particle trajectories Christian L. Vestergaard, Paul C. Blainey, and Henrik Flyvbjerg Phys. Rev. E 89, 022726 – Published 28 February 2014 https://doi.org/10.1103/PhysRevE.89.022726