# Physical background
#-------------------------------------
# Calculate the Hofstadter's Butterfly of G-BN
from __future__ import division # so that 1/2==0.5, and not 0
from math import pi, sqrt,cos,sin
from cmath import exp
import numpy as np
import kwant
#For plotting
from matplotlib import pyplot
#Define the graphnee lattice
sin_30, cos_30=(1 / 2, sqrt(3) / 2)
#Define the scale
a_tb=1.0
gamma=5
L=20*gamma
W=5*gamma*sqrt(3)
ab=4*(pi/sqrt(3))*gamma
En=0.05
graphene=kwant.lattice.general([(1, 0),(sin_30, cos_30)],\
[(0, -1/(4*cos_30)),(0, 1/(4*cos_30))])
a, b = graphene.sublattices
class SimpleNamespace(object):
def __init__(self, **kwargs):
self.__dict__.update(kwargs)
def make_system(L, W, mu_lead=2):
#### Define the scattering region. ###########
# slab scattering region
def slab_shape(pos):
(x, y)=pos
return 0 <= x <= L and 0 <= y <= W
# Define on-site energy
def potentiala(site):
(x, y)=site.pos
da1=cos((-sqrt(3)*x/2+y/2)*ab)+cos(-y*ab)+\
cos((sqrt(3)*x/2+y/2)*ab)
da2=sin((-sqrt(3)*x/2+y/2)*ab)+sin(-y*ab)+\
sin((sqrt(3)*x/2+y/2)*ab)
return (2*pi/gamma)*(da1-sqrt(3)*da2)*En
def potentialb(site):
(x, y)=site.pos
da1=cos((-sqrt(3)*x/2+y/2)*ab)+cos(-y*ab)+\
cos((sqrt(3)*x/2+y/2)*ab)
da2=sin((-sqrt(3)*x/2+y/2)*ab)+sin(-y*ab)+\
sin((sqrt(3)*x/2+y/2)*ab)
return (2*pi/gamma)*(da1+sqrt(3)*da2)*En
def onsitea(site,p):
return potentiala(site)-p.mu
def onsiteb(site,p):
return potentialb(site)-p.mu
sys = kwant.Builder()
sys[a.shape(slab_shape,(1,1))]=onsitea
sys[b.shape(slab_shape,(1,1))]=onsiteb
# specify the hoppings of the graphene lattice in the
# fomat expected by builder.HoppingKind
def Ax(x,y):
result=-En*2*(pi/gamma)*(cos((-sqrt(3)*x/2+y/2)*ab)+\
cos((sqrt(3)*x/2+y/2)*ab)-2*cos(-y*ab))
return result
def Ay(x,y):
result=-En*2*sqrt(3)*(pi/gamma)*(cos((-sqrt(3)*x/2+y/2)*ab)\
-cos((sqrt(3)*x/2+y/2)*ab))
return result
def hopping1(site1,site2,p):
(x1,y1)=site1.pos
result=(-p.t+2*Ax(x1,y1)/3)\
*exp(1j * 2 * pi * 0.4 * x1/(gamma * gamma))
return result
def hopping2(site1,site2,p):
(x1,y1)=site2.pos
result=-p.t+Ay(x1-0.25,y1-sqrt(3)/4)/sqrt(3)\
-Ax(x1-0.25,y1-sqrt(3)/4)/3
return result
def hopping3(site1,site2,p):
(x1,y1)=site2.pos
result=-p.t-Ay(x1+0.25,y1-sqrt(3)/4)/sqrt(3)\
-Ax(x1+0.25,y1-sqrt(3)/4)/3
return result
kind1=kwant.builder.HoppingKind((0,0),a,b)
sys[kind1]=hopping1
kind2=kwant.builder.HoppingKind((0,1),a,b)
sys[kind2]=hopping2
kind3=kwant.builder.HoppingKind((-1,1),a,b)
sys[kind3]=hopping3
#### Define the leads. ####
# left lead0
sym0 = kwant.TranslationalSymmetry(graphene.vec((1, 0)))
def lead0_shape(pos):
x, y = pos
return (W/2-2.5*gamma <= y <= W/2+2.5*gamma)
def lead_hopping(site1,site2,p):
return -p.t
def lead_onsite(site1,p):
return mu_lead
lead0 = kwant.Builder(sym0)
hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
lead0[graphene.shape(lead0_shape, (0, W/2))] = lead_onsite
lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] =lead_hopping
#sys.attach_lead(lead0)
# down lead1
sym1 = kwant.TranslationalSymmetry(graphene.vec((1, -2)))
def lead1_shape(pos):
x, y= pos
return (L/4-2.5*gamma <= x <= L/4 + 2.5*gamma )
lead1 =kwant.Builder(sym1)
lead1[graphene.shape(lead1_shape, (L/4-2.5*gamma+1,0))] = lead_onsite
lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = lead_hopping
#sys.attach_lead(lead1)
# right lead2
lead2=lead0.reversed()
#sys.attach_lead(lead2)
# right-up lead3
sym3= kwant.TranslationalSymmetry(graphene.vec((-1,2)))
def lead3_shape(pos):
x, y = pos
return (0.75*L-2.5*gamma <= x <= 0.75*L +2.5*gamma)
lead3=kwant.Builder(sym3)
lead3[graphene.shape(lead3_shape,(0.75*L-2.5*gamma+1,0))] = lead_onsite
lead3[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]]= lead_hopping
#sys.attach_lead(lead3)
# left-up lead4
lead4=lead1.reversed()
#sys.attach_lead(lead4)
#print("ab",ab)
return sys, [lead0,lead1,lead2,lead3,lead4]
def calculate_sigmas(G):
# reduce by one dimension G -> G[temp, temp]
temp = [0, 1, 3, 4]
G = G[temp, :]
G = G[:, temp]
# invert R = G^-1
# find out whether it is a numpy object
r = np.linalg.inv(G)
# Voltages follow: V = R I[temp]
V = r.dot(np.array([1, -1, 0, 0]))
# Completely solved ihe six terminal system.
# Consider the 2x2 conductance now: Use I = sigma U
#E_x = V[1] - V[0]
E_x=(V[3]-V[2])
#print("V=",V)
#print("V[0]=",V[0],"V[1]=",V[1],"V[2]=",V[2],"V[3]=",V[3])
#E_x=(V[3]-V[2])/L
#E_y = V[1] - V[3]
E_y=V[3]
#E_y=V[3]/W
# formula above
sigma_xx = E_x / (E_x**2 + E_y**2)
sigma_xy = E_y / (E_x**2 + E_y**2)
return sigma_xx, sigma_xy
def plot_Hallconductance(sys,p):
# Compute conductance
Bs = np.linspace(0.5, 1, 10)
num_leads=len(sys.leads)
def G(sys, p):
smatrix = kwant.smatrix(sys, args=[p])
print("trans",smatrix.transmission(0,1))
Bs = np.linspace(0.5, 1, 10)
G = [[smatrix.transmission(i, j) for i in range(num_leads)] for j in range(num_leads)]
print(G)
G -= np.diag(np.sum(G, axis=0))
return calculate_sigmas(G)
sigmasxx, sigmasxy = np.array([G(sys, p) for p.B in Bs]).T
pyplot.figure()
pyplot.plot(Bs,sigmasxy)
pyplot.xlabel("B [t]")
pyplot.ylabel("sigma_{xy} [e^2/h]")
pyplot.show()
def main():
t0=2/(sqrt(3)*ab)
p=SimpleNamespace(t=t0, mu=0.25, B=None)
sys, leads = make_system(L, W, mu_lead=2)
# To highlight the two sublattices of graphene, we plot one with
# a filled, and the other one with an open circle:
def family_colors(site):
return 0 if site.family == a else 1
# Then, plot the system with leads.
#kwant.plot(sys, site_color=family_colors, site_lw=0.1,
# lead_site_lw=0, colorbar=False)
# Attach the leads to the system.
for lead in leads:
sys.attach_lead(lead)
# Then, plot the system with leads.
#kwant.plot(sys, site_color=family_colors, site_lw=0.1,
# lead_site_lw=0, colorbar=False)
# Finized the system
sys=sys.finalized()
# Plot conductance
# Compute conductance
plot_Hallconductance(sys,p)
if __name__ == '__main__':
main()