Hi Christoph,
This is a well-working code that computes conductance for graphene
with 2 leads (conductance is now commented in the end). The issue
appears at line
485.
Loading the large package warning may be due to running from Eric6
Python IDE that I use.
I have now the latest 1.4 Kwant version installed via pip3 , which
added the necessity to remove the default values from "args" parameter
list, as it's no longer supported
As for ppa install, I might be doing something wrong: (Linux Mint 19
= Ubuntu bionic)
sudo add-apt-repository ppa:kwant-project/ppa
sudo apt-get update
sudo apt-get install kwant
returns:
"Unable to locate package kwant"
The test
python3 -c 'import kwant; kwant.test(verbose=False)'
returns 21 warnings: (see below)
.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py::test_map
/home/sergey/.local/lib/python3.6/site-packages/numpy/lib/function_base.py:3652:
RuntimeWarning: Invalid value encountered in percentile
interpolation=interpolation)
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1300:
RuntimeWarning: invalid value encountered in greater
overflow_pct = 100 * np.sum(unmasked_data > new_vmax) /
len(unmasked_data)
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1301:
RuntimeWarning: invalid value encountered in less
underflow_pct = 100 * np.sum(unmasked_data < new_vmin) /
len(unmasked_data)
.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py::test_spectrum
/home/sergey/.local/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py:1069:
UserWarning: Axes3D.figure.canvas is 'None', mouse rotation disabled.
Set canvas then call Axes3D.mouse_init().
"Axes3D.figure.canvas is 'None', mouse rotation disabled. "
/home/sergey/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:83:
RuntimeWarning: invalid value encountered in reduce
return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py::test_density_interpolation
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py::test_current_interpolation
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
/home/sergey/.local/lib/python3.6/site-packages/kwant/plotter.py:1684:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
field_out[field_slice] += magns
.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py::test_current
/home/sergey/.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py:522:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
plotter.current(syst, current, file=out)
/home/sergey/.local/lib/python3.6/site-packages/kwant/tests/test_plotter.py:526:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the
future this will be interpreted as an array index, `arr[np.array(seq)]`,
which will result either in an error or a different result.
plotter.current(syst, current, ax=ax, file=out)
.local/lib/python3.6/site-packages/kwant/tests/test_wraparound.py::test_symmetry
/home/sergey/.local/lib/python3.6/site-packages/kwant/tests/test_wraparound.py:169:
RuntimeWarning: `particle_hole` and `time_reversal` symmetries are set
on the input builder. However they are ignored for the wrapped system,
since Kwant lacks a way to express the existence (or not) of a symmetry
at k != 0.
wrapped = wraparound(syst)
.local/lib/python3.6/site-packages/kwant/tests/test_wraparound.py::test_plot_2d_bands
/home/sergey/.local/lib/python3.6/site-packages/numpy/core/fromnumeric.py:83:
RuntimeWarning: invalid value encountered in reduce
return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
-- Docs: http://doc.pytest.org/en/latest/warnings.html
=================== 282 passed, 21 warnings in 51.40 seconds
===================
### Graphene ballistic edge-state field-effect transistor
from __future__ import division # so that 1/2 == 0.5, and not 0
import math
from numpy import pi, sqrt, tanh, log, sin, cos, exp
import numpy
import os
import functools as ft
import itertools as it
from matplotlib import pyplot
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import scipy.sparse as sparse
#least squres fit
import scipy.optimize as optimize
#import matplotlib.pyplot as plt
#import numpy as np
# For computing eigenvalues
import scipy.sparse.linalg as sla
import kwant
import random
from kwant.digest import gauss
#from colormaps import cmaps
path=os.getcwd()
### consider side contact of length 200 nm to a stripe of width 1000 nm
### we have a highdoping in the metallic stripe of 200 nm width and small doping in the wade graphene
### scale down by a factor of 4.5, so that a = 1 nm
scale=20/2.46
hopping_unit=3.16/scale
Wcontact = 80
narrowing=0 ## narrow the lead, relative to the size of the contact
Wsample=450
ContactDoping = 0.3/hopping_unit
#xsample=20 ###this is 2*(decay length)
xsmoothing=40
xshift=-3
flat_factor=1
xmiddle=80
### in total, the scattering region has size 2*xsample + xmiddle with two symmetric contacts
imp0=0.03/hopping_unit
impmax=0.0
impurity_width=1
edge_potential=0.0
edge_potential_width=0.8
filename="semicircle leftshift cos"
salt='abd'
###in real units, it's xsample*scale Angstroms
## (potential=-ContactDoping)
###Lead1 = contact goes from -infty to 0
###Lead2 = sample goes from xsample to infyt.
### from 0 to xsample is an intermediate matching region with smooth potential ContactDoping
def Dot2(list1,list2):
return list1[0]*list2[0]+list1[1]*list2[1]
def hopping(sitei, sitej, phi,mu,Potmax):
xi, yi = sitei.pos
xj, yj = sitej.pos
return - exp(-0.5j * phi* (xi - xj) * (yi + yj))
### make potential - here it starts after the injection region
### cos potential
#def potential(sitei, phi=0, mu=0, Potmax=0.0):
# xi, yi = sitei.pos
# pottemp=-ContactDoping + (ContactDoping - mu) *(1 - cos(xi*pi/xsample))/2.0
# if -Wcontact/2<yi<Wcontact/2 :
# yfactor=1
# elif -Wcontact/2-xsample<yi<Wcontact/2+xsample:
# shiftx=abs(yi)-Wcontact/2
# yfactor=1-(1 - cos(shiftx*pi/xsample))/2.0
# else:
# yfactor=0
# return yfactor*(pottemp+mu)-mu
### make potential - exponential decay of density with length xsample/4, which is approximately exponential decay of potential with length xsample/2
#def potential(sitei, phi=0, mu=0, Potmax=0.0):
# xi, yi = sitei.pos
# ###here I make symmetric contact:
# xi = (xsample+xmiddle/2)-abs(xi-(xsample+xmiddle/2))
# if xi<xsample:
# pottemp=-ContactDoping + (ContactDoping - mu) *(1 - cos(xi*pi/xsample))/2.0
# else:
# pottemp=-mu
# if -Wcontact/2<yi<Wcontact/2 :
# yfactor=1
# elif -Wcontact/2-xsample<yi<Wcontact/2+xsample:
# shiftx=abs(yi)-Wcontact/2
# yfactor=1-(1 - cos(shiftx*pi/xsample))/2.0
# else:
# yfactor=0
# result=yfactor*(pottemp+mu)-mu ### result without impurities
# #impurity_str=min(imp0/abs(result), impmax)
# impurity_str= impmax
# return result+ impurity_str*gauss(repr(sitei), salt)
def potential(sitei, phi, mu, Potmax):
xi, yi = sitei.pos
###here I make symmetric contact:
xi = (Wcontact/2+xmiddle/2)-abs(xi-(Wcontact/2+xmiddle/2))-xshift
r=sqrt(xi*xi*flat_factor+yi*yi)
r1=sqrt(xi*xi+yi*yi)
if r1>0:
rsmoothing=xsmoothing*r/r1
else:
rsmoothing=xsmoothing
if r<Wcontact/2:
pottemp=-ContactDoping
elif r<Wcontact/2+rsmoothing:
pottemp= -ContactDoping + (ContactDoping - mu) *(1 - cos((r-Wcontact/2)*pi/rsmoothing))/2.0
# -ContactDoping+(r-Wcontact/2)/rsmoothing *(ContactDoping-mu)
else:
pottemp=-mu
###here I make symmetric contact:
if (xi<impurity_width and abs(yi)>Wcontact/2) or abs(yi)>Wsample/2-impurity_width :
impurity_str= impmax
else :
impurity_str= 0
if xi+xshift<=edge_potential_width and abs(yi)>Wcontact/2:
extra=edge_potential
else:
extra=0
return pottemp+ impurity_str*gauss(repr(sitei), salt)+extra
def samplepot(sitei, phi, mu, Potmax):
xi, yi = sitei.pos
return -mu
def addt(tuple1,tuple2):
''' adds the two tuples '''
return tuple(map(lambda x, y: x + y, tuple1, tuple2))
def plot_conductance(sys, lead1,lead2, energies,args,label):
# Compute transmission as a function of energy
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, energy,args=args)
data.append(smatrix.transmission(lead1, lead2))
pyplot.figure()
pyplot.plot(energies*hopping_unit, data)
pyplot.xlabel("sample Fermi level [eV]")
pyplot.ylabel("conductance [e^2/h]")
pyplot.title(label)
return data
# pyplot.show()
def plot_conductance1(sys, lead1,lead2, energies,label):
### Here I plot conductance at zero energy, changing sample mu
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, 0,args=[0, energy, 0])
data.append(smatrix.transmission(lead1, lead2))
pyplot.figure()
pyplot.plot(energies*hopping_unit, data)
pyplot.xlabel("sample Fermi level [eV]")
pyplot.ylabel("conductance per spin [e^2/h]")
pyplot.title(label)
return data
def plot_conductance1s(sys, lead1,lead2, energies,label, smoothing):
### Here I plot conductance at zero energy, changing sample mu
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, 0,args=[0, energy, 0])
data.append(smatrix.transmission(lead1, lead2))
## smoothen data
estep=(energies[2]-energies[1])*hopping_unit
temperaturestep=estep/smoothing
points=numpy.array(range(-int(3/temperaturestep), 1+int(3/temperaturestep)))
smoothlist=1/(2.0 + 2.0 *numpy.cosh(points*temperaturestep)) *temperaturestep
datas=numpy.convolve(data, smoothlist, mode='same')
print(smoothlist)
pyplot.figure()
pyplot.plot(energies*hopping_unit, datas)
pyplot.xlabel("sample Fermi level [eV]")
pyplot.ylabel("conductance per spin [e^2/h]")
pyplot.title(label)
combineddata=numpy.array([(energies*hopping_unit).tolist(), datas.tolist()]).T
numpy.savetxt(path+"/"+filename+" disorder="+str(impmax)+" smoothing=" +str(2*xsmoothing)+" edgepotential="+str(edge_potential), combineddata)
print(combineddata)
print(path+"/"+filename+"disorder="+str(impmax)+" smoothing=" +str(2*xsmoothing)+" edgepotential="+str(edge_potential))
pyplot.figure()
pyplot.plot(energies*hopping_unit, data)
pyplot.xlabel("sample Fermi level [eV]")
pyplot.ylabel("conductance per spin [e^2/h]")
pyplot.title(label)
return datas
def plot_conductanceNormalized(sys, lead1,lead2, energies,args,label):
# Compute transmission as a function of energy
data = []
for energy in energies:
smatrix = kwant.smatrix(sys, energy,args=args)
smatrixNorm = kwant.smatrix(sys, energy,args=[0,args[1],0])
data.append(smatrix.transmission(lead1, lead2)/smatrixNorm.transmission(0, 1))
pyplot.figure()
pyplot.plot(energies, data)
pyplot.xlabel("energy [t]")
pyplot.ylabel("conductance Normalized by direct one [e^2/h]")
pyplot.title(label)
return data
def Eigensystem(bands,k):
mat = bands.hop * complex(math.cos(k), -math.sin(k))
mat += mat.conjugate().transpose() + bands.ham
ev,es = numpy.linalg.eigh(mat)
sortorder=numpy.argsort(numpy.absolute(ev))
return (ev[sortorder], es[:,sortorder])
def plot_wavefunction2D(flead, k, nvector=1,args=''):
"""nvector numbering is from 1 """
bands = kwant.physics.Bands(flead,args=args)
eval, evect = Eigensystem(bands,k)
evecs=evect[:,nvector-1]
wavefunction=numpy.absolute(evecs)
# xcoords = np.array([list(site.pos)[0] for site in flead.sites])[0:flead.cell_size]
ycoords = numpy.array([list(site.pos)[1] for site in flead.sites])[0:flead.cell_size]
# print xcoords, ycoords, zcoords
fig = pyplot.figure()
# ax = fig.add_subplot(1, 1, 1, projection='3d')
# The triangles in parameter space determine which x, y, z points are
# connected by an edge
pyplot.scatter(ycoords, wavefunction) #, triangles=tri.triangles, cmap=plt.cm.Spectral)
# print len(flead.sites)
# print xcoords, xcoords.shape
# pyplot.figure()
# print disorderWidth
# pyplot.scatter(xcoords,wavefunction)
pyplot.xlabel("y (lattice cells)")
pyplot.title("Mode Energy="+str(eval[nvector-1]))
def make_system(ImpurityPotential=0.1, ImpurityPotentialTop=0.05 ):
def central_region(pos): ## -0.1-1/(2*sqrt(3)) <= y <= Ly
x, y = pos
return (0<x <=Wcontact+xmiddle) and (-Wsample/2<=y <= Wsample/2)
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
lat = kwant.lattice.general([(0, 1), (cos_30, sin_30)],
[(0, 0), (1 / sqrt(3), 0)],name=('a','b'))
sys = kwant.Builder()
# types of sublattices - these are extracted by .family method: '0' and '1'
a, b = lat.sublattices
sys[lat.shape(central_region, (0,0))] = potential
## Add one point impurity on a sublattice
## sys[a(0, posimp)] = potential
sys[lat.neighbors()] = hopping
#Now add two leads in x-direction as prolongation of the layers
##chemical potential for conducting lead
# def NumberOfNeighbors(sys,site):
# n=0
# for site1 in sys.sites():
# diff=site.pos-site1.pos
# dist2=Dot2(diff,diff)
# if (dist2<1.0/3.0+0.01) and (dist2>0.00001) and (site1.family==a or site1.family==b):
# n+=1
# return n
# sitesToDelete=[]
# for site in sys.sites():
# neigh=NumberOfNeighbors(sys,site)
# if neigh<3:
# ##print neigh
# sys[site]=ImpurityPotential
# if neigh<3:
# sys[site]=ImpurityPotential
# if neigh<2:
# sitesToDelete.append(site)
# # if abs(site.pos[1]-(Ly+1/sqrt(3)))<0.05:
# # sys[site]=potential
# for site in sitesToDelete:
# del sys[site]
def connect(sys,crit):
for site1, site2 in it.product(sys.sites(), repeat=2):
if crit(site1, site2):
yield (site1, site2)
## site2 = sym.act([1],site2) # act with the symmetry
## if crit(site1, site2):
## yield (site1, site2)
kinds = [kwant.builder.HoppingKind(v, sublat) for v in [(1, 0), (0, 1),(-1,0),(0,-1),(1,-1),(-1,1)] for sublat in [a,b]]
# sys[kinds] = nexthoppingbulk
def connectlead(sys,crit,sym):
for site1, site2 in it.product(sys.sites(), repeat=2):
if crit(site1, site2):
yield (site1, site2)
site2 = sym.act([1],site2) # act with the symmetry
if crit(site1, site2):
yield (site1, site2)
# sys[connect(sys, critTop)]=nexthopping
# sys[connect(sys, critBottom)]=nexthopping
symxp = kwant.TranslationalSymmetry([sqrt(3),0])
leadxp= kwant.Builder(symxp)
def leadp_region(pos):
x, y = pos
return(-(Wcontact/2-narrowing) <= y <= (Wcontact/2-narrowing))
leadxp[lat.shape(leadp_region, (0, 0))] = -ContactDoping
leadxp[lat.neighbors()] = hopping
symxm = kwant.TranslationalSymmetry([-sqrt(3),0])
leadxm= kwant.Builder(symxm)
def leadm_region(pos):
x, y = pos
return (-(Wcontact/2-narrowing) <= y <= (Wcontact/2-narrowing))
leadxm[lat.shape(leadm_region, (0, 0))] = -ContactDoping
leadxm[lat.neighbors()] = hopping
for x in sys.sites():
x.family.norbs=1
for x in leadxp.sites():
x.family.norbs=1
for x in leadxm.sites():
x.family.norbs=1
def family_color(site):
if site.family == a:
return 'black'
elif site.family == b:
return 'blue'
else:
return 'orange'
def edgecolor(site):
return 'red' if site.family == a else 'blue'
def potential_color(site):
pot=-potential(site, phi=0, mu=0.05, Potmax=0)
if pot>ContactDoping:
pot=ContactDoping
if pot<0:
pot=0
return (pot/ContactDoping, pot/ContactDoping, pot/ContactDoping)
#this will work only for non-finalized sys
# plot
sys.attach_lead(leadxp)
sys.attach_lead(leadxm)
sysf=sys.finalized()
kwant.plot(sys,num_lead_cells=3,hop_lw=0.05,site_size=0.1,show=False,site_color=potential_color,
site_edgecolor=family_color,hop_color='black')
pyplot.show(block=False)
# print list(sys.site_value_pairs())
# print [sys[site] for site in sys.sites()]
return (leadxp.finalized(), sysf)
def plot_bandstructure(flead, momenta, args):
global trackzeroes, processors
switch=False
cutoff=0.0005
bands = kwant.physics.Bands(flead,args=args)
energies = [bands(k) for k in momenta]
for i in range(len(momenta)):
for en in energies[i]:
if abs(en)<cutoff:
if switch==False:
trackzeroes.append(momenta[i])
switch=True
else:
switch=False
pyplot.figure()
pyplot.plot(momenta, energies)
pyplot.ylim([-0.5,0.5])
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
# pyplot.show()
def compute_evs(sys,phi,evsnumber,label=''):
# Compute some eigenvalues of the closed system
sparse_mat = sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,args=[phi]))
# rank = matrix_rank(sparse_mat)
# print (sparse_mat.shape)
# mat2=sparse_mat * sparse_mat
# print mat2.shape
evs = sla.eigsh(sparse_mat,k=evsnumber, sigma=0.01, which='LM')[0]
# print evs.real
pyplot.figure()
pyplot.hist(evs.real,bins=100)
pyplot.xlabel("Energy")
pyplot.ylabel("N")
pyplot.title(label)
pyplot.show(block=False)
def density(sys, energy, args, lead_nr,label=""):
wf = kwant.wave_function(sys, energy, args)
dens= (abs(wf(lead_nr))**2).sum(axis=0)
# fig = pyplot.figure()
# Two subplots, the axes array is 1-d
nsites=len(sys.sites)
xcoords = numpy.array([list(site.pos)[0] for site in sys.sites])
ycoords = numpy.array([list(site.pos)[1] for site in sys.sites])
zcoords = numpy.array([list(site.pos)[2] for site in sys.sites])
denslayer=[[dens[i] for i in range(nsites) if (zcoords[i]==j)] for j in range(3)]
xlayer=[[xcoords[i] for i in range(nsites) if (zcoords[i]==j)] for j in range(3)]
ylayer=[[ycoords[i] for i in range(nsites) if (zcoords[i]==j)] for j in range(3)]
fig, axarr = pyplot.subplots(3, sharex='col')
for i in range(3):
# pl=axarr[i].scatter(xlayer[i], ylayer[i], c=denslayer[i], s=50*denslayer[i], cmap=cm.jet)
pl=axarr[i].hist(ylayer[i], weights=denslayer[i])
# fig.colorbar(pl)
axarr[i].set_title(label+str(energy)+" layer "+str(i))
# plt.scatter(x, y, c=y, s=500, cmap='gray')
#kwant.plotter.map(sys, onsite,show=False)
# energy = 0.75
def onsite_plot():
#plot the onsite potential
impuritypot= numpy.array([sys.hamiltonian(site,site, 0 ,"",ndonors,cutoff,ddist,alpha,phiosc,losc ) for site in sitenumbers])
fig = pyplot.figure()
ax = fig.gca(projection='3d')
p=ax.plot_trisurf(xcoords,ycoords,impuritypot,cmap=cm.jet, linewidth=0.2)
ax.set_title("onsite potential")
HoppingWidth=1
chempot=0.0
nexthopping=0.05
## Make a normalization system to eliminate the problems of contacts
#(lead,sysNorm) = make_systemNorm( Lx=Lx, Ly=Ly, HoppingWidth=HoppingWidth, LeadEnergy=1.9, HoppingFromLead=0.05, posimp=0, mu=chempot, nexthopping=nexthopping,nexthoppingbulk=-0.00)
#plot_wavefunction2D(lead, 3.5, nvector=1)
#plot_wavefunction2D(lead, 3.5, nvector=2)
def EdgeShape(x):
s=0
for Defect in EdgeDefects:
s+=Defect[0]*(numpy.arctan((x-Defect[2])/Defect[1])/pi+0.5)
return s
##def EdgeShape(x):
## s=0
## for Defect in EdgeDefects:
## if abs(x-Defect[2])<Defect[1]:
## s+=Defect[0]*sqrt(Defect[1]*Defect[1]-(Defect[2]-x)*(Defect[2]-x))
## return s
(lead, sys) = make_system()
#plot_wavefunction2D(lead, k=-0.2, nvector=1,args=[phi,chempot,0])
#plot_wavefunction2D(lead, k=-0.2, nvector=2,args=[phi,chempot,0])
#plot_wavefunction2D(lead, k=-0.3, nvector=1,args=[phi,chempot,0])
#plot_wavefunction2D(lead, k=-0.3, nvector=2,args=[phi,chempot,0])
#plot_wavefunction2D(lead, k=-0.4, nvector=1,args=[phi,chempot,0])
#plot_wavefunction2D(lead, k=-0.4, nvector=2,args=[phi,chempot,0])
#xcoords = numpy.array([list(site.pos)[0] for site in sys.sites])
#ycoords = numpy.array([list(site.pos)[1] for site in sys.sites])
## Find the number and position of top and bottom edge sites
#momenta = numpy.linspace(1.8,2.5,500)
momentafull = numpy.linspace(0,2*pi,100)
#plot_bandstructure(sys.leads[0], momenta,args=[phi,chempot,0])
trackzeroes=[]
phi=0
chempot=0.0
plot_bandstructure(lead, momentafull,args=[phi,chempot,0])
energies = numpy.linspace(-0.201,0.2,200)
pyplot.show(block=False)
#plot_conductance1s(sys, 0,1, energies,label="", smoothing=0.0015)
J_0 = kwant.operator.Current(sys)
wf = kwant.solvers.default.wave_function(sys, energy=0.01, args=[phi,chempot,0])
psi=wf(1)
current = J_0(psi)
#kwant.plotter.current(sys, current, relwidth=0.05, args=[phi,chempot,0])
pyplot.show()
