# Restart from the previous calculation.
device_configuration = nlread('zero-bias.hdf5', object_id='gID000')[0]

# Array with k-values.
numbers_of_kpoints = range(3,40,2)

# Energies at which to calculate the transmission function.
energies = [-3,0,3]

# Array with transmission values.
T = numpy.zeros([len(energies), len(numbers_of_kpoints)])

# Loop over numbers of k-points.
for i, num_k in enumerate(numbers_of_kpoints):        
    # -------------------------------------------------------------
    # Transmission spectrum
    # -------------------------------------------------------------
    transmission_spectrum = TransmissionSpectrum(
        configuration=device_configuration,
        energies=numpy.array(energies)*eV,
        kpoints=MonkhorstPackGrid(num_k,num_k),
        energy_zero_parameter=AverageFermiLevel,
        infinitesimal=1e-06*eV,
        self_energy_calculator=KrylovSelfEnergy(),
        )

    # Put the k-point averaged transmission into the array.
    T[:,i] = transmission_spectrum.evaluate()

# Convert to numpy array.
numbers_of_kpoints  = numpy.array(numbers_of_kpoints)

# Plot the results.
import pylab

pylab.subplot(311)
pylab.plot(numbers_of_kpoints , T[0,:] ,'ro-',label='E=-3 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()

pylab.subplot(312)
pylab.plot(numbers_of_kpoints , T[1,:] ,'ro-',label='E=0 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()

pylab.subplot(313)
pylab.plot(numbers_of_kpoints , T[2,:] ,'ro-',label='E=3 eV')
pylab.xlabel('# k-points',fontsize=16)
pylab.ylabel('Transmission',fontsize=16)
pylab.legend()

# Show the plot.
pylab.show()