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

v_distribution = VelocityDistribution(md_trajectory)

# Get the bin_centers and the histogram of the velocity distribution.
velocities = v_distribution.velocities().inUnitsOf(Ang/fs)
histogram  = v_distribution.data()

# Plot the data using pylab.
import pylab

pylab.bar(velocities, histogram, label='Velocity distribution', width=0.001)
pylab.xlabel('v (Ang/fs)')
pylab.ylabel('Histogram')
pylab.legend()

pylab.show()