Thank You Adel,

Thank you for helping on previous issue. I have another curiosity here. Right now all incoming flux on lead 0 is all modes from 1 to N. The incoming flux corresponds to summation from 1 to N (i.e  ), I just want first 3 modes only. 

I have attached a picture here what I want to do it. Can I achieve such condition in kwant from finding wavefunction from kwant solver?

​Here I have attached a code to find wavefunction. 

import kwant

import numpy as np

from matplotlib import pyplot

from numpy import sqrt

from math import *

#======================================================================

# Define the shape -------------------

#======================================================================

def shape(pos):

    x, y = pos

    return (np.abs(2.0*x)<=L)&(np.abs(y)<W1/2.0)

#----------------------------------------------------------------------

#======================================================================

def wshape(pos):

    x, y = pos

    return (np.abs(2.0*x)<=L)&(np.abs(y)<W1/2.0)                          

#----------------------------------------------------------------------

a      = 1;

t      = 1;

E0L    = 0*t

L      = 10;  # Length of the

W1   = 30;  # Width on the left

# Define geometry

#--------------------------------------------------------------------------

sys0 = kwant.Builder()

lat  = kwant.lattice.square(a)

#--------------------------------------------------------------------------------------------

# Define onsite energies and couplings

#----------------------------------------------------

sys0[lat.shape(wshape,(0,0))] = E0     # To make all sites the same

sys0[lat.neighbors()]           = t

#--------------------------------------------------------------------------------------------

# Left lead

#----------------------------------------------------------------------------------

left_lead = kwant.Builder(kwant.TranslationalSymmetry([-1,0]))

left_lead[(lat(0,y) for y in range(int(-W1/2+1),int(W1/2)))] = E0L

left_lead[lat.neighbors()] = t

sys0.attach_lead(left_lead);

#--------------------------------------------------------------------------------------------

# Right lead

#---------------------------------------------------------------------------------

right_lead = kwant.Builder(kwant.TranslationalSymmetry([1,0]))

right_lead[(lat(0,y) for y in range(int(-W1/2+1),int(W1/2)))] = E0L

right_lead[lat.neighbors()] = t

sys0.attach_lead(right_lead);

#---------------------------------------------------------------------------------------------

sys = sys0.finalized()

kwant.plot(sys);

wf = kwant.solvers.default.wave_function(sys,E,check_hermiticity=False); wf0 = wf(0); wf1 = wf(1);

def plot_conductance(sys, energies):

    data = []

    for energy in energies:

        smatrix = kwant.smatrix(sys, energy)

        data.append(smatrix.transmission(1, 0))

        #transmission from left 2 leads to right lead

    pyplot.figure()

    pyplot.plot(energies, data)

    pyplot.xlabel("energy [t]")

    pyplot.ylabel("conductance [e^2/h]")

    pyplot.show()

energies=[-3+i*0.02 for i in range(100)]

plot_conductance(sys,energies)