Dear Yinhan, The reason that you obtain 0 is because there are no modes open at the energy at which you are calculating the scattering matrix. When I take your code and plot the bandstructure for the different leads, I see that the bands are in a range of energies between 1.9 and 2.1, however you calculate the scattering matrix at 0 energy (you do not provide an explicit energy value, and the default value is 0). To solve your problem you should either select the energy that you want by providing an explicit value to `kwant.smatrix`, or modify your model so that the bands are open at 0 energy. Happy Kwanting, Joe On 2 June 2016 at 18:36, 张银寒 <zhangyinhan2008@gmail.com> wrote:

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