Calculating diffusion coefficient from Mean Squared Displacement

Question

I am performing a molecular dynamics simulation and I have calculated the velocity autocorrelation function (VACF), $\psi (t)$ and mean square displacement, $\langle (\Delta r)^2 \rangle$. In my simulation, I have $N$ particles interacting with the Lennard-Jones potential function, and I am stepping forward in time from $0$ to $T$, in steps of $\Delta t$.

From here, I see that $$ D = v_0^2 \int_0^{\infty} \psi(t) \mathrm{d}t$$

where $v_0$ is the thermal velocity.

My question is, how do I calculate $D$ if I only know values of $v_0$ and $\psi(t)$ for $t$ in $[0,T]$ (rather than up to $t=\infty$)? Is there a better way to do this?

Any advice would be appreciated.

Answer 1

As with all MD simulations, you have to assume (often wrongly) convergence with finite time. This is fairly easy to do with autocorrelation functions though, because you know that once they become negative, you are already in the random fluctuation (i.e. uncorrelated noise) zone and this is the point where you can stop. This is slightly trickier in some cases still, especially for unthermostatted systems, where you can have deterministic correlated oscillatory behaviour around zero, so this metric is not very good and I am not sure there is a better alternative than just visually inspecting the plot to pick a suitable cutoff. This won't be a problem, since you will typically sample for a longer time than your decorrelation time, mostly because you need multiple datapoints for autocorrelation function convergence. Finally, in some really pathological cases (e.g. in 1D systems), you can get decay to a constant value which is substantially larger than 0. Thankfully, this behaviour isn't very common at all but it is still worth knowing that it is possible (and in these cases I don't even know how one calculates a diffusion coefficient or whether doing so makes any sense in the first place).

The problem with autocorrelation functions is that this oscillation (even the random one) can still impact your results depending on where you decide to stop and that's why I would say that calculating diffusion coefficients with mean-square displacement is much more robust, because you usually get an almost perfect straight line, whose slope you can accurately obtain from linear regression and this is pretty robust to noise, number of datapoints or cutoff choice (!), certainly more so than integrating the autocorrelation function. In cases where you don't get a straight line (e.g. superdiffusion or subdiffusion), this actually tells you that there is something else going on and calculating diffusion coefficients might need additional assumptions for this system. In any case, you get more information using the RMSD method.

So the bottom line is: if you have a choice between integrating and not integrating, (almost) always choose a method that doesn't rely on integration and if this is for coursework, don't be too worried if the VACF method gives you noisier results (although it's always a good idea to be completely sure you aren't doing anything wrong).

Answer 2

Translation diffusion coefficient $D_{t}$ may be calculated by linear fitting of the MSD of center of mass of a molecule:

$$MSD(\tau)=<(\textbf{r}(t+\tau)-\textbf{r}(t))^{2}>$$

$$D_{t}=\frac{MSD(\tau)}{6\tau}$$

where $\tau$ is the lag time between the two positions. One should also be aware of the artifacts from the periodic boundary conditions (is they are imposed).

Rotational diffusion coefficient $D_{r}$ may be estimated through fitting of the correlation functions and estimation of the rotational correlation time $\tau$:

$$\tau=S_{R}^{2}\tau_{R_{s}}+(1-S_{R}^{2})\tau_{R_{f}}$$

with $\tau_{R_{s}}$ and $\tau_{R_{f}}$ being the slow and fast correlation times respectively; $S_{R}^{2}$ is a parameter in fitting of the correlation function. The $D_{r}$ acquires following form:

$$D_{r}=\frac{1}{6\tau}$$

Examples of these calculations can be found here:

Slow-Down in Diffusion in Crowded Protein Solutions Correlates with Transient Cluster Formation

Evaluating Rotational Diffusion from Protein MD Simulations