Dear Christoph,
Here is the minimal example giving a error I described. The system
is a 3-layer graphene and I add a couple of vertical 1D chains with
small hopping to the
main sample to model the STS tip.
Thank you,
Sergey
On 23/09/2020 15:52, Christoph Groth wrote:
Sergey Slizovskiy wrote:
I am having a error with "gauge" when I am adding 1D leads (Chains)
to my 3D system. (I have both Chain leads, where I do not need gauge,
of course,
Could you please provide a minimal self-contained script that allows to
reproduce the problem?
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'))
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 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, 0, -1),a(0, 0, -2)] = tipleadhopping
tipleadb[b(0, 0, -2)] = 0
tipleadb[b(0, 0, -1), b(0, 0, -2)] = tipleadhopping
sys.attach_lead(tipleada)
sys.attach_lead(tipleadb)
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()
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)
gauge = kwant.physics.magnetic_gauge(sys)
