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 = []
muTAll.append([muleft,Tleft])
muTAll.append([muright,Tright])
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;
x=[]
for i1 in xrange(numsite):
x.append(sys.sites[i1].pos)
beg=[]
for i1 in xrange(numsite):
if x[i1] == [1,1]:
beg.append(i1)
end=[]
for i1 in xrange(numsite):
if x[i1] == [2,1]:
end.append(i1)
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.
energies.append(energy)
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)
data.append(t10)
t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) -
Fermi(energy, muTAll[0][1],muTAll[0][0] ))
CurrentLandauer.append(t2)
plt.figure()
plt.plot(energies, CurrentLandauer)
plt.plot(energies, CurrentWave)
plt.xlabel("energy [t]")
plt.ylabel("conductance [e^2/h]")
plt.show()
def BuildSystem(t, W, L, mu, muleft, muright):
sys = kwant.Builder()
#DEFINING THE SYSTEM
a=1.
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
#ATTACHING LEADS
#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
sys.attach_lead(left_lead)
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
sys.attach_lead(right_lead)
sysFinal = sys.finalized()
t2 = time.clock()
return sysFinal