From physical reasoning, I would expect that the current through each link in the direction of the current is the same. So I tried to apply the formula (11) from the NJP-introduction to KWANT. But I get a very different looking answers, comparing to scattering matrix approach, as well as current throught different links in the direction of
Thanks a lot for you package which is rather intuitive in using! Could you please help with understanding of the wave function which KWANT gives? I could not understand what kwant.wave_function gives as an output: -- how to distinguish between the modes which are propagating in the leads and which are decaying in the leads? kwant.wave_function gives an answer at any energy. (1)For example, how do get the density in some transport setup? Is this example correct? def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f def density(sys,energy,lead_nr): wf = kwant.wave_function(sys,energy) return (abs(wf(lead_nr))**2).sum(axis=0) numleads = 2 Tmu = numpy.array([[Tleft,muleft],[Tright,muright]]) #assuming that this range of energies encompasses well the chemical potential range de = 0.05 nenergies = 100 energies=[de * (i - nenergies/2) for i in xrange(nenergies)] for energy in energies: for num_lead in xrange(numleads): d = d + Fermi(energy, Tmu[num_lead]) * density(sys, energy, num_lead) d = d * de (2) Another example which I do not understand properly is about transport through a translationally invariant 2d bar. the current is different. Why are these two answers so different and one is not oscillating, and another is so much depends on energy? (In 1d the comparission works well) def Fermi(energy, T, mu): f = 1./(1. + math.exp ( (energy - mu) / T) ) return f #potentials in the leads muleft = -2. muright = 2. mu = 0 Tleft = 0.1 Tright = 0.1 #potentials and temepratures of all leads muTAll = [] muTAll.append([muleft,Tleft]) muTAll.append([muright,Tright]) energies = [] data = [] nenergies = 300 de = 0.032424 CurrentWave = [] CurrentLandauer = [] t = 1. W = 6 L = 10 numleads = 2 numsite = (L-1)* W sys = BuildSystem(t, W, L,mu,mu,mu) #beginning and end of stm link in terms of program; x=[] for i1 in xrange(numsite): x.append(sys.sites[i1].pos) beg=[] for i1 in xrange(numsite): if x[i1] == [1,1]: beg.append(i1) end=[] for i1 in xrange(numsite): if x[i1] == [2,1]: end.append(i1) print beg, end for imu in xrange(1): sysFinal = BuildSystem(t, W, L, mu, mu, mu) for ienergy in xrange(nenergies): energy = de * ienergy - nenergies * de / 2. energies.append(energy) wf = kwant.wave_function(sysFinal,energy) smatrix = kwant.smatrix(sysFinal, energy) #wave function based calculation which can be used for any correlation function #iterating over number of leads wfInside = numpy.zeros(numsite) #wave function inside of the system for p in xrange( numleads): wave = wf(p) lenwave = wave.shape [0] for ipsi in xrange(lenwave): wfInside = wfInside + Fermi(energy, muTAll[p][1],muTAll[p][0] ) * wave[ipsi] ProbeCurrent = wfInside[beg] * wfInside[end].conjugate() - wfInside[end] * wfInside[beg].conjugate() CurrentWave.append( ProbeCurrent[0].imag) #Getting wave function is done #transmission based calculation t10 = smatrix.transmission(1, 0) data.append(t10) t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) - Fermi(energy, muTAll[0][1],muTAll[0][0] )) CurrentLandauer.append(t2) plt.figure() plt.plot(energies, CurrentLandauer) plt.plot(energies, CurrentWave) plt.xlabel("energy [t]") plt.ylabel("conductance [e^2/h]") plt.show() def BuildSystem(t, W, L, mu, muleft, muright): sys = kwant.Builder() #DEFINING THE SYSTEM a=1. lat = kwant.lattice.square(a) for i in xrange(L-1): for j in xrange(W): 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 #ATTACHING LEADS #lead is pointing away from scattering region sym_left_lead = kwant.TranslationalSymmetry((-a,0)) left_lead = kwant.Builder(sym_left_lead) for j in xrange(W): left_lead[lat(0,j)] = muleft if j > 0: #modeling large band-width in the leads, that is why 10. t left_lead[lat(0,j), lat(0,j-1)] = - t left_lead[lat(1,j), lat(0,j)] = - t sys.attach_lead(left_lead) sym_right_lead = kwant.TranslationalSymmetry((a,0)) right_lead = kwant.Builder(sym_right_lead) for j in xrange(W): numx = L-2 right_lead[lat(numx,j)] = muright if j > 0: #modeling large band-width in the leads, that is why 10. t right_lead[lat(numx,j), lat(numx,j-1)] = - t right_lead[lat(numx,j), lat(numx+1,j)] = - t sys.attach_lead(right_lead) sysFinal = sys.finalized() t2 = time.clock() return sysFinal