Hi, You can calculate the particle current from site "j" to site "i" due to a state "psi": I_ij = -2 Im(psi_i* H_ij psi_j) Where "Im" is the imaginary part, "H_ij" is the Hamiltonian matrix element between sites "i" and "j" and "*" is complex conjugation. H_ij can be obtained from your finalized system, "sys", with `sys.hamiltonian(i, j, args=args)`. The wavefunctions "psi" can be obtained from `kwant.wave_function(sys, args)`. The following Python snippet should give you what you want:: lat = kwant.lattice.square() B = 1.0 # your magnetic field E = 1.0 # energy @ which you want the current sys = make_system() wf = kwant.wave_function(sys, energy=E, arg=(B,)) site_i = lat(4, 5) # whatever sites you want site_j = lat(4, 6) # get the indices of the sites in the wavefunction i = sys.sites.index(site_i) j = sys.sites.index(site_j) def current(psi): H_ij = sys.hamiltonian(i, j, args=(B,)) return -2 * (psi[i].conjugate() * H_ij * psi[j]).imag I_ij = 0 # add the contributions from all open modes in all leads for l, lead in enumerate(sys.leads): for psi in wf(l): I_ij += current(psi) There are several caveats to the above: 1. It will not work if we have sites with more than 1 orbital as we have assumed that there is a 1-1 mapping between sites and elements in "psi". 2. No statistical physics has been included in the above -- the quantity calculated is *not* strictly the particle current, but the "particle current at energy E". To get the actual particle current you need to do the statistical physics right. This will involve integrating over energy and weighting the contributions from each lead by some distribution function (e.g. Fermi-Dirac) to get the full particle current. Depending on what kind of system you are studying you may be able to justify foregoing the energy integration on physical grounds, but this is a question of the physics you are looking at and is independent of kwant. Hope that helps, Joe P.S. I have not tested the above code snippet, and it's entirely possible that I've got a sign wrong, missed a few syntax errors etc. The point is just to give you the gist of how you would go about implementing such a calculation.