Add Bandstructure Features

bands.py

+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)
+ 

momentum.py

+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

tSnS.py

+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