Mailman 3 September 2017 - Kwant-discuss

up/down spin conductance from Greens function like smatrix
by SUDIN GANGULY Sept. 19, 2017

Sept. 19, 2017

Dear Sir, In the new version of Kwant, it is possible to calculate the up and down spin conductance. In order to calculate the different components of the spin conductance (G_x,G_y,G_z), Greens function can be used as shown in the following mail-archive link:https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00483.… If I want to calculate the quantities like G_x^up or G_x^down, is it possible to calculate these with Greens function like the smatrix. Another question I want to ask (forgive me if it is out of the context). We can calculate the zero temperature charge shot-noise using the kwant.physics.two_terminal_shotnoise. Is there any option to calculate the spin shot-noise? With Regards, Sudin Ganguly

2 2

How to plot the local dos for both the scattering region and leads
by kuangyia lee Sept. 19, 2017

Sept. 19, 2017

Dear all, I am learning Kwant through tutorial. My toy model is a two-terminal quantum wire system based on the square lattice. I want to plot the spatial local dos for both the scattering region and leads. However, it does not work for leads part. The following is the code that presents my problem. I am wondering whether it is possible to solve this problem with Kwant. Thanks a lot for any help. Kind regards, kuangyia ============================================================= #!/usr/bin/env python3 import kwant import numpy as np import matplotlib.pyplot as plt def make_system(a=1, t=1.0, W=20, L=50, L_well=10): # Start with an empty tight-binding system and a single square lattice. # `a` is the lattice constant (by default set to 1 for simplicity. lat = kwant.lattice.square(a) sys = kwant.Builder() #### Define the scattering region and its potential profiel. #### def potential(site, pot): (x, y) = site.pos if (L - L_well) / 2 < x < (L + L_well) / 2: return pot else: return 0 def onsite(site, pot=0): return 4 * t + potential(site, pot) sys[(lat(x, y) for x in range(L) for y in range(W))] = onsite sys[lat.neighbors()] = -t #### Define and attach the leads. #### lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0))) lead[(lat(0, j) for j in range(W))] = 4 * t lead[lat.neighbors()] = -t sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys def main(): sys = make_system() # Check that the system looks as intended. kwant.plot(sys) # Finalize the system. sys = sys.finalized() # Calculate ldos at a given energy well_depth = 0.55 ldos = kwant.ldos(sys, energy=0.25, args=[-well_depth]) kwant.plotter.map(sys, ldos, num_lead_cells=10) if __name__ == '__main__': main()

3 3

[kwant] lattice constant and the density of states in four-band Hamiltonian
by Gongwei_Hu Sept. 18, 2017

Sept. 18, 2017

