From ee4b2af1ed80cd833870c5c94ec9b82ed558a63e Mon Sep 17 00:00:00 2001
From: Ammon Fischer <ammon.fischer@rwth-aachen.de>
Date: Mon, 6 Jun 2022 19:08:27 -0400
Subject: [PATCH] Add Bandstructure Features
---
__pycache__/momentum.cpython-38.pyc | Bin 0 -> 2904 bytes
__pycache__/tSnS.cpython-38.pyc | Bin 0 -> 5153 bytes
bands.py | 73 ++++++++++
momentum.py | 113 +++++++++++++++
tSnS.py | 208 ++++++++++++++++++++++++++++
5 files changed, 394 insertions(+)
create mode 100644 __pycache__/momentum.cpython-38.pyc
create mode 100644 __pycache__/tSnS.cpython-38.pyc
create mode 100644 bands.py
create mode 100644 momentum.py
create mode 100644 tSnS.py
diff --git a/__pycache__/momentum.cpython-38.pyc b/__pycache__/momentum.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..f009b97efe8fcf169b2785aed2744b20f5f77b62
GIT binary patch
literal 2904
zcmb7GL66(U6`mQA62;ZZioNY_;vR@52&^EhjkiGy1Ph~%jV=N-Lg6$Q0R}^JxLR5i
zsht_wUV#D%tkL=sz`f)c#K+wG8+z`U*PIsUFX*NH9x1J*It9{M%<K&3y{B*9_vXF%
za(mljXuE&@_R;M&V?WShy;)fN7G1xCK{ClxHsr@l3TZuO(z*~sA^8`qXJ2BU4I&!X
z+m-$Y_!>i3GSV1YwS}QAThft@=X~hMrgWdPVME%m=auf~`6SNFY%+h~wT|N~R$=O|
zzWJZW`qaOdhFbe3_ahaDCidgASt1i{!p!&+67X7KCjCUISk9v4Q5yTx(46+YRz$c~
zjpo$vq3i1y3SJ1o3Kmf3)%nDH;I&$R2vs<VO{_E)`&SK%yJ|EexcU?S`1XhX@u9!(
zKZv7bs`3c9)6keCiv1^XWb#UT&=aLNzKwp^vr1dX>4_o`%Es_LQw{w4*YUxZ8b5Xr
zP9}MFlIZ9(RtNWEJvRCDz$Bwnv;Qc`WRi{cC%KH%ee)E@kLAA3^uc85zy9>Rbb}zt
z5)%YlFslUy)|$Vy@@9XvJJ|3GL7-mw#b8o^xv{I-1z!lm9&_@O6epHJG<+=<maz(c
z!SAz&zh2mwi)wyr;S_ekZ^D<4?HyJ)g{bDv4%=mgg+1p)9I>;H&u$ggcMY<+gZO#{
zdxb@jH6Oyi_iBc9GMR$KPgVFNOmt21)5M(mgJ=FUPr!Sa=Gn+6clM8OLy-bta9(ff
zjom}+4lc>(-}BOu=6o9Wnx)95r8AFJu1n`6%|lZ<D$GXlYjVn_LTrVQZYT%Zr5HRz
zKOfo@b0tqo{<xyA9-1U_V~IE*k)S@sE;e)*18ca<<~yRz+x(X3iuo_Di1`t<4T?ye
zE-;-!{~f#*L<6WXqOdLyx_ivvs(~Ea*!X5&)C9e^a3JNt-B#heO*A-?@3Obq9(%=!
z5)ksfVc#_i`xRd_AnO)3##UtuUyH^~ws0Zg{e_)yX-WIJTe$bx1&4gI(rjYZDq1ht
zn5n-NP3btWFZpi(_K3ZfID0u_^k*@_7T#iO(OyUTl0UqU0CH;Vbe^$z$<?2rzg2iZ
z)UHn;H%?KlnshSRkj<}c1QTSTv9+eb*Uv3!wW$GFXtdWf+zky{ZECd7{_&G>wf=8h
zD~+FyYo!5PBevK!_6zWX>q%E~8%VvvtJvSF#x{*3>lGh?qs8`1arT~cMtsp(=i}jn
zH`oyBo1!DVl{a?5@g7@TD>}us%2FqLXLEmj55MvikZcznL`$c@{}X{rzZL$IBeGkq
zFJj}*rf;P1<wdTpy!Mlf1Om$Kew7VI$&)zqv)SZPtY9QR@y9_e+xPFLqg)|nOpqju
zUJCgnva}XP5d<W%2r5e$nHlmojMRk%IUJ&=NWJ^|e1ck_mbryMtTL{zr(qQLJ!q}f
zg;r|?0<{6D>O+4}0o?fQDtqBD{i;&I0JW>b?fbt+ekX-rN&l`=;km949Obz}89))#
zD}4SYKJ5{l2eFx{>=S?0fQX#-KkZ(fXz820^7IC{NvOuLS~Kw(vM?2;C^RTczBxtJ
z&2?fDC@+2Y2DVn~5_~g7xh$gDs!^$;Bzps0?_n6Bjy&uxxIvYE!FCYPe-&qVF;vO6
z$4d*IrPA8}sN~O(*xi-W*Rrv9l){ZsKAGP6l<ZQ;Z6N*3-@Qt_wPX56DV1-12a}bC
zB6xqqFS&6hE4N${$lqQ58?Dz-cttU_jiI#CI8z;3UZZ6r1wrYkw6k1IN=wfsrB@Z2
zAWh&wDhH(#5~*b)oK9iAbfHqkQrQf3wbkoZHAmefX?GPErArYRfbsGwPPHPR;aL@!
z+zhwt%dor!!-fJmxgOrYUYG_Oo23||o{i6kVmvP!wIj+#P4m#p;-^79uN_M7LPa%B
z^)d9&HTjfvZT@rK<-g>f)#WbV<DPJ_Zd+aZy0|V5`62(<u0~;_^1W|0co*Myc-!vs
zj_}0%+JBMJ;{()Zyks)G!GfUlf&h=6SxWPE5S-zjE_a$iAoC~)lnWs0XEYGF8Zf<R
uUdMtW?udU)vaYZlx9wf8zT#Qtz%kSFno<({F$<-8RKvJKfl<A!T<bp`%jaeQ
literal 0
HcmV?d00001
diff --git a/__pycache__/tSnS.cpython-38.pyc b/__pycache__/tSnS.cpython-38.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..45462ff99d850960856b55c316c47e6591ea6ea7
GIT binary patch
literal 5153
zcmd5AZHya7b-%p!+P=@+M=lNR=9<7ct?zPoX_7`EK3~!#C+M~0t|Bt!mW^j`ym$8Q
zW_RxsTUYW?A=C<~?GGxXC{{%TMF<I02?+^D;zOXq4@y;`Z7G2UMJ_D_1VTln;mvsM
zIKG%4=uXO5GxPS%oA<SE-p=zAEiGXN+CyKyc)HxfFz=9u`o_b<7<8u}2E-sX%cLqZ
zo8pc$h(r8&2Jz>(6pz>`CMJ}Dj+so8vA)&#>;|mHq*%nHIK;v(${{5n0R_%8sQ?P2
zAPSvlQzDXpU&IsKDPwANIxJ0Tt|Kp4vYXLlteH8@we_@Y*`}r0uC6(<IV-zMy5nj{
zhApycAfPieTXi*gJZq*i@(DGkXI<0KRbx;NOKG-FMFC(5`lp~fJ+R_3msrFS{V-o-
z^Nhk#4Cq<I@D#%%9*I7UKrne`jssr>FK~Sfz<d?PT#+pRm$|?tF;8H!KbWs@4^g}X
zm@wH-_7#B!<VFx#i*^xSk$Z?<grXu_LJ|rS*+q`5D=q73iAF%>(-Gc2Z7aYn&|7<%
z2-(|qn+9~n*PtCmV7&t|@J6DO_-pU{cBlB{(J}m?eC@v74|I&##F6;ByRMzPIMIFu
zr`Fz1fA#SrW0-x7y)qm4<XHJ8v@%DCpwo7}GTlg&iQ{9P9%ooyK(%dk#p9C0&__~S
za?}&Ex}j#V$Lq$TC%Lw2IF{*X9y{#`hH2-#K-xAPCnYYL*@c|uahC3J$6^0**nfO9
zB~|dsc*<iZV}d6*T6WeG=a*C)!>**F8_IG@Tqbk!Ta4Wc+Wh;-&=ZbkJ40$NXBxA*
zlg?=N&{55qcTH=^)p5p+pVkee8#tabk(Q17aTj+CXNbh>pta(K6yT}5iqZlsovp-P
zzLSL?-_J(bD8G$;5Pt3KHh#5*?y1tJiASm}06-A>5$Fy%Jb4lt47_X}X10szWeS4J
zU4r8WqmXCw0(refau~VQSXh!C&!~<j&#GzHwCN$yW@ojut1oKu^tep*2dWQ%HEB=L
zk-1Td1#EI0Y?2>rqJhV2%a+F@(}ksOf(Yf{6dLYA&SSG4n@dHqs_W`$O<B~^)RC2e
znC^*IMpv@Aod7%3KCz2kjZ}O^S8Re!lRQyvkdK%)=+yMRphq4u#=Pi6UFIz7vT(-W
zY{HMj0vYK^B~ayr3jV5y1uVEc3O)<&K_PMyiaeY=b`Ek)1z!+=GdIuKPooHmq88MO
z+E6>{K%Hm{>O$S92lb+@s1LUl0!BaDhVDVz(Z|re=sqN)9jG7Ok9MM6=;P=UXn@e<
zAw%>r1rcVF0Rl+sV>qS|EQA1#QA~(pq(T@lyD8=Yiis4WfO(K&;uO<TXa&q5#SBqQ
zTcI5=dnjg@Vmb<)fEl5fQHt48=*kOdFWP&7_{g7Ug#ZIIVjjHK2O+ed;h95BjgIr_
zxF&Q#^iZCG%*6snvB>w*AgpjIbl134rB%6fH<7|uq_`|(!0sYD7xL*sz6{X22q8}8
zahoVls>uL(<XAPe$k#aEXi?STW)@Z7L!*Bi`82Gr^9*Et-}~HMu4!A%r&QCr=xDGp
z-(B^0lPK|{sOpOb>Z%+OHOgstCr#xv&^3|M-qfCk{Tn1_y*(T2w+`PPiJY;!I_5iT
zF%Mb$uGSQ8Nc414pVr{b%ad?F^K&D_ehuEYX?Z#7%Mk1Gvbc$P8E&VyvRw_{H2Y8>
zzom;o`ys=?HyfY>0Ks>=i`fdjE`wE1{QvoVYc7X=PpdXV+-$JaF5jj*SXYZ4nhiG8
z)+SQxY^C4!jkGl{z&-C_sIG*(0GHr{tC!(=3&kU;9oLnCvg*PWxN?h?EmJpKd44bq
z8{*XlaMH$5S{SZ#P@b3L@+gv%q)w5SbT=cXP1A--0E$d{xxSqlY?dThmQO=5>VOIZ
zrse7;6f`m5GwvfNh9)L}FvgWh8Ab2G_6`Ek$87+uK(A{z5;nOf$5=0b+!ZK)MR=n7
zVnSs*xl7-I@b_Y`zVQ<OS_z+Q?L7IJgq6UD;*Z`R|20nF?Z0}gbNil)Kf=$P{_`_y
zFIpvhb@eaj-fiENz(;pXF8uCid<nlhb7pPnACH#sPu@EA)|(SkC7l0p$BCA;QziV1
ztJ$AF_Woo7zx@^K@)v*8Tfzg<Q?GP<|5yo+e(BUdFT4^@;OAcWLI3+)GJy-v_k8X<
zXJ0SjSD#$SmA;`SaOW>SKlAObQwjXSf381U?CVS5m%skk>wkECF@dlA@x8U5es;Ko
z%YQF?dite<34Hm@%wg?~L;~Y;4-f5<jwLYn-htm24vd!Y^pUg4zwHO$Wz{?r9{JN@
zyklnRU;91+nh$^Y^x1Dd53-*70v=IjW)fxMyE3V(%EWJFQrnb6WFYlUnbbSw2pOVe
zAjMm`6^2QhZutl*SqzWRu9m7mBF0tT0bv-Eyw;=x^&@%HkJz?p)AGqCBD`SQ%vo7&
z*$dj5lTj_r3#@8Tz}hldmC~j$w}2s*Y$7t&Q>nAqc`xWkP%2=l3t2tIXVhFuuppdL
zqN^m6${sq8*XCg#549~_Ur7aF9#?Y9K44XWTy@D>@o97ZY*EsrGEyN0fRLUfZ1z-v
zVNk`i8mqn5@e@k@pDkdFa}5SYigJT2cT;3JsED|m0^iDoIDtIfxql!={@1Dy4Uvla
z)-N=`4)FJ~GCRNxkP)hE80m8jigW^IkXtw2iB_!aT{i}0%5!kt80J>D(5m)!4@}0`
z$^i|M{{;*+rzwgTRum||7qSG7D$1D!HS4bg6$P1TMX~$9XEZkKeFP%^PSg#gvbM?J
zDkL<W5DXz+2!wh=VTy&nJ9uZs@dCy|&RVHNA)&70KpD~kJ#?s|%VCmrXo((TL;QaM
DBW$93
literal 0
HcmV?d00001
diff --git a/bands.py b/bands.py
new file mode 100644
index 0000000..da232d4
--- /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 0000000..fe6e8f4
--- /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 0000000..7bfe6be
--- /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
--
GitLab