Dear Adel and Krzysiek, 1. thank you for taking time to look into my code and answering my questions in such detail. Sorry also for the late reply. Finally I did not use wraparound. As suggested by Krzysiek, starting from a 3D Hamiltonian in k-space I rewrote the code into a 2D real space model, keeping k_y as a parameter in the leads (and everywhere). The sytem is finite in z direction and x is the transport direction. Then I let kwant obtain the transmission for each k_y and integrated over k_y to obtain the total transmission T_tot(E)=\int dk/(2\pi)T(k_y,E). I hope this is ok.
I am not so sure why I should use wraparound in my situation when I have such a “simple” translational invariance. In the discussion archive I found examples where someone had repeating structures perpendicular to the transport direction, for example like periodical slits/holes in a potential barrier. Probably this is what wraparound is made for. Or maybe also to avoid rebuilding the system for each ky? 2. Just as a remark on a little problem I experienced and would like to share. As I wanted to create bigger slabs I parallelised the loop I used to go through the k_y values. However the code turned out to not benefit at all from this, due to some, as far as I understand, intrinsic parallelisation in the libraries that kwant uses. For any one who is interested, this can be shut down by os.environ["OMP_NUM_THREADS"] = "1" then each kwant job runs on one core only and the speed increases linearly with the amount of cores again :). 3. One more quick question please. I would like to plot something like E(ky,kz=0) for my system, and would like to compare it to the exact diagonalisation I did independently of kwant for the uncoupled leads/the bulk on each side of the interface. I found the snippet https://kwant-project.org/doc/1/reference/generated/kwant.physics.Bands#kwant.physics.Bands >>> bands = kwant.physics.Bands(some_syst) >>> momenta = numpy.linspace(-numpy.pi, numpy.pi, 101) >>> energies = [bands(k) for k in momenta] >>> pyplot.plot(momenta, energies) >>> pyplot.show() so for my case I can use this piece of code for each ky seperately setting momenta to kz=array([0]) or are there other ways that I overlooked when going through the kwant discussion page/documentation? Thank you again Adel and dziękuję Krzysiek! With best regards, Bogusz > On 8 Jun 2019, at 20:51, Abbout Adel <[email protected]> wrote: > > Dear Bogusz, > > I checked your program and it seems fine and the result should be 'almost' > correct. > However, I would like to make some comments and corrections: > > 1) You are using wraparound and your system if uniform in the y-direction. No > need to use a width W. A one unit cell length is enough. > This may seem at a first glance with no consequences, but it does change the > definition of your ky such that two people using two different widths > will have two different values of T(ky, E) although the integrated value > T(E) should be the same (if integrated on the correct interval in each case). > > 2) Let us first be interested on the sum of T(ky, E) over some chosen > continuous modes ky and discretized open modes in z-direction. Let us also > put the same two leads at the same potential and put the L=1. This will help > us in doing analytic calculation. > The energy of each mode is: > Em=6*t-2*t*cos(ky)-2*t*cos(m*pi/(H+1)) > > There is no inter-mode scattering in your system and the transmission can be > obtained using the Fisher-Lee formula as for a scatterer in 1D: > > T(ky,E)=sum_m 1/(1.+(V0/Gamma)**2) with Gamma=sqrt(4-(E-Em)**2). > > This formula is tested in the below code. > > If you change W>1, it will not work for the reason I mentioned before: some > re normalization of the ky (ky/W) and changes should be done probably. I > tried to figure out how to do it but I failed. > > If you put V0 in one lead (as in your program) a quick calculation can give > you the analytical result in the same way as here (you will need Gamma1, > Gamma2) > > I hope this helps > Adel > > > ############################ > import numpy as np > import matplotlib.pyplot as plt > import scipy as sp > from numpy import * > import kwant > import sys > > #calculate the conductance for a sheet (Dim=2) or a Slab (Dim=3) that is > translational invariant in y direction > > a = 1.0 #lattice constant > V0 = 0.2 #potential step size > > L = 1 #Length x-direction the length is a bit useless in > this code > W = 1 #Width y-direction > H = 4 #Hight z-direction > t = 1 #hopping parameter > > Dim=3 #dimension of the system > # for H=1, Dim=2 > # for H>1, Dim=3 > #params for integration and plotting > Nky=100 #points in ky space > NE=30 #points in energy > lat = kwant.lattice.cubic(a) > > > # Define the scattering region > # Infinite potential plane/Slab in y direction > sys = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((0, W, 0)))) > > > for i in range(L): > for j in range(W): > for k in range(H): > # On-site Hamiltonian > sys[lat(i, j, k)] = 2*Dim * t + V0 #potential > step located at the center between the leads > > sys[lat.neighbors(1)] = -t > > sys = kwant.wraparound.wraparound(sys) > > #Left lead > leadL = kwant.Builder(kwant.TranslationalSymmetry((-a, 0, 0), lat.vec((0, W, > 0)))) #trans inv in transport dir and ky > leadL[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t > #onsite terms in the left lead > leadL[lat.neighbors(1)] = -t > > #Right lead > leadR = kwant.Builder(kwant.TranslationalSymmetry(( a, 0, 0), lat.vec((0, W, > 0)))) > leadR[(lat(0,j,k) for j in range(W) for k in range(H))] = 2*Dim*t > #onsite terms in the right lead > leadR[lat.neighbors(1)] = -t > > leadL = kwant.wraparound.wraparound(leadL, keep=0) #keeping the > translational invariance in x direction > leadR = kwant.wraparound.wraparound(leadR, keep=0) > > sys.attach_lead(leadL) > sys.attach_lead(leadR) > > kwant.plot(sys) > sys = sys.finalized() > > # Calculation of the total transmission for a given Energy > ky_array = np.linspace(-pi,pi,Nky) > dky=ky_array[1]-ky_array[0] > energies = np.linspace(0.0, 5.0,NE) > Ttot=[] > > for energy in energies: > > Tky_array=zeros(Nky) > > for i in arange(len(ky_array)): > ky=ky_array[i] > smatrix = kwant.smatrix(sys, energy, [ky]) > Tky_array[i]=smatrix.transmission(1, 0) > > #Ttot.append((2*pi)/(W)*(sum(Tky_array)-0.5*Tky_array[0]-0.5*Tky_array[Nky-1])*dky) > #cheap trapezodial rule to integrate over ky > Ttot.append(sum(Tky_array)) #cheap trapezodial rule to integrate over ky > > > # Plot transmission > plt.plot(energies, Ttot) > plt.show() > > def T(E,qy,W,H,V0,t=1): > T0=0 > for m in range(1,H+1): > Em=6*t-2*t*cos(qy/W)-2*t*cos(m*pi/(H+1)) > > if abs(E-Em)<=2: # condition for open modes > Gamma=sqrt(4-(E-Em)**2) > T0+=1/(1.+(V0/Gamma)**2) ##use Fisher Lee Formula T=Gamma G > Gamma G^\dagger > return T0 > > > > Trans=[sum(np.array([T(E,qy,W,H,V0,t=1) for qy in ky_array ] )) for E in > energies] > > > plt.plot(energies,Trans,'ro') > plt.plot(energies, Ttot) > plt.show() > > On Thu, Jun 6, 2019 at 7:03 PM Christoph Groth <[email protected] > <mailto:[email protected]>> wrote: > Bogusz Bujnowski wrote: > > > I would like to solve the 3D Schrödinger equation with a 1D potential > > step. With the material in the mailing list as well as from the > > tutorial the 2D case works fine. > > > > Extending to 3D gives the same when just taking one lattice site in > > the z direction (H>1). Increasing the hight of the slab (H>1) changes > > the minimal energy for which transmission is possible thus my code is > > wrong. > > > > Did I define and attach the leads correctly with the wraparound in 3D? > > The short code is below. > > From a quick look at your script, the way you are using wraparound seems > fine. Unfortunately, I don’t have the time to debug your problem, but > here’s a suggestion of what you can do: > > You are using a 3D model with two translational symmetries. One of them > you convert to a momentum parameter using wraparound, and the other one > you use for the leads. > > You could simplify your script to use a 2D model, but still use two > translational symmetries in the same way as now. This way, you can more > easily examine your system by plotting it, and the calculations are > quicker. > > When plotting your wrapped-around system in 2D, you should see hoppings > that "wrap around" (hence the name). This should be equivalent to > creating such a Builder by hand (and it should be quite easy to do in > your case). So you should be able to debug the problem. > > Please let us know when you find a solution. We are particularly > interested if you find some problem with Kwant (or its documentation). > > Owocnego Kwantowania! > > Krzysiek > > > -- > Abbout Adel
