Dear Sir,I am trying to reproduce the Hall conductance in a Graphene-hexagonal boron-nitride structure by using Kwant, referring to (https://arxiv.org/abs/1401.4401).The five-terminal Hall bar geometry is considered to study the Hall conductance. My problem is that the smatrix.transmission always vanishes in my code.How can I fix the problem? I pasted my KWANT code. Any suggestions will be greatly appreciated.ThanksYinhan# Physical background#-------------------------------------# Calculate the Hofstadter's Butterfly of G-BNfrom __future__ import division # so that 1/2==0.5, and not 0from math import pi, sqrt,cos,sinfrom cmath import expimport numpy as npimport kwant#For plottingfrom matplotlib import pyplot#Define the graphnee latticesin_30, cos_30=(1 / 2, sqrt(3) / 2)#Define the scalea_tb=1.0gamma=5L=20*gammaW=5*gamma*sqrt(3)ab=4*(pi/sqrt(3))*gammaEn=0.05graphene=kwant.lattice.general([(1, 0),(sin_30, cos_30)],\[(0, -1/(4*cos_30)),(0, 1/(4*cos_30))])a, b = graphene.sublatticesclass SimpleNamespace(object):def __init__(self, **kwargs):self.__dict__.update(kwargs)def make_system(L, W, mu_lead=2):#### Define the scattering region. ############ slab scattering regiondef slab_shape(pos):(x, y)=posreturn 0 <= x <= L and 0 <= y <= W# Define on-site energydef potentiala(site):(x, y)=site.posda1=cos((-sqrt(3)*x/2+y/2)*ab)+cos(-y*ab)+\cos((sqrt(3)*x/2+y/2)*ab)da2=sin((-sqrt(3)*x/2+y/2)*ab)+sin(-y*ab)+\sin((sqrt(3)*x/2+y/2)*ab)return (2*pi/gamma)*(da1-sqrt(3)*da2)*Endef potentialb(site):(x, y)=site.posda1=cos((-sqrt(3)*x/2+y/2)*ab)+cos(-y*ab)+\cos((sqrt(3)*x/2+y/2)*ab)da2=sin((-sqrt(3)*x/2+y/2)*ab)+sin(-y*ab)+\sin((sqrt(3)*x/2+y/2)*ab)return (2*pi/gamma)*(da1+sqrt(3)*da2)*Endef onsitea(site,p):return potentiala(site)-p.mudef onsiteb(site,p):return potentialb(site)-p.musys = kwant.Builder()sys[a.shape(slab_shape,(1,1))]=onsiteasys[b.shape(slab_shape,(1,1))]=onsiteb# specify the hoppings of the graphene lattice in the# fomat expected by builder.HoppingKinddef Ax(x,y):result=-En*2*(pi/gamma)*(cos((-sqrt(3)*x/2+y/2)*ab)+\cos((sqrt(3)*x/2+y/2)*ab)-2*cos(-y*ab))return resultdef Ay(x,y):result=-En*2*sqrt(3)*(pi/gamma)*(cos((-sqrt(3)*x/2+y/2)*ab)\-cos((sqrt(3)*x/2+y/2)*ab))return resultdef hopping1(site1,site2,p):(x1,y1)=site1.posresult=(-p.t+2*Ax(x1,y1)/3)\*exp(1j * 2 * pi * 0.4 * x1/(gamma * gamma))return resultdef hopping2(site1,site2,p):(x1,y1)=site2.posresult=-p.t+Ay(x1-0.25,y1-sqrt(3)/4)/sqrt(3)\-Ax(x1-0.25,y1-sqrt(3)/4)/3return resultdef hopping3(site1,site2,p):(x1,y1)=site2.posresult=-p.t-Ay(x1+0.25,y1-sqrt(3)/4)/sqrt(3)\-Ax(x1+0.25,y1-sqrt(3)/4)/3return resultkind1=kwant.builder.HoppingKind((0,0),a,b)sys[kind1]=hopping1kind2=kwant.builder.HoppingKind((0,1),a,b)sys[kind2]=hopping2kind3=kwant.builder.HoppingKind((-1,1),a,b)sys[kind3]=hopping3#### Define the leads. ##### left lead0sym0 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))def lead0_shape(pos):x, y = posreturn (W/2-2.5*gamma <= y <= W/2+2.5*gamma)def lead_hopping(site1,site2,p):return -p.tdef lead_onsite(site1,p):return mu_leadlead0 = kwant.Builder(sym0)hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))lead0[graphene.shape(lead0_shape, (0, W/2))] = lead_onsitelead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] =lead_hopping#sys.attach_lead(lead0)# down lead1sym1 = kwant.TranslationalSymmetry(graphene.vec((1, -2)))def lead1_shape(pos):x, y= posreturn (L/4-2.5*gamma <= x <= L/4 + 2.5*gamma )lead1 =kwant.Builder(sym1)lead1[graphene.shape(lead1_shape, (L/4-2.5*gamma+1,0))] = lead_onsitelead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = lead_hopping#sys.attach_lead(lead1)# right lead2lead2=lead0.reversed()#sys.attach_lead(lead2)# right-up lead3sym3= kwant.TranslationalSymmetry(graphene.vec((-1,2)))def lead3_shape(pos):x, y = posreturn (0.75*L-2.5*gamma <= x <= 0.75*L +2.5*gamma)lead3=kwant.Builder(sym3)lead3[graphene.shape(lead3_shape,(0.75*L-2.5*gamma+1,0))] = lead_onsitelead3[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]]= lead_hopping#sys.attach_lead(lead3)# left-up lead4lead4=lead1.reversed()#sys.attach_lead(lead4)#print("ab",ab)return sys, [lead0,lead1,lead2,lead3,lead4]def calculate_sigmas(G):# reduce by one dimension G -> G[temp, temp]temp = [0, 1, 3, 4]G = G[temp, :]G = G[:, temp]# invert R = G^-1# find out whether it is a numpy objectr = np.linalg.inv(G)# Voltages follow: V = R I[temp]V = r.dot(np.array([1, -1, 0, 0]))# Completely solved ihe six terminal system.# Consider the 2x2 conductance now: Use I = sigma U#E_x = V[1] - V[0]E_x=(V[3]-V[2])#print("V=",V)#print("V[0]=",V[0],"V[1]=",V[1],"V[2]=",V[2],"V[3]=",V[3])#E_x=(V[3]-V[2])/L#E_y = V[1] - V[3]E_y=V[3]#E_y=V[3]/W# formula abovesigma_xx = E_x / (E_x**2 + E_y**2)sigma_xy = E_y / (E_x**2 + E_y**2)return sigma_xx, sigma_xydef plot_Hallconductance(sys,p):# Compute conductanceBs = np.linspace(0.5, 1, 10)num_leads=len(sys.leads)def G(sys, p):smatrix = kwant.smatrix(sys, args=[p])print("trans",smatrix.transmission(0,1))Bs = np.linspace(0.5, 1, 10)G = [[smatrix.transmission(i, j) for i in range(num_leads)] for j in range(num_leads)]print(G)G -= np.diag(np.sum(G, axis=0))return calculate_sigmas(G)sigmasxx, sigmasxy = np.array([G(sys, p) for p.B in Bs]).Tpyplot.figure()pyplot.plot(Bs,sigmasxy)pyplot.xlabel("B [t]")pyplot.ylabel("sigma_{xy} [e^2/h]")pyplot.show()def main():t0=2/(sqrt(3)*ab)p=SimpleNamespace(t=t0, mu=0.25, B=None)sys, leads = make_system(L, W, mu_lead=2)# To highlight the two sublattices of graphene, we plot one with# a filled, and the other one with an open circle:def family_colors(site):return 0 if site.family == a else 1# Then, plot the system with leads.#kwant.plot(sys, site_color=family_colors, site_lw=0.1,# lead_site_lw=0, colorbar=False)# Attach the leads to the system.for lead in leads:sys.attach_lead(lead)# Then, plot the system with leads.#kwant.plot(sys, site_color=family_colors, site_lw=0.1,# lead_site_lw=0, colorbar=False)# Finized the systemsys=sys.finalized()# Plot conductance# Compute conductanceplot_Hallconductance(sys,p)if __name__ == '__main__':main()