Dear T.C. Wu, Your problem isn't related to how you define the system, but rather with the units you use for the fields. You put more than one flux quantum of magnetic field per unit cell of your lattice, while all the observable properties have a periodicity of 2 pi in B. Try something like Bs = np.linspace(0, 0.1) Also I recommend you to start from a 2D system, since it's much easier to debug. You can also check out Hall bar simulations over here: http://nbviewer.ipython.org/github/topocm/topocm_content/blob/master/w3_pump... (press the 'show input' button to see the source code). Best, Anton On Sat, Jul 11, 2015 at 6:56 AM, T.C. Wu <tcwu@connect.ust.hk> wrote:
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()