Thank you very much for the prompt reply.
I expect the current change linearly by increasing the bias.
I have attached the code.
Thank you very much for your help.
Best Regards,
Arad
On Thu, May 4, 2017 at 2:20 PM, Joseph Weston <joseph.westo...@gmail.com>
wrote:
> Hi,
>
> > 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.
>
> If you want help you are going to need to be more specific I am afraid.
> What do you expect, and what does your program give you?
>
> I couldn't run the code you dumped into the body of the email just by
> copy/pasting it into a text file because the formatting is all messed
> up. Please submit any scripts as actual attachments which run and
> produce what you wish to show. Only inline *relevant snippets* into the
> body of the email; not a whole script. This greatly reduces the amount
> of effort people have to go through to help out.
>
> Happy Kwanting,
>
> Joe
>
import numpy as np
from matplotlib import pyplot
import kwant
import math
from scipy import exp
from scipy.integrate import quad
t=1
L0 = 30
dletaV = 0.025 #(Steps to change bias)
deltaE = 0.03 #(steps of energy)
beta= 1.0/ 0.025
Itot = []
right_fermies = []
right_fermiesEV = []
current=[]
def make_system(a=1, t=1 , W= 10):
lat = kwant.lattice.square(a)
sys = kwant.Builder()
#by applying bias the energy of each point will change linearly
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)
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)
#the energy of the the right lead will change by bias while the left lead is fixed
def bias_R (site, volt):
(x, y) = site.pos
if x == 0:
return 4 * t -volt
else:
return 0
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 = []
#calculate the conductance for each bias point for energy in energies:
for energy in energies:
smatrix = kwant.smatrix(sys, energy, args=[right_fermi])
conductance.append(smatrix.transmission(1,0))
Iie = 0
#integrate the corresponding conductance from zero(mu left) to energy=right_fermi (mu right).
intvolt= int((right_fermi)/((deltaE )))
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 ))) + 1.0)) - ( 1.0 / (math.exp((1.0/ 0.025)*((deltaE * s) + right_fermi)) + 1.0))
#f_l - fr for conductance[s+1]
F2 = (1.0 / (math.exp((beta)* ((deltaE * (s +1)))) + 1.0))-(1.0 / (math.exp((1.0/ 0.025)*((deltaE * (s +1)) + right_fermi)) + 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()
