import pylab

# Read in a list with the trajectories
trajectories = nlread('trajectory.nc', MDTrajectory)

pylab.figure()
for trajectory in trajectories:
    # Get the time axis of the trajectory
    t = trajectory.times().inUnitsOf(femtoSecond)

    # Calculate the time step
    time_step = t[-1]/trajectory.steps()[-1]

    # Get the energies
    epot = trajectory.potentialEnergies().inUnitsOf(eV)
    ekin = trajectory.kineticEnergies().inUnitsOf(eV)

    # Plot the data
    pylab.plot(t,epot+ekin, label='timestep=' + str(time_step)+' fs')
    pylab.xlabel('time (fs)')
    pylab.ylabel('total energy (eV)')
    pylab.legend()

pylab.show()