Dear KWANT developers, I want to plot the current density for a 3D model with finite width in x and y direction. Below is the energy dispersion of the model Hamiltonian. To calculate the current density profile, I have assumed scattering region with a width of 60 unit cells in z-direction and leads with different onsite and same hoppings (as that of scattering region) are attached. For better understanding the code is attached below. For plotting current density, I have used 2D slices of y-z plane. Now my question is that, while building the second system, what onsite and Hopping terms should I have to use? Because I have already build a 3D system with particular onsite and hopping terms. Can I use onsite =0 for the second system? and what should be the hoppings for the second system? can I use simply lat2.neighbors() for the second system? I tried it in different ways, it gives me different results. I don't know which one correct. please suggest me the correct way. [image: image.png] *import kwantfrom numpy import exp,sin,linspace,sqrt, cos, piimport matplotlib as mplfrom matplotlib import pyplot as pltimport tinyarrayimport pandas as pd* *sigma_0 = tinyarray.array([[1, 0], [0, 1]])sigma_x = tinyarray.array([[0, 1], [1, 0]])sigma_y = tinyarray.array([[0, -1j], [1j, 0]])sigma_z = tinyarray.array([[1, 0], [0, -1]])lat=kwant.lattice.cubic(norbs=2)L=10W=80#H=60def onsite(site): return (1.0 * cos(1.05) + 2 * 1.0) * sigma_z* *def onsite_s(site): return (1.0 * cos(1.05) + 2 * 1.0) * sigma_z - 0.75 * sigma_0def hoppingx(site0, site1): return (-0.5 * 1.0 * sigma_z - 0.5 * 1j * 1.0 * sigma_x)def hoppingy(site0, site1): return -0.5 * 1.0 * sigma_z - 0.5 * 1j * 1.0 * sigma_ydef hoppingz(site0, site1): y = site1.pos[1] return (-0.5 * 1.0 * sigma_z - 0.5 * 1j * 0.25 * sigma_0) * exp(-2 * pi * 1j * 0.001 * 1.0 * (y-40))sys=kwant.Builder()sys[(lat(z,y,x) for z in range(0,10) for y in range(W) for x in range(L))]=onsite_ssys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingxsys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingysys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingzsys[(lat(z,y,x) for z in range(10,70) for y in range(W) for x in range(L))]=onsitesys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingxsys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingysys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingzsys[(lat(z,y,x) for z in range(70,80) for y in range(W) for x in range(L))]=onsite_ssys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingxsys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingysys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz* *# The scattering region 0 < W_z < 10 and 60 < W_z > 70 have onsite energy of leads.* *lead=kwant.Builder(kwant.TranslationalSymmetry((-1,0,0)))lead[(lat(0,y,x) for y in range(W) for x in range(L))] = onsite_slead[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingzlead[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingylead[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingxsys.attach_lead(lead)sys.attach_lead(lead.reversed())kwant.plot(sys)sysf=sys.finalized()Sites=list(sysf.sites)#print(Sites)wave=kwant.wave_function(sysf, energy=0.1)psi=wave(0)[0]#df = pd.DataFrame(Sites)#df.to_csv('data2.csv',index=False, header=False)def Plot(x,psi): def cut(site0,site1): return site0.pos[2]==x and site1.pos[2]==x J_0 = kwant.operator.Current(sysf, where=cut) current = J_0(psi) lat2=kwant.lattice.square(norbs=2) sys2=kwant.Builder() sys2[(lat2(site.pos[0],site.pos[1]) for site in Sites if site.pos[2]==x)]=onsite* * #sys2[(lat2(site.pos[0],site.pos[1]) for site in Sites if site.pos[2]==x)]=0* * sys2[kwant.builder.HoppingKind((1, 0), lat2, lat2)] = hoppingz sys2[kwant.builder.HoppingKind((0, 1), lat2, lat2)] = hoppingy* * #sys2[lat2.neighbors()]=-sigma_0* * #lead2=kwant.Builder(kwant.TranslationalSymmetry((-1, 0))) #lead2[(lat2(site.pos[0], site.pos[1]) for site in Sites if site.pos[2]==x)]=onsite_s #lead2[kwant.builder.HoppingKind((1, 0), lat2, lat2)] = hoppingz #lead2[kwant.builder.HoppingKind((0, 1), lat2, lat2)] = hoppingy #sys2.attach_lead(lead2) #sys2.attach_lead(lead2.reversed()) sys2f=sys2.finalized() #kwant.plot(sys2f) kwant.plotter.current(sys2f, current)for i in range(10): Plot(i,psi)* -- With Best Regards Naveen Yadav Ph.D Research Scholar Deptt. Of Physics & Astrophysics University Of Delhi.