md_trajectory = nlread('alumina_trajectory.nc')[-1]

msd = MeanSquareDisplacement(md_trajectory, atom_selection=Aluminum)

# Get the times in ps and the MSD values in Ang**2.
t        = msd.times().inUnitsOf(ps)
msd_data = msd.data().inUnitsOf(Angstrom**2)

# Plot the data using pylab.
import pylab

pylab.plot(t, msd_data, label='MSD of aluminum')
pylab.xlabel('t (ps)')
pylab.ylabel('MSD(t) (Ang**2)')
pylab.legend()

pylab.show()

# Fit the slope of the MSD to estimate the diffusion coefficient.
# If you discover non-linear behavior at small times, discard this initial part in the fit.
a = numpy.polyfit(t[5:], msd_data[5:], deg=1)

# Calculate the diffusion coefficient in Ang**2/ps.
diffusion_coefficient = a[0]/6.0