Dear Sir, I am trying to plot the energy as a function of magnetic field for a 3D case, but I am getting tedious results. The system is infinite in two directions and has some width in the third direction. Please have a look at the code attached below. I tried a lot but failed. Is the code correct or I am wrong somewhere. Thank you.
*import kwantimport scipy.sparse.linalg as slaimport matplotlib.pyplot as pltimport tinyarrayimport numpy as npfrom numpy import cos, sin, piimport cmathfrom cmath import expsigma_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]])def make_system(a=1, L=30, W=10, H=10, t=1.0, t_x=1.0, t_y=1.0, t_z=1.0, lamda=0.1, beta=1.05): def onsite(site): return (t_z * cos(beta) + 2 * t) * sigma_z def hoppingx(site0, site1): return (-0.5 * t * sigma_z - 0.5 * 1j * t_x * sigma_x) def hoppingy(site0, site1): return -0.5 * t * sigma_z - 0.5 * 1j * t_y * sigma_y def hoppingz(site0, site1, B): y = site1.pos[1] return (-0.5 * t_z * sigma_z - 0.5 * 1j * lamda * sigma_0) * exp(2 * pi * 1j * B * a * (y-40)) syst = kwant.Builder() lat = kwant.lattice.cubic(a) syst[(lat(z, y, x) for z in range(H) for y in range(W) for x in range(L))] = onsite syst[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz syst[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy syst[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx lead1=kwant.Builder(kwant.TranslationalSymmetry((0,-a,0))) lead1[(lat(z,y,x) for z in range(H)for y in range(W)for x in range(L))]=onsite lead1[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz lead1[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy lead1[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx syst.attach_lead(lead1) syst.attach_lead(lead1.reversed()) lead2=kwant.Builder(kwant.TranslationalSymmetry((-a,0,0))) lead2[(lat(z,y,x) for z in range(H)for y in range(W)for x in range(L))]=onsite lead2[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz lead2[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy lead2[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx syst.attach_lead(lead2) syst.attach_lead(lead2.reversed()) syst = syst.finalized() return systdef analyze_system(): syst = make_system() kwant.plot(syst) Bfields = np.linspace(0, 0.002, 100) energies = [] for B in Bfields: ham_mat = syst.hamiltonian_submatrix(params=dict(B=B), sparse=True) ev, evec = sla.eigsh(ham_mat.tocsc(), k=20, sigma=0) energies.append(ev) #print(energies) plt.figure() plt.plot(Bfields, energies) plt.xlabel("magnetic field [${10^-3 h/e}$]") plt.ylabel("energy [t]") plt.ylim(0, 0.11) plt.showdef main(): syst = make_system() analyze_system()main()* -- With Best Regards NAVEEN YADAV Ph.D Research Scholar Deptt. Of Physics & Astrophysics University Of Delhi.