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