Dear Sahu,
What you see is just a flat band. It doesn't mean that there is an error in
the code.
Actually, it is very common to see it in graphene for example if you have
dangling sites in an infinite stripe.
You need of course to avoid that energy in the calculation of the
scattering matrix because it is singular.
I hope this helps,
Adel
On Fri, Oct 14, 2022 at 11:26 AM <sahu.aji...@gmail.com> wrote:
> Dear all,
> I make AGNR in Kwant and introduce the defect in both the scattering
> region and lead.
> When plotting the bandstructure of the lead there is a straight line that
> comes at E=0.
> I am confused that there is some error in the code or not. Please help
>
> Thanking You
> code-
> ________________
> import kwant
> import numpy as np
> import matplotlib.pyplot as plt
>
> pi = np.pi
> a = 1
> d = 1.42;
> a1 = d*np.sqrt(3);
> lat = kwant.lattice.general(
> [[a1, 0], [a1/2, a1*np.sqrt(3)/2]],[[0, 0], [0, a1/np.sqrt(3)]]
> )
> a,b = lat.sublattices
>
>
> def set_pub(): # for plotting the figures
> plt.rcParams.update({
> "font.weight": "bold", # bold fonts
> "font.size":18,
> #"tick.labelsize": 15, # large tick labels
> "lines.linewidth": 2, # thick lines
> "lines.color": "k", # black lines
> "grid.color": "0.5", # gray gridlines
> "grid.linestyle": "-", # solid gridlines
> "grid.linewidth": 0.5, # thin gridlines
> "xtick.major.size":10,
> "xtick.minor.size":5,
> "ytick.major.size":10,
> "ytick.minor.size":5,
> "savefig.dpi": 600, # higher resolution output.
> })
>
>
> def make_sys(w,l):
> def rect(pos):
> x,y = pos
> return -w*a1<x<w*a1 and -l*a1<y<l*a1
> model1 = kwant.Builder()
> model1[lat.shape(rect,(1,1))] = 0
> model1[lat.neighbors()] = 2.46
> return model1
>
>
> def lead_join(model1,w):
> def lead_shape(pos):
> x, y = pos
> return -w*a1<x<w*a1
> Num = 8*(a1*np.sqrt(3)/2)
> sym = kwant.TranslationalSymmetry([0,Num])
> lead = kwant.Builder(sym)
> lead[lat.shape(lead_shape, (1,1))] = 0
> lead[lat.neighbors()] = 2.46
> model1.attach_lead(lead)
> model1.attach_lead(lead.reversed())
> return model1,lead
>
>
> def defect(model1,lead):
> def delet_lead(pos, shift = +a1*np.sqrt(3)):
> x, y = pos
> z=x**2+(y-shift)**2
> return z<(a1)**2
> def delet_str(pos, shift = a1*np.sqrt(3)):
> x, y = pos
> z = x**2 + (y-shift)**2
> return z<(a1)**2#z<(2*a1)**2
> del model1[lat.shape(delet_str, (1,1))]
> del lead[lat.shape(delet_lead, (2,2))]
> kwant.plot(model1)
> set_pub()
> plt.rc('axes', linewidth = 2.0)
> plt.show()
> return model1
>
>
> def plot_bandstructure(flead, momenta):
> bands = kwant.physics.Bands(flead)
> energies = [bands(k) for k in momenta]
> band_gap = 100
> for k in momenta:
> for band in bands(k):
> if(band <0):
> #Lower than fermi band
> min_band=band
> else:
> max_band=band
> if(max_band-min_band<band_gap):
> band_gap=max_band-min_band
> break
> print("Gap Energy= " + str(band_gap) +" ev ")
> kwant.plotter.bands(flead) # bandstructure of lead
> plt.ylim([-3,3])
> set_pub()
> plt.xlabel("K $(\AA^{-1})$",fontweight='bold')
> plt.ylabel("Energy (eV)",fontweight='bold')
> plt.show()
> return band_gap
>
>
> def main():
> w = 5
> l = 8
> sys = make_sys(w,l)
> syst,lead = lead_join(sys,w)
> syst = defect(syst,lead)
> fsys = syst.finalized()
> momenta = [-pi + 0.02 * pi * i for i in range(101)]
> Band_Gap = plot_bandstructure(lead.finalized(), momenta)
>
>
> if __name__ == '__main__':
> main()
>
--
Abbout Adel
