diff --git a/__pycache__/momentum.cpython-38.pyc b/__pycache__/momentum.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..f009b97efe8fcf169b2785aed2744b20f5f77b62
Binary files /dev/null and b/__pycache__/momentum.cpython-38.pyc differ
diff --git a/__pycache__/tSnS.cpython-38.pyc b/__pycache__/tSnS.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..45462ff99d850960856b55c316c47e6591ea6ea7
Binary files /dev/null and b/__pycache__/tSnS.cpython-38.pyc differ
diff --git a/bands.py b/bands.py
new file mode 100644
index 0000000000000000000000000000000000000000..da232d48cea3c167c70bb52fec5ddb2a19137b17
--- /dev/null
+++ b/bands.py
@@ -0,0 +1,73 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from momentum import Momentum
+from tSnS import tSnS
+
+
+def Bandstructure(model, k_points, bandset=2, **kwargs):
+
+ """
+ #Calculate & Plot Bandstructure along irreducible
+ #path with given number of k_points
+ """
+
+ #################################
+ #Calculate Bandstructure
+ #################################
+ #Create path along irreducible path in 1. BZ
+ #High symmetry points
+ S = 0.5*(model.G1+model.G2)
+ X = 0.5*model.G1
+ Y = 0.5*model.G2
+ Gamma = 0.0*model.G1
+
+ irr_path = [Gamma, Y, S, Gamma, X, S]
+ k_path, high_sym_index = BZ.k_path(k_points, irr_path = irr_path)
+
+ band_colors = kwargs.get('band_colors', ['k']*model.N)
+ high_sym_marker = [r'$\Gamma$', r'$Y$', r'$S$', r'$\Gamma$', r'$X$', r'$S$']
+
+ e = []
+ for kpoint in k_path:
+
+ H0 = model.BH_Hamiltonian(kpoint, bandset=bandset)
+ e_temp = np.linalg.eigvalsh(H0)
+ e+= [e_temp]
+ e = np.stack(e)
+
+
+ #################################
+ #Plot Bandstructure
+ #################################
+
+ fig = plt.figure(figsize = (3,2))
+ ax =fig.add_subplot(1,1,1)
+ ax.set_ylabel(r'$\epsilon_b (eV)$')
+ ax.set_ylim(e.min(), e.max())
+
+ #Plot bands
+ x = np.arange(0,k_points,1)
+
+ for i in range (len(e[0])):
+ ax.plot(x, e[:,i], color = band_colors[i], linewidth = 1.7)
+
+ #Indicate high symmetry points
+ for i in range(len(high_sym_index)):
+ ax.axvline(high_sym_index[i], color='k', alpha=0.9, linewidth = 0.5)
+
+ ax.set_xticks(high_sym_index)
+ ax.set_xticklabels(high_sym_marker)
+ ax.set_xlim(high_sym_index[0], high_sym_index[-1])
+
+ plt.tight_layout()
+ plt.show()
+
+
+
+if __name__ == '__main__':
+
+ model = tSnS()
+ BZ = Momentum(model.G1, model.G2)
+ Bandstructure(model, 100, bandset=2)
+
diff --git a/momentum.py b/momentum.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe6e8f417b61bc7edf5fc859d8b6d13aa38e4fea
--- /dev/null
+++ b/momentum.py
@@ -0,0 +1,113 @@
+import numpy as np
+
+
+class Momentum:
+
+ """
+ General Momentum Mesh Class to create equidistant mesh
+ and irreducible path.
+ """
+
+ def __init__(self, G1, G2):
+
+ """
+ Parameters
+ ----------
+ Arguments: G1, G2 - Reciprocal lattice vectors
+ """
+
+ #Reciprocal lattice vectors
+ self.G1, self.G2 = G1, G2
+
+
+ def Bravais_mesh_simple(self, Nx, Ny=None):
+
+ """
+ Create simple Bravais mesh with Nx point along self.G1
+ and Ny points along self.G2
+ """
+ if Ny is None: Ny = Nx
+ mesh = np.zeros((Nx,Ny,3), dtype=float)
+
+ for i in range(Nx):
+ for j in range(Ny):
+ mesh[i,j] = i/Nx*self.G1 + j/Ny*self.G2
+
+ mesh = mesh.reshape(Nx*Ny,3)
+
+ return mesh
+
+
+ def k_path(self, k_points, irr_path):
+
+ """
+ Set up irreducible path along points in irr_path with given number
+ of k_points. Algorithm sets number of kpoints accoridng to actual
+ distance of high-symmetry point in momentum space.
+
+ Arguments: k_points - int, Number of kpoints along irr. path
+ irr_path - List of 3d Arrays with coordinates of high symmetry points
+
+
+ Return: k_path - List of 3d Arrays with momentum points on irr. path
+ marker - List of Int indicating the position of high symmetry points
+ """
+
+ #Compute distances between high symmetry points
+ high_sym = irr_path
+ distances = np.zeros(len(high_sym)-1, dtype = float)
+ for i in range(len(distances)):
+ distances[i] = np.linalg.norm(high_sym[i]-high_sym[i+1])
+ d_tot = np.sum(distances)
+
+ #Distribute points uniformly over the intervals
+ number_points=np.zeros_like(distances, dtype = np.int)
+ ratios = distances/d_tot
+ number_points = np.array([int(x*k_points) for x in ratios])
+ total_number_points = np.sum(number_points)
+
+ while (total_number_points < k_points):
+ number_points[0]+=1
+ total_number_points= np.sum(number_points)
+
+ #Create irr. path
+ kx, ky, kz = [],[],[]
+
+ #Create k_path such that the each segment contains the high symmetry point at
+ #the lower boundary and the last segment contains both high symmetry points
+ for i in range(len(distances)):
+
+ if i == len(distances)-1:
+
+ kx = np.append(kx, np.linspace(high_sym[i][0], high_sym[i+1][0], number_points[i], endpoint= True))
+ ky = np.append(ky, np.linspace(high_sym[i][1], high_sym[i+1][1], number_points[i], endpoint= True))
+ kz = np.append(kz, np.linspace(high_sym[i][2], high_sym[i+1][2], number_points[i], endpoint= True))
+
+ else:
+
+ kx = np.append(kx, np.linspace(high_sym[i][0], high_sym[i+1][0], number_points[i], endpoint= False))
+ ky = np.append(ky, np.linspace(high_sym[i][1], high_sym[i+1][1], number_points[i], endpoint= False))
+ kz = np.append(kz, np.linspace(high_sym[i][2], high_sym[i+1][2], number_points[i], endpoint= False))
+
+ k_path = []
+
+ for i in range(k_points):
+ k_path += [np.array([kx[i], ky[i], kz[i]])]
+
+
+ marker = [0]
+ for i in range(len(distances)):
+
+ if i == len(distances)-1:
+
+ new_marker = marker[-1]+number_points[i] -1
+ marker += [new_marker]
+
+ else:
+
+ new_marker = marker[-1]+number_points[i]
+ marker += [new_marker]
+
+ marker = np.asarray(marker)
+
+ return k_path, marker
diff --git a/tSnS.py b/tSnS.py
new file mode 100644
index 0000000000000000000000000000000000000000..7bfe6bee2a70bedac23f4f0e829702eee6fe48b9
--- /dev/null
+++ b/tSnS.py
@@ -0,0 +1,208 @@
+import sys
+import numpy as np
+
+
+class tSnS:
+
+
+ """
+ Sets up the geometric properties of twisted SnS and generate Bloch Hamiltonian.
+ """
+
+
+ def __init__(self, **kwargs):
+
+ #Initialize supercell vectors
+ self.A1 = np.array([40.718937498,4.292152880,0])
+ self.A2 = np.array([-4.071893750,38.629375919,0])
+ self.A3 = np.array([0,0,29.021099091])
+
+ #Basis transformation matrix (two-dimensional)
+ A_in_x = np.array([self.A1[:-1], self.A2[:-1]])
+ x_in_A = np.linalg.inv(A_in_x)
+ self.T = np.transpose(x_in_A)
+
+ #Reciprocal lattice vectors
+ self.volume = np.linalg.norm(np.cross(self.A1, self.A2))
+ self.G1 = 2*np.pi/self.volume*np.cross(self.A2, np.array([0.,0.,1.]))
+ self.G2 = 2*np.pi/self.volume*np.cross(np.array([0.,0.,1.]), self.A1)
+ self.G3 = -(self.G1+self.G2)
+ self.volume_BZ = np.linalg.norm(np.cross(self.G1, self.G2))
+
+ #Number of effective basis sites (sublattice + spin degrees of freedom)
+ self.N = 4
+
+
+ def phi_lm(self, momentum, l, m):
+
+ """
+ Phase factors of effective TB model
+ """
+
+ lattice_vec = l*self.A1 + m*self.A2
+ return np.exp(-1.j*np.dot(momentum, lattice_vec))
+
+
+ def BH_Hamiltonian(self, k, bandset=2):
+
+ """
+ Set up Bloch Hamiltonian at momentum point k.
+
+ ---------
+ Arguments: k - 3d Array with coordinates of momentum point
+ bandset (optional) - Bandset I/II
+ """
+ assert bandset in [1,2]
+
+ pauli = np.array([[1,0,0,1],[0,1,1,0], [0,-1j,1j,0], [1,0,0,-1]],
+ dtype=complex).reshape(4,2,2)
+
+ ham = np.zeros((2,2,2,2), dtype=complex)
+
+ if bandset == 1:
+
+ popt = [-0.28147216773927675, 0.003966634424731179,
+ 0.0029073602294481787, -4.4983518626626615e-07,
+ -0.004052637332926278, 0.00016072142349528914,
+ -3.7721284264216035e-05, -0.0009141524777857822,
+ -0.0011772772860306072, -0.0014951352789262786,
+ 0.0009438170971867384, -9.088937950755711e-05,
+ -0.000692733427968305, 0.0002266622565665427,
+ 0.0005191476715863555]
+
+ else:
+
+ popt = [-0.32047590155917033, 0.0021542300052420408,
+ 0.001469771721045236, 9.574936436961549e-05,
+ 0.005574114270507838, -0.00026463100214217314,
+ 0.0004064637000665132, 0.0004858014275756446,
+ 0.0006724769177247459, -0.00029328415914492434,
+ -0.0002932514194152098, 1.0106788563322311e-05,
+ 0.0007016914741439542, 0.0021456042488447995,
+ 0.0018535546929481545]
+
+
+
+
+ t_AA_0 = popt[0]
+ t_AB_0, t_AB_1 = popt[1], popt[2]
+ t_AA_x, t_AA_y = popt[3], popt[4]
+ t_AB_mx, t_AB_px = popt[5], popt[6]
+ t_AB_my, t_AB_py = popt[7], popt[8]
+ t_AA_m, t_AA_p = popt[9], popt[10]
+ t_AA_2x, t_AA_2y = popt[11], popt[12]
+ tR_AB_0, tR_AB_1 = popt[13], popt[14]
+
+
+ ##################################################
+ #Nearest-Neighbor/On-Site Coupling (kinetic terms)
+ ##################################################
+
+ #Onsite
+ ham[:,0,:,0] += t_AA_0*pauli[0]
+ ham[:,1,:,1] += t_AA_0*pauli[0]
+
+ #Nearest-Neighbor AB/BA coupling (1. nn)
+ hop = t_AB_0*(1.+self.phi_lm(k,0,-1))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = t_AB_1*(self.phi_lm(k,1,0)+self.phi_lm(k,1,-1))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ #Nearest-Neighbor AA/BB coupling (2. nn)
+ hop = t_AA_x*(self.phi_lm(k,1,0))*pauli[0]
+ ham[:,0,:,0] += hop
+ ham[:,0,:,0] += np.conjugate(hop).T
+ ham[:,1,:,1] += hop
+ ham[:,1,:,1] += np.conjugate(hop).T
+
+ hop = t_AA_y*(self.phi_lm(k,0,1))*pauli[0]
+ ham[:,0,:,0] += hop
+ ham[:,0,:,0] += np.conjugate(hop).T
+ ham[:,1,:,1] += hop
+ ham[:,1,:,1] += np.conjugate(hop).T
+
+
+ #################################################
+ #Higher Neighbor Couplings (kinetic terms)
+ #################################################
+
+ ###############
+ #AA/BB coupling
+ ################
+ #Corner states of first shell (3. nn)
+ hop = t_AA_m*self.phi_lm(k,1,1)
+ ham[:,1,:,1] += hop*pauli[0]
+ ham[:,1,:,1] += np.conjugate(hop)*pauli[0]
+
+ hop = t_AA_m*self.phi_lm(k,-1,1)
+ ham[:,0,:,0] += hop*pauli[0]
+ ham[:,0,:,0] += np.conjugate(hop)*pauli[0]
+
+ hop = t_AA_p*self.phi_lm(k,1,1)
+ ham[:,0,:,0] += hop*pauli[0]
+ ham[:,0,:,0] += np.conjugate(hop)*pauli[0]
+
+ hop = t_AA_p*self.phi_lm(k,-1,1)
+ ham[:,1,:,1] += hop*pauli[0]
+ ham[:,1,:,1] += np.conjugate(hop)*pauli[0]
+
+ #Second shell (5. nn)
+ hop = t_AA_2x*self.phi_lm(k,2,0)*pauli[0]
+ ham[:,0,:,0] += hop
+ ham[:,0,:,0] += np.conjugate(hop).T
+ ham[:,1,:,1] += hop
+ ham[:,1,:,1] += np.conjugate(hop).T
+
+ hop = t_AA_2y*self.phi_lm(k,0,2)*pauli[0]
+ ham[:,0,:,0] += hop
+ ham[:,0,:,0] += np.conjugate(hop).T
+ ham[:,1,:,1] += hop
+ ham[:,1,:,1] += np.conjugate(hop).T
+
+ ###############
+ #AB coupling
+ ################
+ #Edge states of second shell (4. nn)
+ hop = t_AB_mx*(self.phi_lm(k,2,0)+self.phi_lm(k,2,-1))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = t_AB_my*(self.phi_lm(k,1,1)+self.phi_lm(k,1,-2))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = t_AB_px*(self.phi_lm(k,-1,0)+self.phi_lm(k,-1,-1))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = t_AB_py*(self.phi_lm(k,0,1)+self.phi_lm(k,0,-2))*pauli[0]
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+
+ ######################################
+ #SOC Rashba-Coupling
+ ######################################
+ #Nearest-neighbors (1. nn)
+ hop = 1j*tR_AB_0*(-0.5*pauli[1]-0.5*pauli[2])
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = 1j*tR_AB_0*(-0.5*pauli[1]+0.5*pauli[2])*self.phi_lm(k,0,-1)
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = 1j*tR_AB_1*(0.5*pauli[1]-0.5*pauli[2])*self.phi_lm(k,1,0)
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ hop = 1j*tR_AB_1*(0.5*pauli[1]+0.5*pauli[2])*self.phi_lm(k,1,-1)
+ ham[:,0,:,1] += hop
+ ham[:,1,:,0] += np.conjugate(hop).T
+
+ ham = ham.reshape((4,4))
+
+ return ham