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 <abbout.adel@gmail.com> 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 <yu.li@postgrad.manchester.ac.uk> 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