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