Hi dear developers,
My question is about calculating total current for a finite bias in a
nanowire.
I looked at the Datta book and try to implement it in KWANT but I can not
get
a reliable answer, I think I made a mistake somewhere but I can not find it.
I define a nanowire and assume that by applying the bias at the ends of
that, the energy of each point will decrease linearly.
To calculate the current for each bias point I did the following:
For each value of bias I change the energy of the system, calculate
conductance, and integrate the corresponding conductance
from zero(mu left) to energy=bias (mu left).
Would you please help me, I don't know what I am doing wrong.
Thank you so much,
Arad
import numpy as np
from matplotlib import pyplot
import kwant
import math
from scipy import exp
from scipy.integrate import quad
#the effective mass for InAs is o.o23m0 and a=1nm
t=1.65
L0 = 30
dletaV = 0.025
deltaE = 0.03
beta= 1.0/ 0.025
#becouse of the doping the Fermi level is higher then the conductance band:
Ef-Ec= 0.25
delta=0.25
Itot = []
right_fermies = []
right_fermiesEV = []
current=[]
def make_system(a=1, t=1.65 , W= 10):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
def bias (site, volt):
(x, y) = site.pos
if 0 <= x < L0:
return (volt /(L0)) * (x)
elif x == L0:
return (volt)
else:
return 0
def onsite(site, volt=0):
return 4 * t -bias (site, volt)
def bias_R (site, volt):
(x, y) = site.pos
if x == 0:
return 4 * t -volt
else:
return 0
sys[(lat(x, y) for x in range(0,L0) for y in range(W))] = onsite
sys[lat.neighbors()] = -t
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
if j > 0:
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):
right_lead[lat(0, j)] = bias_R
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)
return sys
def cal_current (sys, energies, right_fermi):
conductance = []
for energy in energies:
smatrix = kwant.smatrix(sys, energy, args=[right_fermi])
conductance.append(smatrix.transmission(1,0))
'''
pyplot.figure()
pyplot.plot(energies, conductance)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.show()
'''
Iie = 0
intvolt= int((right_fermi)/((deltaE )*t))
for s in xrange(intvolt):
#0.77*0.0001 is e^2/h
#f_l - fr for conductance[s]
F1 = (1.0 / (math.exp((beta)* ((deltaE * s * t ) - delta)) + 1.0))
- ( 1.0 / (math.exp((1.0/ 0.025)*((deltaE * s * t) + right_fermi - delta))
+ 1.0))
#f_l - fr for conductance[s+1]
F2 = (1.0 / (math.exp((beta)* ((deltaE * (s +1) * t) - delta)) + 1.0
))-(1.0 / (math.exp((1.0/ 0.025)*((deltaE * (s +1) * t ) + right_fermi -
delta)) + 1.0))
H = ((conductance[s]*F1) + (conductance[s +1])*F2)* (deltaE/2) *
0.77 * 0.0001
Iie = H + Iie
right_fermiesEV.append(right_fermi )
Itot.append(Iie)
right_fermies.append(right_fermi)
return Itot
def main():
for j in xrange(20):
sys = make_system()
sys = sys.finalized()
#becouse of the doping the integral is not till o (the left Fermi
level) it is till 0.25
current = cal_current(sys, energies=[((deltaE * i)) for i in xrange(
25) ], right_fermi= ((dletaV * j)))
np.savetxt('Rdatay.dat', Itot)
np.savetxt('Rdatax.dat', right_fermies)
pyplot.figure()
pyplot.plot(right_fermiesEV, current)
pyplot.xlabel("bias [ev]")
pyplot.ylabel("current [A]")
pyplot.show()
if __name__ == '__main__':
main()