I am trying to implement the continuum model Hamiltonian for transition metal dichalcogenide systems as described in the reference: PRL 122, 086402 (2019)
(* Define the constants*)
(* twist angle*)
\[Theta] = 1.2 Degree;
(* For MoTe2 V, psi, w*)
V = 8 (* meV*);
psi = -89.6 ;
w = -8.5 (*meV*);
Vz = -2 (* displacement field meV*);
meff = 0.62 (* valence band effective mass*);
a0 = 3.472 (*angstrom*);
NN = 6(*truncate range*);
\[Theta] =
1.2 Degree; (*\[Theta]\[Pi]/180conversion from degress to radian*)
V = V/1000;
w = w/1000;
Vz = Vz/1000;
psi = psi Degree;
hbar = 4.135667*10^-15/2/\[Pi];
me = 0.51099895*10^6;
meff = meff*me;
kin = -hbar^2/2/meff*9*10^16;
a0 = 3.472 ;
a0 = a0*10^-10;
a1 = a0 {1, 0};
a2 = a0 {Cos[\[Pi]/3], Sin[\[Pi]/3]};
(* Reciprocal lattice vectors*)
G1 = (4 \[Pi])/((Sqrt[3]) a0) {0, 1};
G2 = (4 \[Pi])/((Sqrt[3]) a0) {-Sin[\[Pi]/3], Cos[\[Pi]/3]};
(*G2=(4\[Pi])/((Sqrt[3])a0){Cos[\[Pi]/3],Sin[\[Pi]/3]};*)
g1 = \[Theta] {G1[[2]], -G1[[1]]};
g2 = \[Theta] {G2[[2]], -G2[[1]]};
kD = g1[[1]]/Sqrt[3];
Kb = {-kD*Sqrt[3]/2, -kD/2}; (* K botom layer *)
Kt = {-kD*Sqrt[3]/2,
kD/2} ; (* K top layer *)
(* Define the lattice*)
n = 6;
nlayer = (2 n + 1)^2;
L = Table[0, {i, 1, 2 (2 n + 1)^2}, {j, 1, 2}];
invL = Table[0, {i, 1, 2 n + 1}, {j, 1, 2 n + 1}];
k = 1;
count = 0;
Do[Do[
(*bottom layer*)
L[[k, 1]] = i; L[[k, 2]] = j;
invL[[i + n + 1, j + n + 1]] = count;
count = count + 1;
(* top layer*)
L[[k + nlayer, 1]] = i; L[[k + nlayer, 2]] = j;
k = k + 1;
, {j, -n, n}], {i, -n, n}]
(* check the form of inVL*)
MatrixForm[invL];
(* Hamiltonian implementation *)
Clear[kx, ky]
HHHH[kx_, ky_] :=
Module[{ h, NN},
NN = 6;
siteN = (2 NN + 1)^2;
h = ConstantArray[0, {2 siteN, 2 siteN}];
Do[
ix = L[[i, 1]];
iy = L[[i, 2]];
(* diagnal elements in the Hamiltonian*)
ax = kx - Kb[[1]] + ix g1[[1]] + iy g2[[1]];
ay = ky - Kb[[2]] + ix g1[[2]] + iy g2[[2]];
h[[i, i]] += kin (ax^2 + ay^2) + Vz/2;
ax = kx - Kt[[1]] + ix g1[[1]] + iy g2[[1]];
ay = ky - Kt[[2]] + ix g1[[2]] + iy g2[[2]];
h[[i + siteN, i + siteN]] += kin (ax^2 + ay^2) - Vz/2;
h[[i + siteN, i]] += w;
h[[i, i + siteN]] += w;
If[ix != NN,
j = invL[[ix + 1 + NN + 1, iy + NN + 1]];
h[[j + 1, i]] += V Exp[I*psi];
h[[j + siteN + 1, i + siteN]] += V Exp[-I*psi];
];
If [iy != -NN,
j = invL[[ix + NN + 1, iy - 1 + NN + 1]];
h[[j + 1, i]] += V Exp[I psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[-I psi];
h[[j + 1, i + siteN]] += w;
];
If [ix != -NN && iy != NN,
j = invL[[ix - 1 + NN + 1, iy + 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[-I*psi];
h[[j + siteN + 1, i]] += w;
];
If [ix != -NN,
j = invL[[ix - 1 + NN + 1, iy + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
];
If [iy != NN,
j = invL[[ix + NN + 1, iy + 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
h[[j + siteN + 1, i]] += w;
];
If [ix != NN && iy != -NN,
j = invL[[ix + 1 + NN + 1, iy - 1 + NN + 1]];
h[[j + 1, i]] += V*Exp[-I*psi];
h[[j + siteN + 1, i + siteN]] += V*Exp[I*psi];
h[[j + 1, i + siteN]] += w;];
, {i, 1, siteN}];
Return[h]
]
Num = 50;
KXX1 = N@Subdivide[g1[[1]]/2, 0, Num - 1];
KXX2 = N@Subdivide[0, -g1[[1]]/2, Num - 1];
KXX3 = N@Subdivide[-g1[[1]]/2, -g1[[1]]/2, Num - 1];
KXX4 = N@Subdivide[-g1[[1]]/2, g1[[1]]/2, (2 Num - 1)];
KX = Join[KXX1, KXX2, KXX3, KXX4];
KYY1 = N@Subdivide[-g1[[1]]/(Sqrt [3]/2), 0, Num - 1];
KYY2 = N@Subdivide[0, g1[[1]]/(Sqrt[3]/2), Num - 1];
KYY3 = N@
Subdivide[g1[[1]]/(Sqrt[3]/2), -g1[[1]]/(Sqrt [3]/2), Num - 1];
KYY4 = N@
Subdivide[-g1[[1]]/(Sqrt [3]/2), -g1[[1]]/(Sqrt [3]/2), (2 Num -
1)];
KY = Join[KYY1, KYY2, KYY3, KYY4];
Eigen = Table[0, {i1, Length[KX]}, {j1, 2 siteN}];
Do[Eigen[[k]] = Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1,
Length[KX]}];
out = ParallelTable[
Eigenvalues[HHHH[KX[[k]], KY[[k]]]] // Sort, {k, 1, Length[KX]}];
ListLinePlot[Transpose[out], PlotRange -> {{0, 250}, {0, 0.025}}]
The resulting band structure is look like this
I recently came across the github code which reproduce the band structure from the reference. (https://github.com/zihaophys/twisted_bilayer_graphene/blob/master/TMD/homobilayer/moireTMD.py)
I have studied and trying to reproduce in exact same way. Can someone please point out how to fix the band structure.
Thanks
