Hello, I am trying to simulate the quantum hall effect in a hall bar made of a honeycomb (graphene) lattice. I modified the code submitted by C. Groth as a response to a 2015 thread (link below). As soon as I add the vertical leads, I cant seem to get the expected phenomenon for both longitudinal and hall resistances. Here is my code, if anyone spots an error or a solution, do let me know! Sid Sule -------------------------------------------- from cmath import exp import numpy from matplotlib import pyplot import kwant from kwant.digest import gauss def hopping(sitei, sitej, phi, salt): xi, yi = sitei.pos xj, yj = sitej.pos return -exp(-0.5j * phi * (xi - xj) * (yi + yj)) def onsite(site, phi, salt): return 0.05 * gauss(repr(site), salt) def make_system(L=20, W=30): def central_region(pos): x, y = pos return abs(x) < (W/2) and abs(y) < (L/2) lat = kwant.lattice.honeycomb() sys = kwant.Builder() sys[lat.shape(central_region, (0, 0))] = onsite sys[lat.neighbors()] = hopping vertices = numpy.array([[-15, -6, 5, 10], [-15, -6, -10, -5], [6, 15, 5, 10], [6, 15, -10, -5], [-2, 2, 5, 10], [-2, 2, -10, -5]]) for v in vertices: def in_hole(site): x, y = site.pos return (v[0] <= x <= v[1]) and (v[2] <= y <= v[3]) for site in filter(in_hole, list(sys.sites())): del sys[site] sym = kwant.TranslationalSymmetry((-1, 0)) lead = kwant.Builder(sym) lead[lat.shape(lambda s: abs(s[1]) < (L/4), (0, 0))] = 0 lead[lat.neighbors()] = hopping sys.attach_lead(lead) sys.attach_lead(lead.reversed()) edges = numpy.array([[-6, -2, -4], [2, 6, 4]]) for e in edges: sym = kwant.TranslationalSymmetry([0, numpy.sqrt(3)]) lead = kwant.Builder(sym) def lead_region(pos): x, y = pos return (e[0] < x < e[1]) lead[lat.shape(lead_region, (e[2], 0))] = 0 lead[lat.neighbors()] = hopping sys.attach_lead(lead) sys.attach_lead(lead.reversed()) for x in lead.sites(): x.family.norbs = None return sys.finalized() sys = make_system() kwant.plot(sys, site_lw=0.1, lead_site_lw=0, colorbar=False, show=False) energy = 0.35 reciprocal_phis = numpy.linspace(4, 50, 20) current = numpy.array([-1, 1, 0, 0, 0, 0]) r_long = [] r_hall = [] for phi in 1 / reciprocal_phis: params = {"phi": phi, "salt": ""} smatrix = kwant.smatrix(sys, energy, params=params) cond_mat = smatrix.conductance_matrix() # SOLVES FOR V IN I = G*V voltage = numpy.linalg.solve(cond_mat, current) # DIVIDE BY 2 FOR V = I*R, I = 2 r_long.append(abs(voltage[2]-voltage[4]) / 2) r_hall.append((voltage[2] - voltage[3]) / 2) fig, ax = pyplot.subplots(1, 1) ax.plot(reciprocal_phis, r_long) ax.plot(reciprocal_phis, r_hall) fig.show() --------------------- Link of original thread : https://mail.python.org/archives/list/kwant-discuss@python.org/thread/D3MTMB...