Hello everyone! I want to solve a quantum transport question about the quantum spin Hall bar. The Hamiltonian of my system is The tight-binding model can be written as follow ( is the lattice constant): Onsite item: ; Hopping-x item: ; Hopping-y item: The matrix: My Question: (1) The first one is that I am not sure whether the tight-binding model presented above is right for the case of , such as , especially for the hopping -x and -y item. (2) In order to calculate the density distribution of electron and the conductance, the code presented below cannot exactly match the results in the paper [https://journals.aps.org/prb/abstract/10.1103/PhysRevB.83.081402]. For example the density of states calculated by Kwant near the edge is much thick compared with the paper results. The conductance also cannot give the same result, such as the conductance is 1.1e^2/h for my result but 2.0e^2/h for the paper at Fermi energy 10meV. My kwant code is import math import numpy from numpy import * from matplotlib import pyplot import kwant import tinyarray sigma_0 = tinyarray.array([[1, 0, 0, 0], [0, 1, 0 ,0], [0, 0, 1, 0], [0, 0, 0, 1]]) sigma_1 = tinyarray.array([[0, -1, 0, 0], [-1, 0, 0 ,0], [0, 0, 0, 1], [0, 0, 1, 0]]) sigma_2 = tinyarray.array([[1, 0, 0, 0], [0, -1, 0 ,0], [0, 0, 1, 0], [0, 0, 0, -1]]) sigma_3 = tinyarray.array([[0, -1, 0, 0], [1, 0, 0 ,0], [0, 0, 0, -1], [0, 0, 1, 0]]) def make_system(W_QPC): L1=70 W1=60 L2=70 W=150 A=364.5 B=-686 C=0 D=-512 M=-10 a=2 b=1 def central_region(pos): x, y = pos return -(L1 +W1 + L2)/2 < x < (L1 +W1 + L2)/2 and \ abs(y) < (x <= L1-(L1 +W1 + L2)/2) * abs(0.25 * (2 * W-W_QPC) * (math.sin((x+(L1 +W1 + L2)/2)*math.pi/L1 - math.pi/2) + 1)-W) + (x > L1-(L1 +W1 + L2)/2 and x <= (L1+W1-(L1 +W1 + L2)/2)) * W_QPC/2 + (x > (L1 + W1-(L1 +W1 + L2)/2) and x < (L1+W1+L2-(L1 +W1 + L2)/2)) * abs(0.25 * (2 * W-W_QPC) * (math.sin((x-W1/2)*math.pi/L1 + math.pi/2) + 1)-W) lat = kwant.lattice.square(a) sys = kwant.Builder() sys[lat.shape(central_region, (0, 0))] = (C-4*D/a**2)*sigma_0 + (M-4*B/a**2)*sigma_2 sys[kwant.builder.HoppingKind((b, 0), lat, lat)] = D * sigma_0/a**2+B * sigma_2/a**2+1j * A * sigma_1/a sys[kwant.builder.HoppingKind((0, b), lat, lat)] = D * sigma_0/a**2+B * sigma_2/a**2+ A * sigma_3/a sym = kwant.TranslationalSymmetry((-a, 0)) lead = kwant.Builder(sym) lead[(lat(0, y) for y in range(-int(W/2) + 1, int(W/2)))] = (C-4*D/a**2)*sigma_0 + (M-4*B/a**2)*sigma_2 lead[kwant.builder.HoppingKind((b, 0), lat, lat)] = D * sigma_0/a**2+B * sigma_2/a**2+1j * A * sigma_1/a lead[kwant.builder.HoppingKind((0, b), lat, lat)] = D * sigma_0/a**2+B * sigma_2/a**2+ A * sigma_3/a sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys.finalized() def density(sys, energy): psi = kwant.wave_function(sys, energy) for i in range(2): psi_5 = abs(psi(0)[i]**2)#.sum(axis=0) #psi_5 = abs(psi(0)[5]**2) #numpy.savetxt("psi_n.txt",psi_n); A1, A2 = psi_5[::4], psi_5[1::4]#, psi_n[2::4], psi_n[3::4] D = A1+A2 numpy.savetxt("D_now.txt",D); kwant.plotter.map(sys, D) def plot_conductance(WW_QPC, energy): # Compute conductance data = [] for w_QPC in WW_QPC: sys = make_system(W_QPC=w_QPC) #sys = sys.finalized() smatrix = kwant.smatrix(sys, energy) data.append(smatrix.transmission(1, 0)) numpy.savetxt("WW_QPC.txt",WW_QPC) numpy.savetxt("data.txt",data) pyplot.figure() pyplot.plot(WW_QPC, data) pyplot.xlabel("W_QPC [nm]") pyplot.ylabel("conductance [e^2/h]") pyplot.show() sys = make_system(W_QPC=50) kwant.plot(sys) d = density(sys, energy=10) plot_conductance(WW_QPC = [50*i+50 for i in range(4)], energy = 10) #numpy.savetxt("D.txt",d); #kwant.plotter.map(sys, d)

2 1

Re: [Kwant] 2D Current Density Streamplot of 3D System
by mitchell.greenberg Sept. 17, 2017

Sept. 17, 2017

Hi Anton, Thanks for that. While my lattice is polyatomic I just switched site.pos to site.tag in the code I posted last message to get the tags for my lattice. Doing that I see the tags for the z axis range from 0 to 9 (makes sense as I have 10 layers), which explains why 5.3 worked when I used it but not 10. My function runs now an accurately deletes the sites without any problems. Thank You, Mitchell Greenberg Sent from Samsung Mobile <div>-------- Original message --------</div><div>From: Anton Akhmerov <anton.akhmerov+kd(a)gmail.com> </div><div>Date:17/09/2017 5:11 AM (GMT+10:00) </div><div>To: Mitchell Greenberg <mitchell.greenberg(a)my.jcu.edu.au> </div><div>Cc: Joseph Weston <joseph.weston08(a)gmail.com>,kwant-discuss(a)kwant-project.org </div><div>Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System </div><div> </div>Hi Mitchel, Deleting and adding sites is done by their lattice coordinates (site.tag), not their real space position (site.pos); if you want to access a site close to a certain position, you can use the method Monatomic.closest (see https://kwant-project.org/doc/1/reference/generated/kwant.lattice.Monatomic…). Best, Anton On Sat, Sep 16, 2017 at 2:34 PM, Mitchell Greenberg <mitchell.greenberg(a)my.jcu.edu.au> wrote: > Hi Joseph and Anton, > > > Thank you for the advice with getting the 2d graph up and running, its now > producing results that look good. Sorry that I haven't replied in over 2 > weeks, uni has gotten into full swing. > > > I wasn't sure if I needed to start a new topic for this (if I need to let me > know) but I am now trying to remove lattice points randomly from my sample. > I have created a function that can remove random lattice points (it is > rather messy) however I have noticed that when I try to delete sites for > layers other than at z=0 and z=5.3 I get the following error.: > > > KeyError: Site(kwant.lattice.Monatomic([[4.3763, 0.0, 0.0], [0.0, 3.3136, > 0.0], [0.0, 0.0, 5.3]], [0.0, 0.0, 0.0], '0', 1), array([50, 50, 10])) > > > This error seems to pop up when I try to delete a site that doesn't exist at > the specified coordinates however when I use: > > Positions=[syst.sites[i].pos for i in range(len(syst.sites))] > > > I see that sites definitely exist for z=10 (should be 10.6 but the error > output seems to round). I was therefore wondering if there is something I am > missing when it comes to deleting sites and if there is an easy way to > delete random sites from a 3d sample. > > > Thank You, > > Mitchell Greenberg > > > ________________________________ > From: Joseph Weston <joseph.weston08(a)gmail.com> > Sent: Wednesday, August 30, 2017 8:33:01 PM > To: Mitchell Greenberg; kwant-discuss(a)kwant-project.org > Subject: Re: [Kwant] 2D Current Density Streamplot of 3D System > > Hi Mitchell, > >> Sorry again for the reply, I've had a look at the documentation and it >> seems that relwidth and abswidth control the size of the arrays. What >> I'm stuck on though is that my system has L = 80, W = 450 and H =52, >> I was expecting field to return something of the form (80,450,52,3), I >> must apologise again, I'm pretty new to all things python so my >> knowledge in the more complex areas is lacking. > > The input to 'interpolate_current' is a current defined on the hoppings > of your system. The output is a current defined over a realspace grid. > First the current defined on the hoppings is convoluted with a bump > function to smooth it out, then the resulting continuous vector field is > sampled on a regular grid. The docstring explains how this is done. The > interpolation grid does not in general coincide with the "grid" of sites > that form your system. Typically there will be several interpolation > points between neighboring sites. > > The 'box' returned by 'interpolate_current' allows you to calculate the > realspace position of the interpolation points. > > Also the modifications to 'interpolate_current' so that it functions > correctly for 3D systems just landed on the most recent development > version of kwant [1]. > > > Happy Kwanting, > > Joe > > > [1]: > https://gitlab.kwant-project.org/kwant/kwant/commit/be364e7c64c316caf86869d… > >

1 0

2D Current Density Streamplot of 3D System
by Mitchell Greenberg Sept. 16, 2017

Sept. 16, 2017

Hello, I am currently completing my thesis using Kwant, per discussions with my supervisor we would like to plot the current density as a streamplot of our 3d system for different z-axis 'slices' of the system (our system has an exponentially decaying potential in the z direction, so we'd expect current to flow differently through each layer). I have tried tinkering with the new current functions implemented in 1.3 but the plotting function can only plot 2d systems. Is there a way to obtain the current density streamplot for a 2d 'slice' of our 3d system. Thank You, Mitchell Greenberg

3 13

Understanding the scaling of the energy and wavenumber units of measurement
by George Datseris Sept. 13, 2017

Sept. 13, 2017

Hello, I am having significant trouble understanding which are the units of measurement of the energy and the wavenumber when using kwant *and when using different lattice constants*. I am trying to reproduce the results of the (probably well-known paper) "Scalable Tight Binding Model for Graphene", Phys. Rev. Lett. 114, 036601 (2015). I have read the answers to another question (https://www.mail-archive.com/kwant-discuss@kwant-project.org/msg00069.html) but I simply cannot understand which is the actual unit of measurement used internally so I can translate those to nanometers and electronVolts. Even though I have realized how I can translate internal wavenumbers into nanometers^{-1}, I cannot understand how *and why* to use properly kwants.physics.Bands(). Here is the minimal code that captures my problems: ```python import kwant from matplotlib import pyplot a = 0.142 #inter-carbon distance in nm lc = sqrt(3)*a #lattice constant in nm t = 2.8 #hopping in eV W = 200 #width in UNSCALED lattice constants, as if sf=1. (multiply by lc to get it in nm) def bandstructure(sf): #sf = scaling factor glat = kwant.lattice.honeycomb(sf) #sf*la for lattice constant in nm A, B = glat.sublattices # create an armchair lead def leadshape(pos): #set leads width x, y = pos return (0 <= x < W) sym_lead_vertical = kwant.TranslationalSymmetry(glat.vec((-1*sf,2*sf))) armchair = kwant.Builder(sym_lead_vertical) armchair[glat.shape(leadshape, (0, 0))] = 0 #onsite energy is 0 for Dirac armchair[glat.neighbors()] = -1/sf #hopping. use t/sf for value in eV #function that gives the bands (energies) at a given wave vector: bands = kwant.physics.Bands(armchair.finalized()) #The wavenumbers are measured in units of 1/lattice_const. #that is why they only need to go from -pi to pi (Bloch): wavenums = np.linspace(-pi/10 ,pi/10, 41) #I now want to get the energies, measured in eV. How???? energies = [bands(k/sf)*(t/sf) for k in wavenums] #???? doesn't work no matter what I change! pyplot.figure(figsize=(6,8)) # Rescale the wavevectors so that it is measured in nm instead of lattice constants wavenums /= sf*lc #this works correctly # plot dirac dispersion first: pyplot.plot(wavenums, (3/2)*a*t*wavenums, color="black") #this works correctly! # The band structure does not work. pyplot.plot(wavenums, energies) pyplot.xlabel(r"Wavevector [nm$^{-1}$]", family = "serif") pyplot.ylabel(r"Energy [eV]", family = "serif") pyplot.xlim(-0.4, 0.4) pyplot.ylim(-0.2, 0.2) bandstructure(1) ``` From previous answers, I cannot understand what corresponds to the "the distance units that I have chosen". I have chosen here the distance units such that `lc -> 1`. When I am creating the lattice with the command honeycomb(sf), are my distance units still `lc` or are they now `lc*sf` ? Notice that the above `bandstructure()` should give similar results independently of scaling factor (see the paper for details on why). Therefore, the only problem that may remain is the scaling of energy and/or wavenumber. Let me summarize my questions as clearly as possible, to help you giving me an easy answer: 1. Which is the unit of measurement of length? I always thought that it is simply `1`! If I write `honeycomb(32.123)` is the unit of measurement of length still `1` or does it become 32.123, irrespectively of how many nanometers 1 actually corresponds to? When I get the lattice sites, I know that they are assumed (i,j) TIMES the lattice constant, so the actual unit of distance is still (i,j)*32.123*1, and `1` corresponds to whatever units I choose. 2. Now, is the wavenumber represented in units of `1`, or in units of 1/32.123? If I have a value of `k=3` and my `1` corresponds to 1 nanometer, does this value of `k` correspond to `3`nm^{-1} or to `3/32.123` nm^{-1} ? 3. What is the input given to the function `bands(k)`, obtained from `kwant.physics.Bands` ? Does it have to *always be within the interval -pi/2 to pi/2 *, meaning that it is always expected in units of 1/lattice const. /irrespectively of which are the units of measurement of the wavenumber? / 4. What is the unit of energy, that is returned by the function `bands()`? Does it depend on the value that I give to the hopping, or it depends on how I define `1` ? E.g. when I write `armchair[glat.neighbors()] = -1/sf` is the unit of energy internally defined to be `1` or `1/sf` ??? It cannot possibly be the latter, because I could also assign some number to the on-site potential. So it must be the first, right? But then, how do you answer my question number 4: // 5. Given that I know the length of the lattice constant in nanometers, and I know the energy of the hopping value in eV, how can I get the bands in eV versus nm^{-1} ? I would also like to comment that the documentation string of `bands()` would solve all of these questions by only having 2 more sentences in it. The sentences that answer questions 2-3-4 would probably be enough. I am sure that the answers to my questions are almost trivial, but after many hours of trying to understand what is going on, I am not so confused that I had to ask. So, sorry in advance for the very easy question! Best, George Datseris MPI For Dynamics & Self-Organization

2 5

A weird result of 1D quantum wire
by Liu Feng Sept. 6, 2017

Sept. 6, 2017

Dear Adel Many thanks for your informative reply. Now I think I have resolved the problem here. The result kwant obtained is correct. One thing is that for odd L, the transmission is always unity, not extinction. This quite counter-intuitive. As even the hopping is given as for example J_0(2.4048) (zeros of Bessel function), it still gives unity transmission. For analytic details please refer to PRB 78, 045428. Best regards, feng Abbout Adel <abbout.adel(a)gmail.com>于2017年9月6日周三上午6:29写道： > Dear Liu, > > I think that there is no problem in your result. It is the expected > physical result! > > This behaviour is due to the fact that you are near the middle of the > conduction band. > > To see it better, put the energy exactly =0. you will find that the > conductance is zero for odd numbers and oscillating for even. > > In fact, for Ef=0, you have the wavenumber k=pi/2 and Lambda=2pi/k=4 > (E=-2t cos(k), whatever is the value of t, when E=0, k=pi/2 ) > > The system is like a tunneling barrierand therefore we will have maximum > transmission if L= N *lambda/2 (like for a quantum point contact) with N > Integer and expect extinction if N is half integer. > > example: L=30, lambda =4 ===> N=15 oscillating at E=0. L=31 ====>N=15+1/2 > extinction. > > Since you are using an energy E=0.0001 a little bit different from 0 so > you do not see the full extinction. This is due to the fact that wavelength > lambda is not integer anymore > > If you use an energy far from the middle of the band, even and odd systems > will behave the same way. > > take more values of das=np.linspace(0,10,201) to see correctly the > oscillating behaviour of the conductance. > > I hope this helps > Adel > > On Tuesday, September 5, 2017, Liu Feng <ruserzzz(a)gmail.com> wrote: > >> Dear All, >> >> I have met a strange problem in computing conductance of 1D quantum >> wire. >> The setup is as follows: >> The hopping of the quantum wire is given by the first kind of >> Bessel function, for example, J0(A). >> The coupling of the lead is set as 1. >> As we know in this case the transmission should be proportional to >> J0(A)^2. However, I get very different results for wire of even and odd >> sites numbers as shown in the attachment... >> Is this a bug or do I miss something else? >> >> source code >> ------------------------------------------ >> """ >> @author: liu >> """ >> #import the plot function >> from matplotlib import pyplot >> import numpy as np >> #turn on kwant >> import kwant >> import scipy.special as sp #Bessel function >> >> def family_color(site): >> if site.tag[0] in defect: >> return 'magenta' >> else: >> return 'black' >> >> def conduct(mu,da): >> >> >> sys=kwant.Builder() >> a=1 >> lat=kwant.lattice.square(a) >> t=sp.jn(0,da) >> W=1 >> L=31 >> #da=1.0 >> # mu=0.0 >> for i in range(L): >> for j in range(W): >> # sys[lat(i,j)]=4*t >> if i in defect: >> sys[lat(i,j)]=mu+da >> else: >> sys[lat(i,j)]=mu >> >> if j>0: >> sys[lat(i, j), lat(i, j-1)]=-t >> if i>0: >> sys[lat(i, j), lat(i-1, j)]=-t >> >> #build the lead >> sys_left_lead=kwant.TranslationalSymmetry((-a,0)) >> left_lead=kwant.Builder(sys_left_lead) >> >> #define the hopping in the lead, it is enough to only add one unit cell >> hopping >> #translational symmetry takes care the remainings. >> >> for j in range(W): >> left_lead[lat(-1, j)]=mu >> if j>0: >> left_lead[lat(0, j), lat(0, j-1)]=-1.0 >> left_lead[lat(-1, j), lat(-2, j)]=-1.0 >> >> #attach the left lead >> sys.attach_lead(left_lead) >> >> #similarly we define the right lead >> sym_right_lead=kwant.TranslationalSymmetry((a,0)) >> right_lead=kwant.Builder(sym_right_lead) >> for j in range(W): >> right_lead[lat(L,j)]=mu >> if j>0: >> right_lead[lat(L,j),lat(L,j+1)]=-1.0; >> right_lead[lat(L,j),lat(L+1,j)]=-1.0; >> sys.attach_lead(right_lead) >> >> #plot system >> # kwant.plot(sys,site_color=family_color) >> # kwant.plotter.bands(left_lead.finalized()) >> >> #finalize the system for conductivity computation >> sys=sys.finalized() >> return energies,kwant.smatrix(sys,0.0001).transmission(1,0) >> >> das=np.linspace(0,10,51) >> >> ts=[] >> defect=[1000] >> #energies,data,t=conduct(mu=0.0,da=0.0) >> #pyplot.plot(energies,data) >> for tda in das: >> energies,t=conduct(mu=0,da=tda) >> ts.append(t) >> >> #plot the conductance >> >> >> >> >> >> #pyplot.ylim([0,1]) >> pyplot.show() >> pyplot.figure() >> pyplot.plot(das,ts) >> pyplot.xlabel('A') >> pyplot.ylabel('conductance[e^2/h]') >> >> >> Best regards, >> Feng Liu >> >> > > -- > Abbout Adel >

1 0

A weird result of 1D quantum wire
by Liu Feng Sept. 6, 2017

Sept. 6, 2017

Dear Joe, Many thanks for your reply. After replacing all the manually given hoppings by lat.neighbors() method, I still obtain the same result. Also, the system is plotted for the length L of even and odd numbers, which gives corrected tight-binding system. So I think as Adel said in his reply, the result is the expected physical one. After checking some materials, I found it is a known effect in clean 1D quantum wires which is called even-odd oscillation, for example, discussed in PRB 78, 045428 and PRB 41, 8553. Best regards, feng Joseph Weston <joseph.weston08(a)gmail.com>于2017年9月6日周三上午1:08写道： > Dear Liu, > > > > I have met a strange problem in computing conductance of 1D > > quantum wire. > > The setup is as follows: > > The hopping of the quantum wire is given by the first kind of > > Bessel function, for example, J0(A). > > The coupling of the lead is set as 1. > > As we know in this case the transmission should be proportional > > to J0(A)^2. However, I get very different results for wire of even and > > odd sites numbers as shown in the attachment... > > Is this a bug or do I miss something else? > > > > This is probably not a bug in Kwant. I see in your script that you do > not construct your system in a very idiomic manner (you manually add > nearest-neighbor hoppings etc). Are you sure you did not make a mistake > during system construction? Did you plot the system in the two cases > (even/odd) and see that everything is as you expect? > > Happy Kwanting, > > Joe > >

1 0

the density of states in four-band Hamiltonian
by Gongwei_Hu Sept. 6, 2017

Sept. 6, 2017

Dear Joseph ! Thanks for your teams creating the Kwant to calculate the quantum transport, and I have read some of the kwant-discussions that is relevant with my questions. However, there are still several confuse and I hope that you can help me for your free time. (1) I am not sure whether the energy used in function (such as kwant.smatrix(sys, energy), kwant.wave_function(sys, energy) or kwant.ldos(sys, energy).) is Fermi energy in a given system. My previous insight about this energy is the energy that one electron is incoming from one lead. (2) In some literatures such as topological insulator, the Hamiltonian has the form blocks(four-band Hamiltonian). For example, in ref [Papaj M, Cywiński Ł, Wróbel J, et al. Physical Review B, 2016, 93(19): 195305.], the Hamilton takes the following form , . If we use the function ( local_ldos =)kwant.ldos(sys, energy) to calculate local the density of states (LDOS) for this system , the return array local_ldos is . So, my question is that whether we can think the following Schrodinger equation holds for site- . Other relevant questions: (a) Dose the eigenvalue correspond to the energy inputted into kwant.ldos(sys, energy)? (b) For site-, whether the LDOS of electron with spin up and down is , respectively? If not, how can we obtain it?. (3) When we use the function ( wf =)kwant.wave_function(sys,enengy) to obtain the wave function of electron, the array wf(0)(with elements) for electron incoming from lead 1 can give M models (channels) and oribits (each orbit has four wave function) For the with model , the Schrodinger equation is The question: (a) Whether the eigenvalue is all the same for . In other words, I am not sure that whether all the states for is degenerate in ? (b) In order to calculate the density of states for , we can use the following relation Does this give the total density of states? If we only want to know the case of spin up or down, how can we calculate the density of states? Is the following correct? (c) What is the difference of the density of state obtained from above method (kwant.wave_function)and LDOS(kwant.ldos)? The code for calculating the density of states is the follow: import math import numpy as np import numpy from numpy import * from matplotlib import pyplot import kwant import tinyarray def make_system(W = 197, L = 400, F = 190, de = 100): A, B, D, M, Dz, De, Dh, Ee = 0.365, -0.686, -0.511, -0.01, 0.0016, -0.0128, 0.0211, 0.0 p = 1 sigma_1 = tinyarray.array([[1, 0, 0, 0], [0, -1, 0 ,0], [0, 0, 1, 0], [0, 0, 0, -1]]) sigma_2 = tinyarray.array([[1, 0, 0, 0], [0, 0, 0 ,0], [0, 0, 1, 0], [0, 0, 0, 0]]) sigma_3 = tinyarray.array([[0, 0, 0, 0], [0, 1, 0 ,0], [0, 0, 0, 0], [0, 0, 0, 1]]) sigma_33 = tinyarray.array([[0, 0, 0, -1], [0, 0, 1 ,0], [0, 1, 0, 0], [-1, 0, 0, 1]]) sigma_4 = tinyarray.array([[0, A, De+1j*Ee, 0], [A, 0, 0 ,Dh], [De+1j*Ee, 0, 0, -A], [0, Dh, -A, 0]]) sigma_5 = tinyarray.array([[0, 1j*A, 1j*De+Ee, 0], [-1j*A, 0, 0 ,-1j*Dh], [-1j*De-Ee, 0, 0, 1j*A], [0, 1j*Dh, -1j*A, 0]]) def central_region(pos): x, y = pos return -L/2 < x < L/2 and \ abs(y) < L/2 - F*math.exp((-x**2)/(de**2)) lat = kwant.lattice.square(p) sys = kwant.Builder() sys[lat.shape(central_region, (0, 0))] = M*sigma_1-4*(B+D)*sigma_2+4*(B-D)*sigma_3+Dz*sigma_33 sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = (B+D)*sigma_2-(B-D)*sigma_3-1j*sigma_4 sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = (B+D)*sigma_2-(B-D)*sigma_3-1j*sigma_5 sym = kwant.TranslationalSymmetry((-p, 0)) lead = kwant.Builder(sym) lead[(lat(0, y) for y in range(-W + 1, W))] = M*sigma_1-4*(B+D)*sigma_2+4*(B-D)*sigma_3+Dz*sigma_33 lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = (B+D)*sigma_2-(B-D)*sigma_3-1j*sigma_4 lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = (B+D)*sigma_2-(B-D)*sigma_3-1j*sigma_5 sys.attach_lead(lead) sys.attach_lead(lead.reversed()) return sys.finalized() def plot_ldos(sys, energy): local_ldos = kwant.ldos(sys, energy) sx1 = local_ldos[::4] sx2 = local_ldos[1::4] sx3 = local_ldos[2::4] sx4 = local_ldos[3::4] kwant.plotter.map(sys, sx1+sx2) def density(sys, energy): wf = kwant.wave_function(sys, energy) wf_0 = abs(wf(0)**2).sum(axis=0) A1, A2, A3, A4 = wf_0[::4], wf_0[1::4], wf_0[2::4], wf_0[3::4] D = A1 + A2# + A3 + A4 #numpy.savetxt("D_D.txt",D); kwant.plotter.map(sys, D) sys = make_system() kwant.plot(sys) # density of state for spin up d = density(sys, energy=-0.02) # local density of state for spin up plot_ldos(sys, energy=-0.02)

1 0

A weird result of 1D quantum wire
by Liu Feng Sept. 5, 2017

Sept. 5, 2017

Dear All, I have met a strange problem in computing conductance of 1D quantum wire. The setup is as follows: The hopping of the quantum wire is given by the first kind of Bessel function, for example, J0(A). The coupling of the lead is set as 1. As we know in this case the transmission should be proportional to J0(A)^2. However, I get very different results for wire of even and odd sites numbers as shown in the attachment... Is this a bug or do I miss something else? source code ------------------------------------------ """ @author: liu """ #import the plot function from matplotlib import pyplot import numpy as np #turn on kwant import kwant import scipy.special as sp #Bessel function def family_color(site): if site.tag[0] in defect: return 'magenta' else: return 'black' def conduct(mu,da): sys=kwant.Builder() a=1 lat=kwant.lattice.square(a) t=sp.jn(0,da) W=1 L=31 #da=1.0 # mu=0.0 for i in range(L): for j in range(W): # sys[lat(i,j)]=4*t if i in defect: sys[lat(i,j)]=mu+da else: sys[lat(i,j)]=mu if j>0: sys[lat(i, j), lat(i, j-1)]=-t if i>0: sys[lat(i, j), lat(i-1, j)]=-t #build the lead sys_left_lead=kwant.TranslationalSymmetry((-a,0)) left_lead=kwant.Builder(sys_left_lead) #define the hopping in the lead, it is enough to only add one unit cell hopping #translational symmetry takes care the remainings. for j in range(W): left_lead[lat(-1, j)]=mu if j>0: left_lead[lat(0, j), lat(0, j-1)]=-1.0 left_lead[lat(-1, j), lat(-2, j)]=-1.0 #attach the left lead sys.attach_lead(left_lead) #similarly we define the right lead sym_right_lead=kwant.TranslationalSymmetry((a,0)) right_lead=kwant.Builder(sym_right_lead) for j in range(W): right_lead[lat(L,j)]=mu if j>0: right_lead[lat(L,j),lat(L,j+1)]=-1.0; right_lead[lat(L,j),lat(L+1,j)]=-1.0; sys.attach_lead(right_lead) #plot system # kwant.plot(sys,site_color=family_color) # kwant.plotter.bands(left_lead.finalized()) #finalize the system for conductivity computation sys=sys.finalized() return energies,kwant.smatrix(sys,0.0001).transmission(1,0) das=np.linspace(0,10,51) ts=[] defect=[1000] #energies,data,t=conduct(mu=0.0,da=0.0) #pyplot.plot(energies,data) for tda in das: energies,t=conduct(mu=0,da=tda) ts.append(t) #plot the conductance #pyplot.ylim([0,1]) pyplot.show() pyplot.figure() pyplot.plot(das,ts) pyplot.xlabel('A') pyplot.ylabel('conductance[e^2/h]') Best regards, Feng Liu

3 2