Dear all, I need to access the matrix diagonal element in the scattering region in my program. Is there any way to access the matrix diagonal element in the following program or do I have to use numpy? Best regard, import kwant import matplotlib.pyplot as plt import tinyarray import numpy as np import math from functools import partial from kwant.physics import dispersion #mport matplotlib d=1.42; a1=d*math.sqrt(3); latt = kwant.lattice.general([(a1,0),(0.5*a1,a1*math.sqrt(3)/2)], [(0,0),(0,a1/math.sqrt(3))]) a,b = latt.sublattices syst= kwant.Builder() #................................................................................... def rectangle(pos): x, y = pos z=x**2+y**2 return -4*3*a1<y<4.5*3*a1 and -4.5*d<=x<4.5*d syst[a.shape(rectangle, (2,2))] = 0 syst[b.shape(rectangle, (2,2))] = 0 #nearest neighbors............................................................. syst[[kwant.builder.HoppingKind((0,0),a,b)]] = -2.7 syst[[kwant.builder.HoppingKind((0,1),a,b)]] = -2.7 syst[[kwant.builder.HoppingKind((-1,1),a,b)]] = -2.7 ax=kwant.plot(syst); sym = kwant.TranslationalSymmetry(latt.vec((-1,0))) sym.add_site_family(latt.sublattices[0], other_vectors=[(-1, 2)]) sym.add_site_family(latt.sublattices[1], other_vectors=[(-1, 2)]) lead = kwant.Builder(sym) def lead_shape(pos): x, y = pos return -4*3*a1<y<4.5*3*a1 lead[a.shape(rectangle, (2,2))] = 0 lead[b.shape(rectangle, (2,2))] = 0 #nearest neighbors............................................................. lead[[kwant.builder.HoppingKind((0,0),a,b)]] = -2.7 lead[[kwant.builder.HoppingKind((0,1),a,b)]] = -2.7 lead[[kwant.builder.HoppingKind((-1,1),a,b)]] = -2.7 syst.attach_lead(lead) syst.attach_lead(lead.reversed()) ax=kwant.plot(syst); ###############################################################################\ fsys = syst.finalized() Sites= list(fsys.sites) #list of all the sites in the scattering region number_of_sites=len(Sites) bands = dispersion.Bands(fsys.leads[0]) momenta = np.linspace(-np.pi,np.pi, 100) energies=[bands(k) for k in momenta] np.savetxt('band1.txt',bands(-np.pi)) plt.plot(momenta,energies) plt.xlabel(r'$K$') plt.ylabel(r'$band structure (eV)$') plt.ylim((-0.5,0.5)) plt.savefig('bandsarmBiStrain.pdf') plt.show()