Hi Christoph,
Here is my code.
I did all as you described, you may try to run it yourself. Must take
around a minute.
Best wishes,
Sergey
On 25/09/2020 10:51, Christoph Groth wrote:
Sergey,
I’m a bit puzzled by the fact that your profile contains traces of
Cython cdef functions. They shouldn’t be visible to the Cython
profiler. Would you mind sharing the script that you profiled (or
something equivalent, but then with another profile)?
Cheers
Christoph
### Graphene ballistic edge-state field-effect transistor
from __future__ import division # so that 1/2 == 0.5, and not 0
import math
import numpy as np
from numpy import pi, sqrt, tanh, log, sin, cos, exp
import os
#import functools as ft
import itertools as it
from matplotlib import pyplot as plt
#from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import scipy.sparse as sparse
import scipy.sparse.linalg as sla
import kwant
#import random
from kwant.digest import gauss
#from colormaps import cmaps
from joblib import Parallel, delayed
path=os.getcwd()
sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
lat = kwant.lattice.general([(1, 0, 0), (sin_30,cos_30, 0), (0,0, 2)],
[(0, 0, 0), (0,1 / sqrt(3), 0), (0, 0, 1), (1 /(2), 1/(2*sqrt(3)), 1)],name=('a','b','c','d'))
## to make it more clear, I tilt the layers, preserving the same hoppings. Now, g1 will look skew
#lat = kwant.lattice.general([(0, 1, 0), (cos_30, sin_30, 0), (0, 0, 1)],
# [(0, 0, 0), (1 / sqrt(3), 0, 0)],name=('a','b','d','c'))
fakelat= kwant.lattice.Monatomic([(1, 0, 0), (sin_30,cos_30, 0), (0,0, 2)], [0,0,0.1])
sys = kwant.Builder()
# types of sublattices - these are extracted by .family method: '0' and '1'
a, b, c, d = lat.sublattices
#delta=0
#g2=0
#g3=0
#g4=0
def Dot2(list1,list2):
return list1[0]*list2[0]+list1[1]*list2[1]
def g0hopping(sitei, sitej, peirels):
xi, yi, zi = sitei.pos
xj, yj ,zj= sitej.pos
return -g0 * peirels(sitei,sitej)
def g3hopping(sitei, sitej, peirels):
xi, yi, zi = sitei.pos
xj, yj ,zj= sitej.pos
return -g3 * peirels(sitei,sitej)
def g4hopping(sitei, sitej, peirels):
xi, yi, zi = sitei.pos
xj, yj ,zj= sitej.pos
return g4 * peirels(sitei,sitej)
def addt(tuple1,tuple2):
''' adds the two tuples '''
return tuple(map(lambda x, y: x + y, tuple1, tuple2))
def Eigensystem(bands,k):
mat = bands.hop * complex(math.cos(k), -math.sin(k))
mat += mat.conjugate().transpose() + bands.ham
ev,es = np.linalg.eigh(mat)
sortorder=np.argsort(np.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=np.absolute(evecs)
# xcoords = np.array([list(site.pos)[0] for site in flead.sites])[0:flead.cell_size]
ycoords = np.array([list(site.pos)[1] for site in flead.sites])[0:flead.cell_size]
# print xcoords, ycoords, zcoords
fig = plt.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
plt.scatter(ycoords, wavefunction) #, triangles=tri.triangles, cmap=plt.cm.Spectral)
# print len(flead.sites)
# print xcoords, xcoords.shape
# plt.figure()
# print disorderWidth
# plt.scatter(xcoords,wavefunction)
plt.xlabel("y (lattice cells)")
plt.title("Mode Energy="+str(eval[nvector-1]))
eps=0.001
def TipPotential(sitei, depth, width,gate):
x, y, z=sitei.pos
if (sitei.family==a) or (sitei.family==c):
ad=DeltaAB+g2-Ef
else:
ad = g2-Ef
if z==0:
return depth*exp(-((x-shift)**2+y*y)/(2*width*width)) + ad
else:
if z==2:
return gate + ad
else:
return gate/3 + ad +U2/1000.0 + depth*exp(-((x-shift)**2+y*y)/(2*width*width))/3
def LeadPot(sitei, gate):
x, y, z=sitei.pos
if (sitei.family==a) or (sitei.family==c):
ad=DeltaAB+g2-Ef
else:
ad = g2-Ef
if z==0:
return ad
else:
if z==2:
return gate + ad
else:
return gate/3 + ad
def make_system(size: 50):
def connect(sys, crit):
for site1, site2 in it.product(sys.sites(), repeat=2):
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 critInterlayer(site1, site2): ## gamma1 hopping integral
diff = site1.pos - site2.pos
[x, y] = diff
if (0.01 < abs(x) + abs(y) < 0.1): # and (site1.tag!=site2.tag):
return True
else:
return False
## armchair stripe
sys = kwant.Builder()
def sys_region(pos):
x, y, z = pos
return (x * x + y * y <= size*size) and (0 <= z <= Nl-1)
sys[lat.shape(sys_region, (0, 0,0))] = TipPotential
hopping0 = (((0, 0,0), a, b), ((0, 1,0), a, b), ((-1, 1,0), a, b), ((0, 0,0), d, c),((0, -1,0), d, c),((-1, 0,0), d, c))
hoppingg1 = (((0, 0, 0), a, c),((0, 0, 1), a, c))
hoppingg4 = (((0, 0, 0), b, c),((0, 0, 1), b, c), ((1, -1, 0), b, c),((1, -1, 1), b, c), ((0, -1, 0), b, c), ((0, -1, 1), b, c),
((0, 0, 0), a, d),((0, 0, 1), a, d), ((1, 0, 0), a, d),((1, 0, 1), a, d), ((0, 1, 0), a, d), ((0, 1, 1), a, d))
hoppingg3 = (((0, 0, 0), b, d),((0, 0, 1), b, d), ((1, 0, 0), b, d),((1, 0, 1), b, d), ((1, -1, 0), b, d),((1, -1, 1), b, d))
hoppingg2 = (((0, 0, 1), b, b),((0, 0, 1), d, d))
hoppingg5 = (((0, 0, 1), c, c), ((0, 0, 1), a, a))
print("setting hoppings")
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hopping0]] = g0hopping
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg1]] = g1
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg4]] = g4hopping
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg3]] = g3hopping
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg2]] = g2/2
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg5]] = g5/2
### add tip sites and hoppings from them
tiphopping=0.001
tipleadhopping=1
sys[a(0, 0, -1)] = 0
sys[a(0,0,-1),a(0,0,0)] = tiphopping
sys[b(0, 0, -1)] = 0
sys[b(0, 0, -1), b(0, 0, 0)] = tiphopping
tipleada=kwant.Builder(kwant.TranslationalSymmetry([0,0,-2]))
tipleadb = kwant.Builder(kwant.TranslationalSymmetry([0, 0, -2]))
tipleada[a(0,0,-2)]=0
tipleada[a(0, 1, -2)] = 0
tipleada[a(0, 1, -2),a(0, 0, -2)] = 0
tipleada[a(0, 0, -1),a(0, 0, -2)] = tipleadhopping
tipleadb[b(0, 0, -2)] = 0
tipleadb[b(0, 0, -1), b(0, 0, -2)] = tipleadhopping
tipleadb[b(0, 1, -2)] = 0
tipleadb[b(0, 1, -2),b(0, 0, -2)] = 0
Wlead=size/2-1.5
symvecs=[[1, 0, 0], [sin_30,cos_30, 0], [1-sin_30,-cos_30, 0], [-1, 0, 0], [-sin_30,-cos_30, 0], [-1+sin_30,cos_30, 0]]
leadsymm = [kwant.TranslationalSymmetry(vec) for vec in symvecs]
leads = [kwant.Builder(sym) for sym in leadsymm]
leadregions=[lambda pos: ((-Wlead <= vec[0]*pos[1]-vec[1]*pos[0] <= Wlead) and (0 <= pos[2] <= Nl-1)) for vec in symvecs]
sys.attach_lead(tipleada)
sys.attach_lead(tipleadb)
for l in range(6):
leads[l][lat.shape(leadregions[l],(0,0,0))]=LeadPot
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hopping0]] = g0hopping
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg1]] = g1
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg4]] = g4hopping
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg3]] = g3hopping
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg2]] = g2 / 2
leads[l][[kwant.builder.HoppingKind(*hopping) for hopping in hoppingg5]] = g5 / 2
sys.attach_lead(leads[l].substituted(peirels='peirels_lead'+str(l)))
# sys.attach_lead(leads[1])
print("hoppings ready")
for x in sys.sites():
x.family.norbs = 1
def family_color(site):
if site.family == a:
return 'black'
elif site.family == b:
return 'blue'
elif site.family==c:
return 'gray'
else:
return 'cyan'
def edgecolor(site):
if site.family == a:
return 'black'
elif site.family == b:
return 'blue'
elif site.family==c:
return 'gray'
else:
return 'cyan'
# axis=matplotlib.axes.Axes()
kwant.plot(sys, num_lead_cells=3, hop_lw=0.02, site_size=0.1, show=False, site_color=family_color, site_edgecolor=edgecolor,
hop_color='black') # , ax=axis)
#kwant.plot(leads[1], num_lead_cells=3, hop_lw=0.02, site_size=0.1, show=False, site_color=family_color,
# site_edgecolor=edgecolor, hop_color='black') # , ax=axis)
plt.show(block=False)
# print list(sys.site_value_pairs())
# print [sys[site] for site in sys.sites()]
return sys.finalized()
def gauss(x,broad):
if np.abs(x)>6*broad:
return 0
else:
return exp(-x*x/(2*broad*broad))/(broad*sqrt(2*pi))
def compute_ldos(sys,phi, depth, width,gate):
# Compute some eigenvalues of the closed system
Ham_mat = sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gate':gate})
#sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gap':gap}))
#this gives the indices of A and B sites in the middle of the circular region
siteA= sys.id_by_site[a(0, 0, 0)]
siteB = sys.id_by_site[b(0, 0, 0)]
ev, es = sla.eigsh(Ham_mat,k=4000,sigma=0.001, which='LM')
weightA=pow(abs(es[siteA]),2)
weightB=pow(abs(es[siteB]),2)
epoints=np.linspace(-0.1,0.1,400)
DOSA= [sum([w*gauss(e-ev,0.0015) for (w,ev) in zip(weightA,ev)]) for e in epoints]
DOSB = [sum([w * gauss(e - ev, 0.0015) for (w, ev) in zip(weightB, ev)]) for e in epoints]
# print(ev)
# print(weightA)
# selection=weightA>0.005
# print(np.stack((ev,weightA)).transpose()[selection])
# histA,bin_edges = np.histogram(ev[selection],bins=50,range=(-0.5,0.5),weights=weightA[selection])
# histB,bin_edges = np.histogram(ev[selection], bins=50, range=(-0.5, 0.5), weights=weightB[selection])
#print(epoints,DOSA)
return [DOSA,DOSB]
def compute_ldos1(sys,phi, depth, width,gate):
# Compute some eigenvalues of the closed system
Ham_mat = sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gate':gate})
#sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gap':gap}))
#this gives the indices of A and B sites in the middle of the circular region
siteA= sys.id_by_site[a(0, 0, 0)]
siteB = sys.id_by_site[b(0, 0, 0)]
def LDOS(e):
ev, es = sla.eigsh(Ham_mat,k=30,sigma=e, which='LM')
weightA=pow(abs(es[siteA]),2)
weightB=pow(abs(es[siteB]),2)
return [sum([w * gauss(e - ev, 0.0015) for (w, ev) in zip(weightA, ev)]),sum([w * gauss(e - ev, 0.0015) for (w, ev) in zip(weightB, ev)])]
epoints=np.linspace(-0.1,0.1,100)
DOS = [LDOS(e) for e in epoints]
return np.array(DOS).T
# print(ev)
# print(weightA)
# selection=weightA>0.005
# print(np.stack((ev,weightA)).transpose()[selection])
# histA,bin_edges = np.histogram(ev[selection],bins=50,range=(-0.5,0.5),weights=weightA[selection])
# histB,bin_edges = np.histogram(ev[selection], bins=50, range=(-0.5, 0.5), weights=weightB[selection])
#print(epoints,DOSA)
return [DOSA,DOSB]
def compute_ldos2(sys,phi, depth, width,gate):
# Compute some eigenvalues of the closed system
Ham_mat = sys.hamiltonian_submatrix(sparse=False,params={'phi':phi,'depth':depth,'width':width,'gate':gate})
#sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gap':gap}))
#this gives the indices of A and B sites in the middle of the circular region
siteA= sys.id_by_site[a(0, 0, 0)]
siteB = sys.id_by_site[b(0, 0, 0)]
ev, es = np.linalg.eigh(Ham_mat)
weightA=pow(abs(es[siteA]),2)
weightB=pow(abs(es[siteB]),2)
epoints=np.linspace(-0.1,0.1,400)
DOSA= [sum([w*gauss(e-ev,0.0015) for (w,ev) in zip(weightA,ev)]) for e in epoints]
DOSB = [sum([w * gauss(e - ev, 0.0015) for (w, ev) in zip(weightB, ev)]) for e in epoints]
def compute_ldos3(sys, phi, depth, width, gate,epoints:np.linspace(-0.1, 0.1, 400)):
# Compute some eigenvalues of the closed system
print("computing LDOS for fi="+str(phi))
#Ham_mat = sys.hamiltonian_submatrix(sparse=True, params={'phi': phi, 'depth': depth, 'width': width, 'gate': gate})
# sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gap':gap}))
# this gives the indices of A and B sites in the middle of the circular region
LeadsB=[[0,0,phi] for i in range(6)]
peirels_sysLeadList=list(gauge([0,0,phi],*LeadsB))
siteA = sys.id_by_site[a(0, 0, 0)]
siteB = sys.id_by_site[b(0, 0, 0)]
lead_params={'peirels_lead'+str(i):peirels_sysLeadList[i+1] for i in range(6)}
params = {'peirels':peirels_sysLeadList[0], 'depth': depth, 'width': width, 'gate': gate, **lead_params}
DOS = [kwant.ldos(sys,e,params=params)[[siteA,siteB]] for e in epoints]
return np.array(DOS).T
# print(ev)
# print(weightA)
# selection=weightA>0.005
# print(np.stack((ev,weightA)).transpose()[selection])
# histA,bin_edges = np.histogram(ev[selection],bins=50,range=(-0.5,0.5),weights=weightA[selection])
# histB,bin_edges = np.histogram(ev[selection], bins=50, range=(-0.5, 0.5), weights=weightB[selection])
#print(epoints,DOSA)
return [DOSA,DOSB]
def compute_conductance(sys, phi, depth, width, gate,epoints:np.linspace(-0.1, 0.1, 400)):
# Compute some eigenvalues of the closed system
print("computing LDOS for fi="+str(phi))
#Ham_mat = sys.hamiltonian_submatrix(sparse=True, params={'phi': phi, 'depth': depth, 'width': width, 'gate': gate})
# sparse.csc_matrix(sys.hamiltonian_submatrix(sparse=True,params={'phi':phi,'depth':depth,'width':width,'gap':gap}))
# this gives the indices of A and B sites in the middle of the circular region
LeadsB=[[0,0,0],[0,0,0]]+[[0,0,phi] for i in range(6)]
peirels_sysLeadList=list(gauge([0,0,phi],*LeadsB)) ### gauges for system and leads
siteA = sys.id_by_site[a(0, 0, 0)]
siteB = sys.id_by_site[b(0, 0, 0)]
lead_params={'peirels_lead'+str(i):peirels_sysLeadList[i+3] for i in range(6)}
params = {'peirels':peirels_sysLeadList[0], 'depth': depth, 'width': width, 'gate': gate, **lead_params}
dat=[]
for en in epoints:
smatrix = kwant.smatrix(sys, energy=en, params=params, in_leads=[0,1],out_leads=[3])
dat.append([smatrix.transmission(lead_out=3, lead_in=0), smatrix.transmission(lead_out=3, lead_in=1)])
return np.array(dat)
# print(ev)
# print(weightA)
# selection=weightA>0.005
# print(np.stack((ev,weightA)).transpose()[selection])
# histA,bin_edges = np.histogram(ev[selection],bins=50,range=(-0.5,0.5),weights=weightA[selection])
# histB,bin_edges = np.histogram(ev[selection], bins=50, range=(-0.5, 0.5), weights=weightB[selection])
#print(epoints,DOSA)
return [DOSA,DOSB]
path = os.getcwd()
scaling=30
Btophi=0.0000919035*pow(scaling,2)
g0=3.16/scaling
g1=0.39
g4=0.14/scaling
g3=0.315/scaling
g2=-0.020
g5=0.044
DeltaAB=0.02# 0.02
Ef=-0.0213
Nl=3
sys=make_system(size=5)
print("system is ready")
gauge = kwant.physics.magnetic_gauge(sys)
print("gauge is ready")
