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()