Good morning, Joe,
  Thank you, then it looks as I have done it right,  but I was confused as I got slightly asymmetric current plot from symmetric geometry.  Then, the problem might be in asymmetric choice of starting points for streamplot.
I am trying to manually add the symmetric set of starting points.  My code is:

J_0 = kwant.operator.Current(sys)
seed_points = numpy.array([ [x, y]  for y in numpy.nditer(numpy.linspace(-Wsample/2, Wsample/2, 40)) for x in numpy.nditer(numpy.linspace(5 , Wcontact+xmiddle-5, 3))])

for en in [0.001, 0.01, 0.02]:
    wf = kwant.solvers.default.wave_function(sys, energy=en, args=[phi,chempot,0])
    currents = [J_0(p,  args=[phi,chempot,0]) for p in psi]
    kwant.plotter.current(sys, current, relwidth=0.05, start_points=seed_points)

I have two small questions here:
1)  This code gives a strange error:

The debugged program raised the exception unhandled TypeError
"streamplot() got an unexpected keyword argument 'start_points'"
File: /home/sergey/.local/lib/python3.6/site-packages/kwant/, Line: 2164

which is strange since start_points  is in the streamplot documentation

2)  How to make several plots be generated without manual closure of the previous one?  My command  does not seem to work in this case.


On 19/10/18 10:01, Joseph Weston wrote:
Good morning,

Hi Joe,  Thank you, now it works!

  But I am a bit confused about the physical meaning of the current,
created by modes in lead n.   How to plot the real current between the
two leads?  Imagine,  as the simplest case, that my geometry is
inversion-symmetric,  so that lead 0 goes into lead 1 under
inversion.  Then I expect the current to be inversion-symmetric as well.
Summing the current contributions for *all* the scattering states
originating from lead L at an energy E gives you the current you would
measure if the fermi levels in all the leads are at energy E and you
apply an infinitesimal voltage to lead L.

For example if we want to see the current profile if we add an
infinitesimal voltage dV to lead 0:

    wfs = kwant.wave_function(syst, energy=E, params=...)
    J_0 = dV * sum(J(psi) for psi in wfs(0))

This is analogous to the transmission obtained by kwant.smatrix:

    smatrix = kwant.smatrix(syst, energy=E, params=...)
    I_10 = dV * smatrix.transmission(1, 0)

I_10 is the current we would measure in lead 1 after applying dV to lead 0.

In the above I have elided the e^2/h factor for brevity, but hopefully I
have been clear enough to get the point across.

Hope that helps,