md_trajectory = nlread('ag-md.hdf5', class_type=MDTrajectory)[-1]

# Get the vibrational density of states.
vdos = VibrationalDensityOfStates(md_trajectory)

# Get the frequencies associated with the DOS.
frequencies = vdos.frequencies()

# Convert the frequencies to energies in meV.
energies = (frequencies*planck_constant).inUnitsOf(meV)

# Plot the data using pylab.
import pylab

pylab.plot(energies, vdos.vibrationalDensityOfStates(), label='Vibrational density of states')
pylab.xlabel('E (meV)')
pylab.ylabel('S(E)')
pylab.xlim(0,50)
pylab.legend()

pylab.show()