Dear kwant developers, I am now working on a 3D cuboid system which is similar to the previous email. I would like to calculate the Rxx and Rxy of my system in the presence of a magnetic field. However, a strange periodicity of Rxy appears and thus I am writing to seek your help to check if I have made hoppings correctly. Becuase of the translational invariance of the hopping for the leads, I used two guages, namely A = (By,0,0) and A = (0,-Bx,0). At the boundary, I interpolate the two guages by using hopping3 and 4. Am I using the right way? Thanks. Best, TC #---------------------------------Here is my code: ------------------------------# import numpy as np import matplotlib import kwant from matplotlib import pyplot as plt from matplotlib import pyplot space = 5 L = 20 W = 20 H = 20 t = 1.0 energy = 0.5 Bs = [4*i for i in xrange(10)] def cube(L, W, H,space): """A cuboid region of BHZ material with 3 leads attached. """ lat = kwant.lattice.general(np.identity(3)) sys = kwant.Builder() sym_xy = kwant.TranslationalSymmetry((0, 0, -1)) sym_yz = kwant.TranslationalSymmetry((-1, 0, 0)) sym_zx = kwant.TranslationalSymmetry((0, -1, 0)) lead_xy = kwant.Builder(sym_xy) lead_yz = kwant.Builder(sym_yz) lead_zx = kwant.Builder(sym_zx) def shape_lead_yz(pos): (x, y, z) = pos return (space <= z < H-space) and (space <= y < W-space) def shape_lead_xy(pos): (x, y, z) = pos return (space <= x < L-space) and (space <= y < W-space) def shape_lead_zx(pos): (x, y, z) = pos return (space <= x < L-space) and (space <= z < H-space) def shape(pos): (x, y, z) = pos return (0 <= z < H) and (0 <= y < W) and (0 <= x < L) def hopping1(site1, site2, phi = 0): xt, yt, zt = site1.pos xs, ys, zs = site2.pos return -t * np.exp(-1j * phi * yt) #return -t*np.exp(-0.5j*phi*(xt-xs)*(yt+ys)) def hopping2(site1, site2, phi = 0): xt, yt, zt = site1.pos xs, ys, zs = site2.pos return -t * np.exp(-1j * phi * (-xt)) #return -t*np.exp(-0.5j*phi*(xt+xs)*(yt-ys)) # --------------- interpolate between different guages------------------# def hopping3(site1, site2 , phi = 0): xt, yt, zt = site1.pos xs, ys, zs = site2.pos if (2< yt <= 3 or W-3 > yt >= (W-4)): return -t*(np.exp(-1j * phi *(0.8*yt))) elif (1< yt <= 2 or W-2 > yt >=(W-3)): return -t*(np.exp(-1j * phi * (0.6*yt))) elif (0< yt <= 1 or W-1 > yt >= (W-2)): return -t*(np.exp(-1j * phi *(0.4*yt))) elif ( -1< yt <= 0 or W> yt >= (W-1)): return -t*(np.exp(-1j * phi *(0.2*yt))) elif ( yt <= -1 or yt >= W): return -t else: return -t * np.exp(-1j * phi * (yt)) def hopping4(site1, site2 , phi = 0): xt, yt, zt = site1.pos xs, ys, zs = site2.pos if (2< yt <= 3 or W-3 > yt >= (W-4)): return -t*(np.exp(-1j * phi *(0.2*(-xt)))) elif (1< yt <= 2 or W-2 > yt >=(W-3)): return -t*(np.exp(-1j * phi * (0.4*(-xt)))) elif (0< yt <= 1 or W-1 > yt >= (W-2)): return -t*(np.exp(-1j * phi *(0.6*(-xt)))) elif ( -1< yt <= 0 or W> yt >= (W-1)): return -t*(np.exp(-1j * phi *(0.8*(-xt)))) elif ( yt <= -1 or yt >= W): return -t * np.exp(-1j * phi * (-xt)) else: return -t # scattering system sys[lat.shape(shape, (0,0,0))] = 4*t sys[(lat(-1,y,z) for y in range(space,L-space) for z in range(space,L-space))] = 4*t sys[(lat(L,y,z) for y in range(space,L-space) for z in range(space,L-space))] = 4*t sys[(lat(x,-1,z) for x in range(space,L-space) for z in range(space,L-space))] = 4*t sys[(lat(x,L,z) for x in range(space,L-space) for z in range(space,L-space))] = 4*t sys[(lat(x,y,-1) for x in range(space,L-space) for y in range(space,L-space))] = 4*t sys[(lat(x,y,L) for x in range(space,L-space) for y in range(space,L-space))] = 4*t #sys[lat.neighbors()] = -t sys[kwant.HoppingKind((1,0,0), lat)] = hopping3 sys[kwant.HoppingKind((0,1,0), lat)] = hopping4 sys[kwant.HoppingKind((0,0,1), lat)] = -t # leads no. 0 and 1 # leads along z axis lead_xy[lat.shape(shape_lead_xy, (space,space,space))] = 4*t lead_xy[kwant.HoppingKind((1,0,0), lat)] = hopping1 lead_xy[kwant.HoppingKind((0,1,0), lat)] = -t lead_xy[kwant .HoppingKind((0,0,1), lat)] = -t # leads no. 2 and 3 # leads along x axis lead_yz[lat.shape(shape_lead_yz, (space,space,space))] = 4*t lead_yz[kwant.HoppingKind((1,0,0), lat)] = hopping1 lead_yz[kwant.HoppingKind((0,1,0), lat)] = -t lead_yz[kwant.HoppingKind((0,0,1), lat)] = -t # leads no. 4 and 5 # leads along y axis lead_zx[lat.shape(shape_lead_zx, (space,space,space))] = 4*t lead_zx[kwant.HoppingKind((1,0,0), lat)] = -t lead_zx[kwant.HoppingKind((0,1,0), lat)] = hopping2 lead_zx[kwant.HoppingKind((0,0,1), lat)] = -t # attaching leads (this line determines the lead #) sys.attach_lead(lead_xy) sys.attach_lead(lead_xy.reversed()) sys.attach_lead(lead_yz) sys.attach_lead(lead_yz.reversed()) sys.attach_lead(lead_zx) sys.attach_lead(lead_zx.reversed()) return sys.finalized() def plot_conductance(sys, energy, Bs): #Bs is the magnetic field n = 6 #for temp in xrange(2): data1 = [] data2 = [] data3 = [] for B in Bs: smatrix = kwant.smatrix(sys, energy, args=[B]) s = smatrix # First we calculate the conductance matrix given the scattering matrix S. cond = np.array([[s.transmission(i, j) for j in xrange(n)] for i in xrange(n)]) # n is the number of leads required cond -= np.diag(cond.sum(axis=0)) print (cond) cm = cond[:-1,:-1] # cm stands for conductance matrix nonlocal_resistance0 = np.linalg.solve(cm, [0, 0, 1,-1,0])[0] nonlocal_resistance1 = np.linalg.solve(cm, [0, 0, 1,-1,0])[1] nonlocal_resistance2 = np.linalg.solve(cm, [0, 0, 1,-1,0])[2] nonlocal_resistance3 = np.linalg.solve(cm, [0, 0, 1,-1,0])[3] nonlocal_resistance4 = np.linalg.solve(cm, [0, 0, 1,-1,0])[4] data1.append(-nonlocal_resistance4) #rxy data2.append(nonlocal_resistance3-nonlocal_resistance2) #rxx data3.append(nonlocal_resistance1-nonlocal_resistance0) #rxz pyplot.figure() rxy = pyplot.plot(Bs, data1) rxx = pyplot.plot(Bs, data2) rxz = pyplot.plot(Bs, data3) pyplot.xlabel("magnetic field [T]") pyplot.ylabel("R") plt.setp(rxy,'color','r') plt.setp(rxx,'color','g') plt.setp(rxz,'color','b') plt.savefig('try_R_'+'L,W,H='+repr(L)+','+repr(H)+','+repr(H)+'_space='+repr(space)+'_E='+repr(energy)+'_2.png') #pyplot.show() def main(): sys = cube(L, W, H,space) kwant.plot(sys) plot_conductance(sys,energy,Bs) if __name__ == '__main__': main()