Sorry, I just realized I made a silly mistake. I should replace sym_lead_vertical = kwant.TranslationalSymmetry((0,-1,0)) by sym_lead_vertical = kwant.TranslationalSymmetry((0,1,0)) although I don't understand why it is the case. I apologize for inconvenience caused. Best, T.C. On 20 July 2015 at 00:03, T.C. Wu <tcwu@ust.hk> wrote:

Dear all,

I am a beginner in kwant and I am trying to reproduce the quantum Hall effect demostrated in

http://nbviewer.ipython.org/github/topocm/topocm_content/blob/master/w3_pump...

But in my case, I modified the sample code a little bit such that a 3D pill-box like system with thickness of 1 atomic layer is built (see attached figure). I changed nothing but replaced (xs,ys) by (xs,ys,zs) and so on. However, the Hall conductance simulated in my system is completely different from the 2D case ( no plateaus formed, see the attached figure).

Could you please give me some advice?

# My code#

import numpy as np import matplotlib import kwant from matplotlib import pyplot as plt

L = 100 W = 60 H = 1 w_lead = W-2 w_vert_lead = L/2-1 t = 1.0 mu = 0.3 mu_lead = 0.3 Bs = np.linspace(.02, .15, 200)

class SimpleNamespace(object): def __init__(self, **kwargs): self.__dict__.update(kwargs)

## Quantum hall bar codes def qhe_hall_bar(L=50, W=10,H=2, w_lead=10, w_vert_lead=None): """Create a hall bar system.

Square lattice, one orbital per site. Returns finalized kwant system.

Arguments required in onsite/hoppings: t, mu, mu_lead """

L = 2*(L//2) W = 2*(W//2) H = 2*(H//2) w_lead = 2*(w_lead//2) if w_vert_lead is None: w_vert_lead = w_lead else: w_vert_lead = 2*(w_vert_lead//2)

# bar shape def bar(pos): (x, y,z) = pos return (x >= -L/2 and x <= L/2) and (y >= -W/2 and y <= W/2) and (z>= -H/2 and z <= H/2)

# Onsite and hoppings def onsite(site, par): (x,y,z) = site.pos return 4*par.t - par.mu

def hopping_Ax(target, source, par): xt, yt, zt = target.pos xs, ys, zs = source.pos return -par.t * np.exp(-0.5j * par.B * (xt + xs) * (yt - ys))

def make_lead_hop_y(x0): def hopping_Ay(target, source, par): print (x0) xt, yt, zt = target.pos xs, ys, zs = source.pos return -par.t * np.exp(-1j * par.B * x0 * (yt - ys)) return hopping_Ay

def lead_hop_vert(target, source, par): return -par.t

# Building system lat = kwant.lattice.general(np.identity(3)) sys = kwant.Builder() print('sys built')

sys[lat.shape(bar, (0, 0,0))] = onsite sys[lat.neighbors()] = hopping_Ax print('sys hopping done') #kwant.plot(sys) # Attaching leads sym_lead = kwant.TranslationalSymmetry((-1, 0,0))

def lead_shape(pos): (x, y,z) = pos print('lead_shape') return (-w_lead / 2 <= y <= w_lead / 2) and (z>= -H/2 and z <= H/2)

lead_onsite = lambda site, par: 4 * par.t - par.mu_lead

sym_lead_vertical = kwant.TranslationalSymmetry((0,-1,0)) lead_vertical1 = kwant.Builder(sym_lead_vertical) lead_vertical2 = kwant.Builder(sym_lead_vertical) print('vert lead done') def lead_shape_vertical1(pos): (x, y,z) = pos return (-L/4 - w_vert_lead/2 <= x <= -L/4 + w_vert_lead/2) and (H/2 >= z>= -H/2)

def lead_shape_vertical2(pos): (x, y,z) = pos return (+L/4 - w_vert_lead/2 <= x <= +L/4 + w_vert_lead/2) and (H/2 >= z>= -H/2)

lead_vertical1[lat.shape(lead_shape_vertical1, (-L/4 , 0,0))] = lead_onsite print('onsite') lead_vertical1[lat.neighbors()] = lead_hop_vert lead_vertical2[lat.shape(lead_shape_vertical2, (L/4,0, 0))] = lead_onsite lead_vertical2[lat.neighbors()] = lead_hop_vert

print('lead defined')

sys.attach_lead(lead_vertical1) sys.attach_lead(lead_vertical2)

sys.attach_lead(lead_vertical1.reversed()) sys.attach_lead(lead_vertical2.reversed()) lead = kwant.Builder(sym_lead) lead[lat.shape(lead_shape, (-1, 0,0))] = lead_onsite lead[lat.neighbors()] = make_lead_hop_y(-L/2)

sys.attach_lead(lead)

lead = kwant.Builder(sym_lead) lead[lat.shape(lead_shape, (-1, 0,0))] = lead_onsite lead[lat.neighbors()] = make_lead_hop_y(L/2)

sys.attach_lead(lead.reversed()) print('lead attached') return sys.finalized()

def calculate_sigmas(G): # smat = smat.data # smat = abs(smat)**2 # reduce by one dimension G -> G[temp, temp] temp = [0, 1, 3, 4, 5] 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([0, 0, 0, -1, 1])) # Completely solved the six terminal system. # Consider the 2x2 conductance now: Use I = sigma U E_x = V[1] - V[0] E_y = V[1] - V[3] #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 ,-1*E_x,-1*E_y print('hi') par = SimpleNamespace(t=t, mu=mu, mu_lead=mu_lead, B=None) sys = qhe_hall_bar(L=L, W=W,H=H ,w_lead=w_lead, w_vert_lead=w_vert_lead)

kwant.plot(sys) sigmasxx = [] sigmasxy = [] rxx = [] rxy = [] #Bs = np.linspace(.02, .15, 200)

for par.B in Bs: s = kwant.smatrix(sys, args=[par]) G = np.array([[s.transmission(i, j) for i in range(len(sys.leads))] for j in range(len(sys.leads))]) G -= np.diag(np.sum(G, 0)) sigma = calculate_sigmas(G) sigmasxx.append(sigma[0]) sigmasxy.append(sigma[1]) rxx.append(sigma[2]) rxy.append(sigma[3])

fig = plt.figure() ax = plt.subplot(1, 1, 1)

ax.plot(1/Bs, sigmasxx, 'k') ax.plot(1/Bs, sigmasxy, 'r'); ax.set_xticks([]) ax.set_xlabel('$B^{-1} [a.u.]$') ax.set_yticks(range(4)) ax.set_yticklabels(['${0}$'.format(i) for i in range(4)]) ax.set_ylabel('$\sigma_{xx}, \sigma_{xy}\,[e^2/h]$') plt.show()

plt.plot(Bs,rxx) plt.plot(Bs,rxy) plt.show() #ax.set_ylim([-.1, 3.3]);