Calculating self-diffusion coefficient from MSD programmatically
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.