I'm trying to simulate a scattering region confined between two reservoirs(leads) for that purpose I defines a wire with linearly increasing onsite potentials. Now I want to find the local current density at each site on the wire but I don't know how to do it. Can anyone guide me through this? :)
def make_system(a=1, t=1.0, W = 60, L = 100, pot=-0.2):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
def hop(site1, site2, B=0): # The magnetic field is controlled by the parameter B x1, y1 = site1.pos x2, y2 = site2.pos return -t * exp(0.5j * B * (y1+y2)* (x1-x2))#exp(0.5j * B * (y1-x1))
def potential_prof(i,L,pot): #potential profile linearly determined here return (L-i)*pot/L # Define the scattering region
for i in xrange(L): for j in xrange(W): # On-site Hamiltonian sys[lat(i, j)] = 4 * t+potential_prof(i,L,pot)
sys[kwant.builder.HoppingKind((1, 0), lat, lat)] = hop
sys[kwant.builder.HoppingKind((0, 1), lat, lat)] = hop
sym_left_lead = kwant.TranslationalSymmetry((-a, 0)) left_lead = kwant.Builder(sym_left_lead)
for j in xrange(W): left_lead[lat(0, j)] = 4 * t+pot if j > 0: 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): right_lead[lat(0, j)] = 4 * t if j > 0: right_lead[lat(0, j), lat(0, j - 1)] = -t right_lead[lat(1, j), lat(0, j)] = -t
sys.attach_lead(right_lead) # Modify the scattering region del sys[lat(85, 2)] del sys[lat(87, 2)] del sys[lat(86, 3)]
del sys[lat(25, 37)] del sys[lat(27, 37)] del sys[lat(26, 38)]
sys = sys.finalized()
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,
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.
I noticed in the example code you posted that you have some magnetic field in the scattering region of your system but none in the leads. I don't know anything about what precisely you may be wanting to simulate, but in general when there are abrupt changes in the Hamiltonian you'll generate some unphysical effects in your results. If you want to be able to put magnetic field in the leads as well then read on, else ignore me :).
For your case of uniform magnetic field you can choose a gauge such that the vector potential, and hence Peierls phase, is invariant by translation in the symmetry direction of your leads (if you have leads which are not all parallel to each other things become more complicated). You can then set the hoppings in the leads in exactly the same way as in the central scattering region.
It seems that the Landau gauge chosen in your example code is the correct one to use when you have leads in the x-direction (as you do).
Thanks a lot. I did it on purpose. :)