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

q_number = QNumbers(
    md_trajectory,
    element_selection=[Oxygen, Silicon, Oxygen],
    fuzz_factor=1.2,
    calculate_distribution=False,
)


# Get the histogram and the associated the coordination numbers.
times = q_number.times().inUnitsOf(fs)
q_number_evolution = q_number.data()

# Plot the data using pylab.
import pylab

# Make a plot for each q-number.
n_q_numbers = q_number_evolution.shape[1]
for i in range(n_q_numbers):
    label = 'Q%i O-Si-O' % i
    pylab.plot(times, q_number_evolution[:, i], label=label)
    pylab.xlabel('Time (fs)')
    pylab.ylabel('Q-number count')


pylab.legend()

pylab.show()
