Dear sir,

As you said I tried the wraparound module but I got errors while diagonalizing the system. I don't know how to use this module for 3D system.

Please help. The code is attached below

```
UserCodeError: Error occurred in user-supplied value function "f".
See the upper part of the above backtrace for more information.
```

sys=kwant.Builder(kwant.TranslationalSymmetry((0,0,-1)))

lat=kwant.lattice.cubic(norbs=2)

sys[(lat(x,y,0) for x in range(L) for y in range(W))]=onsite

sys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx

sys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

sys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz

#lead=kwant.Builder(kwant.TranslationalSymmetry((0,0,-1)))

#lead[(lat(x,y,z) for x in range(L) for y in range(W)for z in range(H))]=onsite

#lead[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx

#lead[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

#lead[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz

#sys.attach_lead(lead)

#sys.attach_lead(lead.reversed())

sys = kwant.wraparound.wraparound(sys, keep=None)

kwant.plot(sys)

sysf=sys.finalized()

ham_mat = sysf.hamiltonian_submatrix()

ev = sla.eigsh(ham_mat, k=31, which='SM')

evecs = ev[1]

prob_dens = np.abs(evecs[:, 30])**2

#print(prob_dens)

#Sites=list(sysf.sites)

#rho = kwant.operator.Density(sysf)

#density = rho(psi)

#wf_sqr = sum(rho(psi) for psi in wf(0))

lat=kwant.lattice.cubic(norbs=2)

sys[(lat(x,y,0) for x in range(L) for y in range(W))]=onsite

sys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx

sys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

sys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz

#lead=kwant.Builder(kwant.TranslationalSymmetry((0,0,-1)))

#lead[(lat(x,y,z) for x in range(L) for y in range(W)for z in range(H))]=onsite

#lead[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingx

#lead[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

#lead[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingz

#sys.attach_lead(lead)

#sys.attach_lead(lead.reversed())

sys = kwant.wraparound.wraparound(sys, keep=None)

kwant.plot(sys)

sysf=sys.finalized()

ham_mat = sysf.hamiltonian_submatrix()

ev = sla.eigsh(ham_mat, k=31, which='SM')

evecs = ev[1]

prob_dens = np.abs(evecs[:, 30])**2

#print(prob_dens)

#Sites=list(sysf.sites)

#rho = kwant.operator.Density(sysf)

#density = rho(psi)

#wf_sqr = sum(rho(psi) for psi in wf(0))

On Fri, May 8, 2020 at 3:30 AM Abbout Adel <abbout.adel@gmail.com> wrote:

Sorry, Naveen. I don't have an example ready to use.The idea behind is that the wraparound module helps you in getting the homiltonian for each K point: H(kx,ky,kz)Diagonalizing it will give you E(kx,ky,kz) which means a point (or few points for a multiband system)the eigenvectors will help you to calculate the density for a given mode.I hope this helps,AdelOn Thu, May 7, 2020 at 9:46 PM Naveen Yadav <naveengunwal72@gmail.com> wrote:Dear sir,Could you please provide me an working example of this type?Best Regards

Naveen Yadav

Research Scholar

Department of Physics & Astrophysics

University of Delhi

New Delhi-110007On Thu, May 7, 2020, 20:37 Abbout Adel <abbout.adel@gmail.com> wrote:Hi again,I f you want translational symmetry, you need to use wraparound module and you will need to find your result only on one unit cell.You will have also to do integration on the Brillouin zone.I hope this helpsAdelOn Thu, May 7, 2020 at 5:12 PM Naveen Yadav <naveengunwal72@gmail.com> wrote:Dear sir,I understand what you have said. But how can I maintain translation symmetry ?because my system is 3D and for plotting current I access the sites of 3D system usinglist(sys.sites)and plot the current for 2D slice. Please suggest me.Best Regards

Naveen Yadav

Research Scholar

Department of Physics & Astrophysics

University of Delhi

New Delhi-110007On Thu, May 7, 2020, 19:06 Abbout Adel <abbout.adel@gmail.com> wrote:Dear Naveen,What you get is what is expected. You do not have translational symmetry.Rewrite your code by keeping only two dimensions and you will see why.I hope this helps,AdelOn Thu, May 7, 2020 at 12:52 PM Naveen Yadav <naveengunwal72@gmail.com> wrote:Dear KWANT Developers,I am trying to plot the current density. The procedure is straightforward. I have attached leads to the scattering region (leads have same onsite and hopping as of the scattering region) assys[(lat(z,y,x) for z in range(H) for y in range(W)for x in range(L))]=onsite

sys[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz

sys[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

sys[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx

lead=kwant.Builder(kwant.TranslationalSymmetry((1,0,0)))

lead[(lat(z,y,x) for z in range(H) for y in range(W)for x in range(L))]=onsite

lead[kwant.builder.HoppingKind((1, 0, 0), lat, lat)] = hoppingz

lead[kwant.builder.HoppingKind((0, 1, 0), lat, lat)] = hoppingy

lead[kwant.builder.HoppingKind((0, 0, 1), lat, lat)] = hoppingx

sys.attach_lead(lead, add_cells=80)

sys.attach_lead(lead.reversed())

As you can see from the plot there is a discontinuity in the plot at (W=20). Why it is so even if the scattering region and leads have same onsite and hopping and there is translation symmetry throughout ?

--Best Regards,Naveen YadavResearch Scholar

Department of Physics & AstrophysicsUniversity Of DelhiNew Delhi-110007--Abbout Adel--Abbout Adel--Abbout Adel

--

Best Regards,

Department of Physics & Astrophysics

Naveen Yadav

Research ScholarDepartment of Physics & Astrophysics

University Of Delhi

New Delhi-110007