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

nearest_neighbor = NearestNeighbor(md_trajectory,
                                   start_time=10000.0*fs,
                                   end_time=50000.0*fs,
                                   cutoff_radius=2.5*Angstrom)

# Get the times in ps and the nearest neighbor values.
t       = nearest_neighbor.times().inUnitsOf(ps)
nn_data = nearest_neighbor.data()

# Plot the nearest neighbor number and its variance using pylab.
import pylab

pylab.plot(t, nn_data[:, 0], label='Number of nearest neighbors')
pylab.plot(t, nn_data[:, 1], label='Variance nearest neighbors')
pylab.xlabel('t (ps)')
pylab.ylabel('<NN>')
pylab.legend()

pylab.show()