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. Thanks Yinhan # Physical background #------------------------------------- # Calculate the Hofstadter's Butterfly of G-BN from __future__ import division # so that 1/2==0.5, and not 0 from math import pi, sqrt,cos,sin from cmath import exp import numpy as np import kwant #For plotting from matplotlib import pyplot #Define the graphnee lattice sin_30, cos_30=(1 / 2, sqrt(3) / 2) #Define the scale a_tb=1.0 gamma=5 L=20*gamma W=5*gamma*sqrt(3) ab=4*(pi/sqrt(3))*gamma En=0.05 graphene=kwant.lattice.general([(1, 0),(sin_30, cos_30)],\ [(0, -1/(4*cos_30)),(0, 1/(4*cos_30))]) a, b = graphene.sublattices class SimpleNamespace(object): def __init__(self, **kwargs): self.__dict__.update(kwargs) def make_system(L, W, mu_lead=2): #### Define the scattering region. ########### # slab scattering region def slab_shape(pos): (x, y)=pos return 0 <= x <= L and 0 <= y <= W # Define on-site energy def potentiala(site): (x, y)=site.pos da1=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)*En def potentialb(site): (x, y)=site.pos da1=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)*En def onsitea(site,p): return potentiala(site)-p.mu def onsiteb(site,p): return potentialb(site)-p.mu sys = kwant.Builder() sys[a.shape(slab_shape,(1,1))]=onsitea sys[b.shape(slab_shape,(1,1))]=onsiteb # specify the hoppings of the graphene lattice in the # fomat expected by builder.HoppingKind def 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 result def 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 result def hopping1(site1,site2,p): (x1,y1)=site1.pos result=(-p.t+2*Ax(x1,y1)/3)\ *exp(1j * 2 * pi * 0.4 * x1/(gamma * gamma)) return result def hopping2(site1,site2,p): (x1,y1)=site2.pos result=-p.t+Ay(x1-0.25,y1-sqrt(3)/4)/sqrt(3)\ -Ax(x1-0.25,y1-sqrt(3)/4)/3 return result def hopping3(site1,site2,p): (x1,y1)=site2.pos result=-p.t-Ay(x1+0.25,y1-sqrt(3)/4)/sqrt(3)\ -Ax(x1+0.25,y1-sqrt(3)/4)/3 return result kind1=kwant.builder.HoppingKind((0,0),a,b) sys[kind1]=hopping1 kind2=kwant.builder.HoppingKind((0,1),a,b) sys[kind2]=hopping2 kind3=kwant.builder.HoppingKind((-1,1),a,b) sys[kind3]=hopping3 #### Define the leads. #### # left lead0 sym0 = kwant.TranslationalSymmetry(graphene.vec((1, 0))) def lead0_shape(pos): x, y = pos return (W/2-2.5*gamma <= y <= W/2+2.5*gamma) def lead_hopping(site1,site2,p): return -p.t def lead_onsite(site1,p): return mu_lead lead0 = 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_onsite lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] =lead_hopping #sys.attach_lead(lead0) # down lead1 sym1 = kwant.TranslationalSymmetry(graphene.vec((1, -2))) def lead1_shape(pos): x, y= pos return (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_onsite lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = lead_hopping #sys.attach_lead(lead1) # right lead2 lead2=lead0.reversed() #sys.attach_lead(lead2) # right-up lead3 sym3= kwant.TranslationalSymmetry(graphene.vec((-1,2))) def lead3_shape(pos): x, y = pos return (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_onsite lead3[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]]= lead_hopping #sys.attach_lead(lead3) # left-up lead4 lead4=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 object r = 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 above sigma_xx = E_x / (E_x**2 + E_y**2) sigma_xy = E_y / (E_x**2 + E_y**2) return sigma_xx, sigma_xy def plot_Hallconductance(sys,p): # Compute conductance Bs = 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]).T pyplot.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 system sys=sys.finalized() # Plot conductance # Compute conductance plot_Hallconductance(sys,p) if __name__ == '__main__': main()