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.
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 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 = []
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;
for i1 in xrange(numsite):
for i1 in xrange(numsite):
    if x[i1] == [1,1]:
for i1 in xrange(numsite):
    if x[i1] == [2,1]:
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.
        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)
        t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) - Fermi(energy, muTAll[0][1],muTAll[0][0] ))
plt.plot(energies, CurrentLandauer)
plt.plot(energies, CurrentWave)
plt.xlabel("energy [t]")
plt.ylabel("conductance [e^2/h]")

def BuildSystem(t, W, L, mu, muleft, muright):

    sys = kwant.Builder()

    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


#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

    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
    sysFinal = sys.finalized()
    t2 = time.clock()
    return sysFinal