How to calculate diffusion coefficient from MSD graph using GROMACS?

Question

In GROMACS I have used the command below to calculate mean square displacement (MSD) over time,

gmx msd -f ../run01/run01.trr -s ../run01/run01.tpr -s ../setup/index.ndx -o msd01.xvg -tu ns -type z

From the equation,

$$D = \frac{1}{6t} \langle(r(t) - r(0))^{2}\rangle$$

I suppose that in GROMACS the default resetting time for $r(0)$ is $\pu{10 ps}$ thus $t=\pu{10 ps}$. However my most confusing part is selecting the correct linear region. Is there proper mathematical method to find correct linear region?

In GROMACS, the args -beginfit -endfit defaulted to fitting the curve between 10% and 90%. That means I should select linear region as between 10% to 90% along the x axis?

An example graph,

Answer 1

MSD accuracy decreases as a function of the lag-time $t$, because for large $t$ there are fewer pairs of time slices of interval $t$ to average over. Think of the limit in calculating MSD for the largest possible lag-time interval of a simulation. In that case you can only use the initial and final frame to calculate the MSD at that $t$. Whereas for small $t$, you can average this over many sets of nearby time frames.

The above is why, beyond $\ce{15 ns}$, your MSD curve becomes highly non-linear. You simply don't have good statistics here. Very short times are also not ideal to use for $D$, as displacements are ballistic prior to a particle's first collision and the MSD slope is usually very steep at short $t$. Hence, choose the large smooth region in the middle of your MSD plot to calculate the slope. As long as your simulation is not changing it's diffusion properties over time, and you have sufficient data to average over, this region of the MSD should have a constant slope.

In your case, it looks like you either don't quite have enough data to generate a perfectly smooth MSD beyond $\ce{6 ns}$, or your simulation may be changing in some way over time. However, the medium-time MSD is smooth enough to calculate a diffusion constant. Just use the region from about $\ce{1 ns}$ to $\ce{5 ns}$.

By the way, you can improve MSD statistics by using a shorter reset time (this also increases computational cost).

Answer 2

If you can obtain the raw data used in that plot, as a purely statistical and numerical exercise you can determine the interval which will give you minimum uncertainty in the estimate of the gradient.

A short line segment will inherently have higher uncertainty in estimated gradient simply because you are using fewer sample points. But, for your data, a too-long line segment will also begin to have higher uncertainty in estimated gradient because the segment will start including clearly non-linear, noisy data. This implies that there will be some optimum segment of the line to use. Find that optimum and use it, and as a byproduct you also get an estimate of the uncertainty for free.

However, this is only correct as a matter of abstract statistics. The uncertainty will be drastically underestimated (relative to the ground truth of "if I ran a completely independent simulation would I get a similar gradient"), because molecular dynamics snapshots are correlated instead of being statistically independent. If you need to be very accurate in stating the uncertainty (and you may not need to) then you should be referring to this recent paper and choosing a good methodology from it: Optimal estimates of self-diffusion coefficients from molecular dynamics simulations. Note that the paper authors have made available a python script which implements their recommendation, the generalised least-squares approach.