Calculating self-diffusion coefficient from MSD programmatically

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

In 1D diffusion along the $x$ axis by Brownian motion, let the mean jump lengths be $a$ each taking a time $\tau$ and this should be large compared to the time between successive collisions, which are $\approx 10$ fs, i.e. $a$ is the sum of all displacements during time $\tau$. Assuming this means that there is no correlation between collisions at time $t$ and $t+\tau$. From this the diffusion equation (Ficks Law) can be derived and the diffusion coefficient defined by $D_x=\langle \Delta x^2\rangle/2\tau$ where $\langle \Delta x^2\rangle$ is the variance of the jump lengths $a$ or $\langle\Delta x^2\rangle = ( x-\langle x\rangle)^2$. In 3D as there is no correlation between the motions on each axis, $D=\langle \Delta x^2\rangle/6\tau$. This equation is subtly different from your equation and so you now need to connect these terms with what you have calculated. You should find that your diffusion coefficient settles down to a constant value (plus or minus statistical fluctuations) as the calculation proceeds.