--------------------------------------------------------------------------------
Diffusion coefficient calculations

D = <(delta r)^2>/(6 t)

Use qbox_msd.py <species> md1.r md2.r ... > msd.dat

Use the fit function of gnuplot to fit a*x+b
<(delta r)^2> = a*t + b

a = slope of MSD in [Bohr^2/step]

Conversion to cm2/s:
D [cm2/s] = a [Bohr^2/step] * 0.192946 / dt [a.u.]

The constant 0.192946 is (41341/6)*(0.529177)^2
where 1 ps = 41341 au_t and 1 Bohr = 0.529177 Ang
--------------------------------------------------------------------------------