Dear Yu Li,
It seems that you deal within Slater Koster framework. I would like to
share with you a enclosed example (I usual give to my student) which
explains how to get the band structure of a real case .
Let as assume that we have a flat 2d shape of silicene (buckling is not
considered for simplicity) and we consider that we know the : Vppσ, as well
as Vppπ. To sum up let us say (we have Es, px, py and pz like-orbitals)
For the consideration of Fm and NM we have to generalize this example to
your case.
Hop this will help
best, Adel
## ------ IMPORTING PACKAGES -------------------------
import numpy as np
import tinyarray
from numpy import pi,sqrt,arccos
import kwant
from kwant.wraparound import wraparound,plot_2d_bands
from matplotlib import pyplot
import scipy
###
-----------------------------------------------------------------------------------
### -------------------------- THE MAIN USEFUL
FUNCTIONS-------------------------------
###
-----------------------------------------------------------------------------------
def wraparound_syst(syst):
return kwant.wraparound.wraparound(syst).finalized() ## or
wraparound(syst).finalized()
def get_A(syst):
B = np.asarray(syst._wrapped_symmetry.periods).T
#print("============ np.linalg.pinv(B).T =============")
#print(np.linalg.pinv(B).T)
return np.linalg.pinv(B).T
def momentum_to_lattice(k, syst):
"""Transform momentum to the basis of reciprocal lattice vectors.
See
https://en.wikipedia.org/wiki/Reciprocal_lattice#Generalization_of_a_dual_lattice
"""
A = get_A(syst)
#print("============ A =============")
#print(A)
#print(scipy.linalg.lstsq(A, k))
k, residuals = scipy.linalg.lstsq(A, k)[:2]
if np.any(abs(residuals) > 1e-7):
raise RuntimeError("Requested momentum doesn't correspond"
" to any lattice momentum.")
#print("============ k =============")
#print(k)
return k
### for 2D, we omit k_z: k =(kx, ky). for 3D, we set k =(kx, ky, kz)
def energies(k, syst):
"""
THE SOLUTION OF THE EIGENVALUE PROBLEM
"""
## ----------- FOR 2D symmetry --------------
k_x, k_y = momentum_to_lattice(k, syst)
params = {'k_x': k_x, 'k_y': k_y}
## ----------- FOR 3D symmetry --------------
## we use
#k_x, k_y, k_z = momentum_to_lattice(k, syst)
#params = {'k_x': k_x, 'k_y': k_y, 'k_z': k_z}
H = syst.hamiltonian_submatrix(params=params)
energies = np.sort(np.linalg.eigvalsh(H))
return energies ## or np.array(energies)
###
-----------------------------------------------------------------------------------
### ------- THE TIGHT-BINDING PARAMETERS WITHIN SLATER FRAMEWORK
----------------------
###
-----------------------------------------------------------------------------------
## --------------------- For example silicene (buckling is not considered
for simplicity reason)
Es_a = -4.0497;
Es_b = -4.0497;
## ---------------------
Ep_a = 1.0297;
Ep_b = 1.0297;
## ---------------------
Vss_sig_aa = -0.0000
Vss_sig_bb = -0.0000
Vss_sig_ab = -2.0662
## ---------------------
Vsp_sig_aa = +0.0000
Vsp_sig_bb = +0.0000
Vsp_sig_ab = +2.0850
## --------------------
Vpp_sig_aa = +0.8900
Vpp_sig_bb = +0.8900
Vpp_sig_ab = +3.1837
## --------------------
Vpp_pi_aa = -0.3612
Vpp_pi_bb = -0.3612
Vpp_pi_ab = -0.9488
## ---------------------------------------------
## ====== Onsites for the system ===============
def Onsite_Mo(site):
"""
This function is used when we take into account masse term wher Ep_a =
- Ep_b
"""
if site.family == a :
## -----------------------------------------------------------------
Epx, Epy, Epz = Ep_a*np.ones(3)
Es = Es_a
Onsite_Mat = tinyarray.array([ [ Es, 0, 0, 0],
[ 0, Epx, 0, 0], ## Epx =
Ep
[ 0, 0, Epy, 0], ## Epy =
Ep
[ 0, 0, 0, Epz], ## Epz =
Ep
])
else:
##
------------------------------------------------------------------
## If we include the masse terme we set Epx, Epy, Epz =
-1*Ep_b*np.ones(3)
Epx, Epy, Epz = Ep_b*np.ones(3)
Es = Es_b
Onsite_Mat = tinyarray.array([ [ Es, 0, 0, 0],
[ 0, Epx, 0, 0], ## Epx =
Ep
[ 0, 0, Epy, 0], ## Epy =
Ep
[ 0, 0, 0, Epz], ## Epz =
Ep
])
return Onsite_Mat
## ------------------------------------------------
## ====== Hopping for the system ===============
def Hopp_Mo(site_1, site_2):
## ----------------------------------------------------
"""
the hopping is between A-A B-B and A-B and given by slater koster
parameters
"""
x_1, y_1 = site_1.pos[0], site_1.pos[1]
x_2, y_2 = site_2.pos[0], site_2.pos[1]
Norm_Vec = np.sqrt( (x_2 - x_1)**2 + (y_2 - y_1)**2)
## in cease of 3d we set
#x_1, y_1, z_1 = site_1.pos[0], site_1.pos[1], site_1.pos[2]
#x_2, y_2, z_2 = site_2.pos[0], site_2.pos[1], site_2.pos[2]
#Norm_Vec = np.sqrt( (x_2 - x_1)**2 + (y_2 - y_1)**2+ (z_2 - z_1)**2)
## ---------------------
## Cosinus Directions
l = (x_2 - x_1)/Norm_Vec
m = (y_2 - y_1)/Norm_Vec
## For 3D we set : n = (z_2 - z_1)/Norm_Vec
## In the actual case wedeal with 2d so n=0
n = 0
## ----------------------------------------------------
if (site_1.family == a and site_2.family == b):
## ---------------------------------------------------
#print("===============")
#print(l)
#print(m)
#print(n)
S_S = Vss_sig_ab
S_px = l*Vsp_sig_ab
S_py = m*Vsp_sig_ab
S_pz = n*Vsp_sig_ab
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_ab + (1-l**2)*Vpp_pi_ab
px_py = (l*m)*(Vpp_sig_ab - Vpp_pi_ab)
px_pz = (l*n)*(Vpp_sig_ab - Vpp_pi_ab)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_ab + (1-m**2)*Vpp_pi_ab
py_pz = (m*n)*(Vpp_sig_ab - Vpp_pi_ab)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_ab + (1-n**2)*Vpp_pi_ab
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ##
-pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ##
pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ##
pza_pxb = pxa_pzb
])
elif (site_1.family == a and site_2.family == a):
## ---------------------------------------------------
S_S = Vss_sig_aa
S_px = l*Vsp_sig_aa
S_py = m*Vsp_sig_aa
S_pz = n*Vsp_sig_aa
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_aa + (1-l**2)*Vpp_pi_aa
px_py = (l*m)*(Vpp_sig_aa - Vpp_pi_aa)
px_pz = (l*n)*(Vpp_sig_aa - Vpp_pi_aa)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_aa + (1-m**2)*Vpp_pi_aa
py_pz = (m*n)*(Vpp_sig_aa - Vpp_pi_aa)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_aa + (1-n**2)*Vpp_pi_aa
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ##
-pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ##
pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ##
pza_pxb = pxa_pzb
])
elif (site_1.family == b and site_2.family == b):
## ----------------------------------------------------
S_S = Vss_sig_bb
S_px = l*Vsp_sig_bb
S_py = m*Vsp_sig_bb
S_pz = n*Vsp_sig_bb
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_bb + (1-l**2)*Vpp_pi_bb
px_py = (l*m)*(Vpp_sig_bb - Vpp_pi_bb)
px_pz = (l*n)*(Vpp_sig_bb - Vpp_pi_bb)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_bb + (1-m**2)*Vpp_pi_bb
py_pz = (m*n)*(Vpp_sig_bb - Vpp_pi_bb)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_bb + (1-n**2)*Vpp_pi_bb
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ##
-pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ##
pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ##
pza_pxb = pxa_pzb
])
return Hop_Mat
## -------------------------------- Make a lattice
-------------------------------------
##
=====================================================================================
def make_2dHonycomb():
a_sisi = 1.42;
"""Create a builder for a periodic infinite sheet of graphene."""
a1 = a_sisi*sqrt(3)
lat = kwant.lattice.general([[a1, 0], [a1/2, a1*sqrt(3)/2]], #
lattice vectors
[[ 0, 0], [ 0, a1/sqrt(3)]], #
Positions
name = ['a', 'b'])
# sublatice'names
a, b = lat.sublattices
## ---------- For unit cell
----------------------------------------------------------------
sym = kwant.TranslationalSymmetry(lat.vec((1, 0)), lat.vec((0, 1)))
sys = kwant.Builder(sym)
sys[a.shape(lambda p: True, (0, 0))] = Onsite_Mo
sys[b.shape(lambda p: True, (0, 0))] = Onsite_Mo
sys[a.neighbors(1)] = Hopp_Mo
sys[b.neighbors(1)] = Hopp_Mo
hoppings_ab = (((0, 0), a, b), ((0, +1), a, b), ((-1, +1), a, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings_ab]] =
Hopp_Mo
#sys[lat.neighbors(2)] =
##
-------------------------------------------------------------------------------------------
return sys, (a, b)
### ------------------- DEFINING THE UNIT CELL --------------
### =========================================================
def hop_lw(site1, site2):
return 0.01 if site1.family == site2.family else 0.03
def hopping_color(site1,site2):
return 'r' if site1.family==site2.family else 'k'
syst, (a, b) = make_2dHonycomb()
kwant.plot(syst, site_size=0.1, site_lw=0.012, hop_color=hopping_color,
hop_lw=hop_lw)
pyplot.show()
### ===============================
unit_cell = wraparound_syst(syst)
## =================================== High Symmetry
Points=========================================
## ----- G ------- -------- K -------- ------ M -----
------- G -------
## x- 0 2*pi/(3*sqrt(3)*a_cc) 0
0
## y- 0 2*pi/(3*a_cc) 2*pi/(3*a_cc)
0
##
=================================================================================================
a_sisi = 1.42;
a_cc = a_sisi
ac = sqrt(3)*a_cc
d_GK = sqrt((2*pi/(3*sqrt(3)*a_cc)-(0))**2 + ( 2*pi/(3*a_cc)-0)**2)
d_KM = sqrt(((0)-2*pi/(3*sqrt(3)*a_cc))**2 + (
2*pi/(3*a_cc)-2*pi/(3*a_cc))**2)
d_MG = sqrt(((0)-(0))**2 + ((0) - 2*pi/(3*a_cc))**2)
print("the distance between G-K is " , d_GK )
print("the distance between K-M is " , d_KM )
print("the distance between M-G is ", d_MG )
## -------------------------------------------------
## we set the nb for np.linspace(, ,nb) for each path
number = 400 ## we set this for the long path which is 1.7 for K1-M
nb_GK = number ### number*(d_GK/d_GK)= number## since number sinse GK is
the longest path
nb_KM = number*(d_KM/d_GK)
nb_MG = number*(d_MG/d_GK)
## ---------------------------------------------------------------------
kx_KG = np.linspace( 2*pi/(3*sqrt(3)*a_cc), 0,
int(nb_GK) )
ky_KG = np.linspace( 2*pi/(3*a_cc), 0,
int(nb_GK) )
## ---------------------------------------------------------------------
kx_GM = np.linspace( 0, 2*pi/(3*a_cc),
int(nb_MG) )
ky_GM = np.linspace( 0, 2*pi/(3*a_cc),
int(nb_MG) )
## ---------------------------------------------------------------------
kx_MK = np.linspace( 2*pi/(3*a_cc), 2*pi/(3*sqrt(3)*a_cc),
int(nb_KM) )
ky_MK = 2*pi/(3*a_cc)*np.ones(int(nb_KM))
## ---------------------------------------------------------------------
## For Kx
Kx_path=[]
Kx_path = np.append(Kx_path, kx_KG)
Kx_path = np.append(Kx_path, kx_GM)
Kx_path = np.append(Kx_path, kx_MK)
## For Kx
Ky_path=[]
Ky_path = np.append(Ky_path, ky_KG)
Ky_path = np.append(Ky_path, ky_GM)
Ky_path = np.append(Ky_path, ky_MK)
##
--------------------------------------------------------------------------------
## Extracting the energy bands along high symmetry point
==========================
Energy_bands = [energies((kx, ky), unit_cell) for (kx, ky) in zip(Kx_path,
Ky_path)]
## for HSP plot, we can just set a random vectors wich has the size of
Energy_bands
## Here will will show G-K-M-K so we do not need the exact values of
reciprocal vector
size_HSP = np.size(np.array(Energy_bands)[:,0])
HSP_For_Plot = np.linspace(0, 1, size_HSP)
###
---------------------------------------------------------------------------------------------
###====================== SETTING THE PLOT PARAMETERS FOR G-K-M-G
===============================
### -----------------------
----------------------------------------------------------------------
fig = pyplot.figure()
#axes.axis('equal')
axes = fig.add_subplot(1, 1, 1)
for i in range(np.size(np.array(Energy_bands)[0])):
axes.plot(HSP_For_Plot, np.array(Energy_bands)[:,i])
## -------------------------------------------------------
## Setting the Horizontal lines for the high symmetry path
axes.axvline(x = HSP_For_Plot[0], linewidth=2, color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)], linewidth=2, color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)+ np.size(kx_GM)], linewidth=2,
color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)+ np.size(kx_GM)+
np.size(kx_MK)-1], linewidth=2, color='k')
axes.axhline(y = 0, color='k')
##
----------------------------------------------------------------------------
## for HSP plot, we can just set a random vecto wich has the size of
Energy_bands
## the normalized position for each high symmetry ponits are defined as:
dK = HSP_For_Plot [0]
dKG = HSP_For_Plot[np.size(kx_KG)]
dKGM = HSP_For_Plot[np.size(kx_KG)+np.size(kx_GM)]
dKGMK = HSP_For_Plot[np.size(kx_KG)+np.size(kx_GM)+np.size(kx_MK)-1]
HighSymmetry_positions = [dK , dKG, dKGM, dKGMK]
HighSymmetry_names= ('$K$',r'$\Gamma$', '$M$', '$K$')
## calling xticks and plotting horizontal lines related to G, K, M, and G
from matplotlib.pyplot import xticks
locs, labels = xticks()
axes.xaxis.set_visible(True)
axes.yaxis.set_visible(True)
xticks(HighSymmetry_positions, HighSymmetry_names)
Le mer. 5 févr. 2020 à 19:56, Abbout Adel <[email protected]> a écrit :
> Dear Yu,
>
> In kwant you simulate models. So it is up to you to decide what
> corresponds best to your real system.
>
> Yes you can change the value of the hopping depending on the justification
> you may find. It is also possible to put it as a parameter and do your
> study as a function of the value of "t". You can also choose it as a matrix
> in order to take into account the different orbitals in your system
>
> I hope this helps,
> Adel
>
> On Wed, Feb 5, 2020 at 11:22 AM Yu Li <[email protected]>
> wrote:
>
>> Dear all,
>>
>> Now I’m studying some transport property based on multilayered system,
>> with nonmagnetic layer/ferromagnetic layer (NM/FM) stack.
>>
>> May I ask when creating a tight binding system, and setting hoppings by
>> the following function:
>>
>> syst[(lat(1, 0), lat(1, 1))] = t,
>>
>> how to associate the hopping with a realistic system? Such as hoppings
>> between FM and FM atoms, or between FM and NM atoms? I found in many
>> literatures people just simply put t=1, but I assume the value should
>> depend on the type of atoms.
>>
>> My another question is, is it possible to define more than one type of
>> hopping between two sites?
>> For the tight-binding approximation of p electrons, there are two types
>> of hopping integrals: Vppσ, and Vppπ, so I’m wondering if it’s possible to
>> simulate the effect of both hoppings?
>>
>> Thanks a lot for reading my questions!
>> Cheers,
>> Yu Li
>>
>
>
> --
> Abbout Adel
>
Le mer. 5 févr. 2020 à 09:22, Yu Li <[email protected]> a
écrit :
> Dear all,
>
> Now I’m studying some transport property based on multilayered system,
> with nonmagnetic layer/ferromagnetic layer (NM/FM) stack.
>
> May I ask when creating a tight binding system, and setting hoppings by
> the following function:
>
> syst[(lat(1, 0), lat(1, 1))] = t,
>
> how to associate the hopping with a realistic system? Such as hoppings
> between FM and FM atoms, or between FM and NM atoms? I found in many
> literatures people just simply put t=1, but I assume the value should
> depend on the type of atoms.
>
> My another question is, is it possible to define more than one type of
> hopping between two sites?
> For the tight-binding approximation of p electrons, there are two types of
> hopping integrals: Vppσ, and Vppπ, so I’m wondering if it’s possible to
> simulate the effect of both hoppings?
>
> Thanks a lot for reading my questions!
> Cheers,
> Yu Li
>
## ------ IMPORTING PACKAGES -------------------------
import numpy as np
import tinyarray
from numpy import pi,sqrt,arccos
import kwant
from kwant.wraparound import wraparound,plot_2d_bands
from matplotlib import pyplot
import scipy
### -----------------------------------------------------------------------------------
### -------------------------- THE MAIN USEFUL FUNCTIONS-------------------------------
### -----------------------------------------------------------------------------------
def wraparound_syst(syst):
return kwant.wraparound.wraparound(syst).finalized() ## or wraparound(syst).finalized()
def get_A(syst):
B = np.asarray(syst._wrapped_symmetry.periods).T
#print("============ np.linalg.pinv(B).T =============")
#print(np.linalg.pinv(B).T)
return np.linalg.pinv(B).T
def momentum_to_lattice(k, syst):
"""Transform momentum to the basis of reciprocal lattice vectors.
See https://en.wikipedia.org/wiki/Reciprocal_lattice#Generalization_of_a_dual_lattice
"""
A = get_A(syst)
#print("============ A =============")
#print(A)
#print(scipy.linalg.lstsq(A, k))
k, residuals = scipy.linalg.lstsq(A, k)[:2]
if np.any(abs(residuals) > 1e-7):
raise RuntimeError("Requested momentum doesn't correspond"
" to any lattice momentum.")
#print("============ k =============")
#print(k)
return k
### for 2D, we omit k_z: k =(kx, ky). for 3D, we set k =(kx, ky, kz)
def energies(k, syst):
"""
THE SOLUTION OF THE EIGENVALUE PROBLEM
"""
## ----------- FOR 2D symmetry --------------
k_x, k_y = momentum_to_lattice(k, syst)
params = {'k_x': k_x, 'k_y': k_y}
## ----------- FOR 3D symmetry --------------
## we use
#k_x, k_y, k_z = momentum_to_lattice(k, syst)
#params = {'k_x': k_x, 'k_y': k_y, 'k_z': k_z}
H = syst.hamiltonian_submatrix(params=params)
energies = np.sort(np.linalg.eigvalsh(H))
return energies ## or np.array(energies)
### -----------------------------------------------------------------------------------
### ------- THE TIGHT-BINDING PARAMETERS WITHIN SLATER FRAMEWORK ----------------------
### -----------------------------------------------------------------------------------
## --------------------- For example silicene (buckling is not considered for simplicity reason)
Es_a = -4.0497;
Es_b = -4.0497;
## ---------------------
Ep_a = 1.0297;
Ep_b = 1.0297;
## ---------------------
Vss_sig_aa = -0.0000
Vss_sig_bb = -0.0000
Vss_sig_ab = -2.0662
## ---------------------
Vsp_sig_aa = +0.0000
Vsp_sig_bb = +0.0000
Vsp_sig_ab = +2.0850
## --------------------
Vpp_sig_aa = +0.8900
Vpp_sig_bb = +0.8900
Vpp_sig_ab = +3.1837
## --------------------
Vpp_pi_aa = -0.3612
Vpp_pi_bb = -0.3612
Vpp_pi_ab = -0.9488
## ---------------------------------------------
## ====== Onsites for the system ===============
def Onsite_Mo(site):
"""
This function is used when we take into account masse term wher Ep_a = - Ep_b
"""
if site.family == a :
## -----------------------------------------------------------------
Epx, Epy, Epz = Ep_a*np.ones(3)
Es = Es_a
Onsite_Mat = tinyarray.array([ [ Es, 0, 0, 0],
[ 0, Epx, 0, 0], ## Epx = Ep
[ 0, 0, Epy, 0], ## Epy = Ep
[ 0, 0, 0, Epz], ## Epz = Ep
])
else:
## ------------------------------------------------------------------
## If we include the masse terme we set Epx, Epy, Epz = -1*Ep_b*np.ones(3)
Epx, Epy, Epz = Ep_b*np.ones(3)
Es = Es_b
Onsite_Mat = tinyarray.array([ [ Es, 0, 0, 0],
[ 0, Epx, 0, 0], ## Epx = Ep
[ 0, 0, Epy, 0], ## Epy = Ep
[ 0, 0, 0, Epz], ## Epz = Ep
])
return Onsite_Mat
## ------------------------------------------------
## ====== Hopping for the system ===============
def Hopp_Mo(site_1, site_2):
## ----------------------------------------------------
"""
the hopping is between A-A B-B and A-B and given by slater koster parameters
"""
x_1, y_1 = site_1.pos[0], site_1.pos[1]
x_2, y_2 = site_2.pos[0], site_2.pos[1]
Norm_Vec = np.sqrt( (x_2 - x_1)**2 + (y_2 - y_1)**2)
## in cease of 3d we set
#x_1, y_1, z_1 = site_1.pos[0], site_1.pos[1], site_1.pos[2]
#x_2, y_2, z_2 = site_2.pos[0], site_2.pos[1], site_2.pos[2]
#Norm_Vec = np.sqrt( (x_2 - x_1)**2 + (y_2 - y_1)**2+ (z_2 - z_1)**2)
## ---------------------
## Cosinus Directions
l = (x_2 - x_1)/Norm_Vec
m = (y_2 - y_1)/Norm_Vec
## For 3D we set : n = (z_2 - z_1)/Norm_Vec
## In the actual case wedeal with 2d so n=0
n = 0
## ----------------------------------------------------
if (site_1.family == a and site_2.family == b):
## ---------------------------------------------------
#print("===============")
#print(l)
#print(m)
#print(n)
S_S = Vss_sig_ab
S_px = l*Vsp_sig_ab
S_py = m*Vsp_sig_ab
S_pz = n*Vsp_sig_ab
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_ab + (1-l**2)*Vpp_pi_ab
px_py = (l*m)*(Vpp_sig_ab - Vpp_pi_ab)
px_pz = (l*n)*(Vpp_sig_ab - Vpp_pi_ab)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_ab + (1-m**2)*Vpp_pi_ab
py_pz = (m*n)*(Vpp_sig_ab - Vpp_pi_ab)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_ab + (1-n**2)*Vpp_pi_ab
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ## -pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ## pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ## pza_pxb = pxa_pzb
])
elif (site_1.family == a and site_2.family == a):
## ---------------------------------------------------
S_S = Vss_sig_aa
S_px = l*Vsp_sig_aa
S_py = m*Vsp_sig_aa
S_pz = n*Vsp_sig_aa
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_aa + (1-l**2)*Vpp_pi_aa
px_py = (l*m)*(Vpp_sig_aa - Vpp_pi_aa)
px_pz = (l*n)*(Vpp_sig_aa - Vpp_pi_aa)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_aa + (1-m**2)*Vpp_pi_aa
py_pz = (m*n)*(Vpp_sig_aa - Vpp_pi_aa)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_aa + (1-n**2)*Vpp_pi_aa
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ## -pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ## pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ## pza_pxb = pxa_pzb
])
elif (site_1.family == b and site_2.family == b):
## ----------------------------------------------------
S_S = Vss_sig_bb
S_px = l*Vsp_sig_bb
S_py = m*Vsp_sig_bb
S_pz = n*Vsp_sig_bb
## ----------------------------------------------------
px_px = (l**2)*Vpp_sig_bb + (1-l**2)*Vpp_pi_bb
px_py = (l*m)*(Vpp_sig_bb - Vpp_pi_bb)
px_pz = (l*n)*(Vpp_sig_bb - Vpp_pi_bb)
## ----------------------------------------------------
py_py = (m**2)*Vpp_sig_bb + (1-m**2)*Vpp_pi_bb
py_pz = (m*n)*(Vpp_sig_bb - Vpp_pi_bb)
## ----------------------------------------------------
pz_pz = (n**2)*Vpp_sig_bb + (1-n**2)*Vpp_pi_bb
## ----------------------------------------------------
Hop_Mat = tinyarray.array([ [ S_S, S_px, S_py, S_pz],
[-S_px, px_px, px_py, px_pz], ## -pya_Sb = -Sa_pxb
[-S_py, px_py, py_py, py_pz], ## pya_pxb = pxa_pyb
[-S_pz, px_pz, py_pz, pz_pz], ## pza_pxb = pxa_pzb
])
return Hop_Mat
## -------------------------------- Make a lattice -------------------------------------
## =====================================================================================
def make_2dHonycomb():
a_sisi = 1.42;
"""Create a builder for a periodic infinite sheet of graphene."""
a1 = a_sisi*sqrt(3)
lat = kwant.lattice.general([[a1, 0], [a1/2, a1*sqrt(3)/2]], # lattice vectors
[[ 0, 0], [ 0, a1/sqrt(3)]], # Positions
name = ['a', 'b']) # sublatice'names
a, b = lat.sublattices
## ---------- For unit cell ----------------------------------------------------------------
sym = kwant.TranslationalSymmetry(lat.vec((1, 0)), lat.vec((0, 1)))
sys = kwant.Builder(sym)
sys[a.shape(lambda p: True, (0, 0))] = Onsite_Mo
sys[b.shape(lambda p: True, (0, 0))] = Onsite_Mo
sys[a.neighbors(1)] = Hopp_Mo
sys[b.neighbors(1)] = Hopp_Mo
hoppings_ab = (((0, 0), a, b), ((0, +1), a, b), ((-1, +1), a, b))
sys[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings_ab]] = Hopp_Mo
#sys[lat.neighbors(2)] =
## -------------------------------------------------------------------------------------------
return sys, (a, b)
### ------------------- DEFINING THE UNIT CELL --------------
### =========================================================
def hop_lw(site1, site2):
return 0.01 if site1.family == site2.family else 0.03
def hopping_color(site1,site2):
return 'r' if site1.family==site2.family else 'k'
syst, (a, b) = make_2dHonycomb()
kwant.plot(syst, site_size=0.1, site_lw=0.012, hop_color=hopping_color, hop_lw=hop_lw)
pyplot.show()
### ===============================
unit_cell = wraparound_syst(syst)
## =================================== High Symmetry Points=========================================
## ----- G ------- -------- K -------- ------ M ----- ------- G -------
## x- 0 2*pi/(3*sqrt(3)*a_cc) 0 0
## y- 0 2*pi/(3*a_cc) 2*pi/(3*a_cc) 0
## =================================================================================================
a_sisi = 1.42;
a_cc = a_sisi
ac = sqrt(3)*a_cc
d_GK = sqrt((2*pi/(3*sqrt(3)*a_cc)-(0))**2 + ( 2*pi/(3*a_cc)-0)**2)
d_KM = sqrt(((0)-2*pi/(3*sqrt(3)*a_cc))**2 + ( 2*pi/(3*a_cc)-2*pi/(3*a_cc))**2)
d_MG = sqrt(((0)-(0))**2 + ((0) - 2*pi/(3*a_cc))**2)
print("the distance between G-K is " , d_GK )
print("the distance between K-M is " , d_KM )
print("the distance between M-G is ", d_MG )
## -------------------------------------------------
## we set the nb for np.linspace(, ,nb) for each path
number = 400 ## we set this for the long path which is 1.7 for K1-M
nb_GK = number ### number*(d_GK/d_GK)= number## since number sinse GK is the longest path
nb_KM = number*(d_KM/d_GK)
nb_MG = number*(d_MG/d_GK)
## ---------------------------------------------------------------------
kx_KG = np.linspace( 2*pi/(3*sqrt(3)*a_cc), 0, int(nb_GK) )
ky_KG = np.linspace( 2*pi/(3*a_cc), 0, int(nb_GK) )
## ---------------------------------------------------------------------
kx_GM = np.linspace( 0, 2*pi/(3*a_cc), int(nb_MG) )
ky_GM = np.linspace( 0, 2*pi/(3*a_cc), int(nb_MG) )
## ---------------------------------------------------------------------
kx_MK = np.linspace( 2*pi/(3*a_cc), 2*pi/(3*sqrt(3)*a_cc), int(nb_KM) )
ky_MK = 2*pi/(3*a_cc)*np.ones(int(nb_KM))
## ---------------------------------------------------------------------
## For Kx
Kx_path=[]
Kx_path = np.append(Kx_path, kx_KG)
Kx_path = np.append(Kx_path, kx_GM)
Kx_path = np.append(Kx_path, kx_MK)
## For Kx
Ky_path=[]
Ky_path = np.append(Ky_path, ky_KG)
Ky_path = np.append(Ky_path, ky_GM)
Ky_path = np.append(Ky_path, ky_MK)
## --------------------------------------------------------------------------------
## Extracting the energy bands along high symmetry point ==========================
Energy_bands = [energies((kx, ky), unit_cell) for (kx, ky) in zip(Kx_path, Ky_path)]
## for HSP plot, we can just set a random vectors wich has the size of Energy_bands
## Here will will show G-K-M-K so we do not need the exact values of reciprocal vector
size_HSP = np.size(np.array(Energy_bands)[:,0])
HSP_For_Plot = np.linspace(0, 1, size_HSP)
### ---------------------------------------------------------------------------------------------
###====================== SETTING THE PLOT PARAMETERS FOR G-K-M-G ===============================
### ---------------------------------------------------------------------------------------------
fig = pyplot.figure()
#axes.axis('equal')
axes = fig.add_subplot(1, 1, 1)
for i in range(np.size(np.array(Energy_bands)[0])):
axes.plot(HSP_For_Plot, np.array(Energy_bands)[:,i])
## -------------------------------------------------------
## Setting the Horizontal lines for the high symmetry path
axes.axvline(x = HSP_For_Plot[0], linewidth=2, color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)], linewidth=2, color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)+ np.size(kx_GM)], linewidth=2, color='k')
axes.axvline(x = HSP_For_Plot[np.size(kx_KG)+ np.size(kx_GM)+ np.size(kx_MK)-1], linewidth=2, color='k')
axes.axhline(y = 0, color='k')
## ----------------------------------------------------------------------------
## for HSP plot, we can just set a random vecto wich has the size of Energy_bands
## the normalized position for each high symmetry ponits are defined as:
dK = HSP_For_Plot [0]
dKG = HSP_For_Plot[np.size(kx_KG)]
dKGM = HSP_For_Plot[np.size(kx_KG)+np.size(kx_GM)]
dKGMK = HSP_For_Plot[np.size(kx_KG)+np.size(kx_GM)+np.size(kx_MK)-1]
HighSymmetry_positions = [dK , dKG, dKGM, dKGMK]
HighSymmetry_names= ('$K$',r'$\Gamma$', '$M$', '$K$')
## calling xticks and plotting horizontal lines related to G, K, M, and G
from matplotlib.pyplot import xticks
locs, labels = xticks()
axes.xaxis.set_visible(True)
axes.yaxis.set_visible(True)
xticks(HighSymmetry_positions, HighSymmetry_names)
