Mathematica Asked on August 27, 2021
I am trying to plot band structure of twisted bilayer graphene.
This is a ten bands model studied in this work reference.
Everything is working fine but the bands are taking jump among each other for no reason.
I have tried some solutions but none of them have worked. If anyone have a nicer approach that also takes optimum time to plot these bands.
Because with these bad visualization it takes almost one minute to plot it.
This is the code.
The other approaches could be calling, Method -> {"Arnoldi", "Criteria" -> "RealPart"}
. But this only works for 8 bands not for 10 bands.
taz = 0;
tp = 0.003;
tpp = 0.004;
tmp = 0.004;
tk = 0;
tkb = 0;
tpppz = 0.016;
tmppz = 0.016;
tpkp = 0.016;
tmkp = -0.016;
[Mu]az = -0.1 - 6*taz;
[Mu]aa = 0.11 - 4 (tk + tkb);
[Mu]dw = 0.11 - 4 (tk + tkb);
A = 0.110;
b = 0.033;
c = 0.033;
d = 0.573;
g = [Pi]/8;
te = I;
AA[px_, py_, a_] :=
te*Exp[-I*g] (1 + [Phi]1[px, py, a, 1] + [Phi]2[px, py, a, -1]);
AA1[px_, py_, a_] := Conjugate[AA[px, py, a]]
[Phi]1[px_, py_, a_, l_] := Exp[-I*l*Sqrt[3]/2 a (-px + Sqrt[3] py)]
[Phi]2[px_, py_, a_, l_] := Exp[-I*l*Sqrt[3]/2 a (px + Sqrt[3] py)]
[Phi]3[px_, py_, a_, l_] := Exp[-I*l*3 a*py]
W[l_] := Exp[I*l*(2*[Pi])/3]
[Xi][l_] := Exp[I*l*(2*[Pi])/6]
H1[px_, py_, a_] := ([Phi]2[px, py, a, 1] + [Phi]3[px, py, a, 1] + [Phi]1[px,
py, a, 1] + [Phi]2[px, py, a, -1] + [Phi]3[px, py,
a, -1] + [Phi]1[px, py, a, -1])
R1[px_, py_, a_] := {0, 0, 0, 0, 0,
0, -(W[1] + [Phi]3[px, py, a, 1]*W[-1] + [Phi]1[px, py, a,
1]) [Xi][-1]*
A, (1 + [Phi]3[px, py, a, 1]*W[-1] + [Phi]1[px, py, a, 1]*
W[1]) [Xi][-1]*
A, (W[1] + [Phi]2[px, py, a, 1]*W[-1] + [Phi]3[px, py, a,
1]) [Xi][1]*
A, -(W[1] + [Phi]2[px, py, a, 1] + [Phi]3[px, py, a, 1]*
W[-1]) [Xi][1]*A}
R2[px_, py_, a_] := {0, 0, 0, 0, 0,
0, (W[-1] + [Phi]3[px, py, a, 1]*W[1] + [Phi]1[px, py, a,
1]) [Xi][1]*
b, (1 + [Phi]3[px, py, a, 1] + [Phi]1[px, py, a,
1]) c, (W[-1] + [Phi]2[px, py, a, 1]*W[1] + [Phi]3[px, py, a,
1]) [Xi][-1]*
b, (1 + [Phi]2[px, py, a, 1] + [Phi]3[px, py, a, 1]) c}
R3[px_, py_, a_] := {0, 0, 0, 0, 0,
0, (1 + [Phi]3[px, py, a, 1] + [Phi]1[px, py, a,
1]) c, (1 + [Phi]3[px, py, a, 1]*W[1] + [Phi]1[px, py, a, 1]*
W[-1]) [Xi][1]*
b, (1 + [Phi]2[px, py, a, 1] + [Phi]3[px, py, a,
1]) c, (W[-1] + [Phi]2[px, py, a, 1] + [Phi]3[px, py, a, 1]*
W[1]) [Xi][-1]*b}
R4[px_, py_, a_] := {0, 0, 0, 0, 0,
0, -I*d*[Phi]2[px, py, a, -1], -I*d*[Phi]2[px, py, a, -1], I*d,
I*d}
R5[px_, py_, a_] := {0, 0, 0, 0, 0, 0, -I*d*W[1], -I*d*W[-1],
I*d*W[1], I*d*W[-1]}
R6[px_, py_, a_] := {0, 0, 0, 0, 0,
0, -I*W[-1]*d*[Phi]1[px, py, a, 1], -I*W[1]*
d*[Phi]1[px, py, a, 1], I*W[-1]*d, I*W[1]*d}
R7[px_, py_,
a_] := {-(W[-1] + [Phi]3[px, py, a, -1]*W[1] + [Phi]1[px, py,
a, -1]) [Xi][
1] A, (W[1] + [Phi]3[px, py, a, -1]*W[-1] + [Phi]1[px, py,
a, -1]) b*[Xi][-1], (1 + [Phi]3[px, py, a, -1] + [Phi]1[px,
py, a, -1]) c, I*d*[Phi]2[px, py, a, 1], I*W[-1]*d,
I*W[1]*d*[Phi]1[px, py, a, -1], 0, 0, 0, 0}
R8[px_, py_,
a_] := {(1 + [Phi]3[px, py, a, -1]*W[1] + [Phi]1[px, py, a, -1]*
W[-1]) [Xi][
1] A, (1 + [Phi]3[px, py, a, -1] + [Phi]1[px, py,
a, -1]) c, (1 + W[-1]*[Phi]3[px, py, a, -1] +
W[1]*[Phi]1[px, py, a, -1]) [Xi][-1]*b,
I*d*[Phi]2[px, py, a, 1], I*W[1]*d,
I*W[-1]*d*[Phi]1[px, py, a, -1], 0, 0, 0, 0}
R9[px_, py_,
a_] := {(W[-1] + [Phi]2[px, py, a, -1]*W[1] + [Phi]3[px, py,
a, -1]) [Xi][-1] A, (W[1] + [Phi]2[px, py, a, -1]*
W[-1] + [Phi]3[px, py, a, -1]) b*[Xi][
1], (1 + [Phi]2[px, py, a, -1] + [Phi]3[px, py, a, -1]) c, -I*
d, -I*W[-1]*d, -I*W[1]*d, 0, 0, 0, 0}
R10[px_, py_,
a_] := {-(W[-1] + [Phi]2[px, py, a, -1] + [Phi]3[px, py, a, -1]*
W[1]) [Xi][-1] A, (1 + [Phi]2[px, py, a, -1] + [Phi]3[px,
py, a, -1]) c, (W[1] + [Phi]2[px, py,
a, -1] + [Phi]3[px, py, a, -1]*W[-1]) b*[Xi][1], -I*d, -I*
W[1]*d, -I*W[-1]*d, 0, 0, 0, 0}
r1[px_, py_,
a_] := {[Mu]az, (-I*
tpppz ([Phi]1[px, py, a, -1] + [Phi]3[px, py, a, 1]*
W[-1] + [Phi]2[px, py, a, -1] W[1]) - (-1) I*
tmppz ([Phi]1[px, py, a, 1] + [Phi]3[px, py, a, -1]*
W[-1] + [Phi]2[px, py, a, 1] W[1])), (-(-1) I*
tpppz ([Phi]1[px, py, a, 1] + [Phi]3[px, py, a, -1]*
W[1] + [Phi]2[px, py, a, 1] W[-1]) - (-1) (-1) I*
tmppz ([Phi]1[px, py, a, -1] + [Phi]3[px, py, a, 1]*
W[1] + [Phi]2[px, py, a, -1] W[-1])), 0, 0, 0, 0, 0, 0, 0}
r2[px_, py_,
a_] := {(I*
tpppz ([Phi]1[px, py, a, 1] + [Phi]3[px, py, a, -1]*
W[1] + [Phi]2[px, py, a, 1] W[-1]) -
I*tmppz ([Phi]1[px, py, a, -1] + [Phi]3[px, py, a, 1]*
W[1] + [Phi]2[px, py, a, -1] W[-1])),
tp*H1[px, py, a] + [Mu]aa,
tpp ([Phi]1[px, py,
a, -1] + [Phi]3[px, py, a, 1] W[-1] + [Phi]2[px, py,
a, -1] W[1]) +
tmp ([Phi]1[px, py, a,
1] + [Phi]3[px, py, a, -1] W[-1] + [Phi]2[px, py, a, 1] W[
1]), tpkp*[Phi]2[px, py, a, 1] - tmkp*[Phi]3[px, py, a, 1],
tpkp*[Phi]3[px, py, a, 1]*W[1] - tmkp*W[1],
tpkp*W[-1] - tmkp*[Phi]2[px, py, a, 1]*W[-1], 0, 0, 0, 0}
r3[px_, py_,
a_] := {(-I*
tpppz ([Phi]1[px, py, a, -1] + [Phi]3[px, py, a, 1]*
W[-1] + [Phi]2[px, py, a, -1] W[1]) - (-1) I*
tmppz ([Phi]1[px, py, a, 1] + [Phi]3[px, py, a, -1]*
W[-1] + [Phi]2[px, py, a, 1] W[1])),
tpp ([Phi]1[px, py, a,
1] + [Phi]3[px, py, a, -1] W[1] + [Phi]2[px, py, a,
1] W[-1]) +
tmp ([Phi]1[px, py,
a, -1] + [Phi]3[px, py, a, 1] W[1] + [Phi]2[px, py,
a, -1] W[-1]), tp*H1[px, py, a] + [Mu]aa,
tpkp*[Phi]3[px, py, a, 1] - tmkp*[Phi]2[px, py, a, 1],
tpkp*W[-1] - tmkp*[Phi]3[px, py, a, 1]*W[-1],
tpkp*[Phi]2[px, py, a, 1]*W[1] - tmkp*W[1], 0, 0, 0, 0}
r4[px_, py_, a_] := {0,
tpkp*[Phi]2[px, py, a, -1] - tmkp*[Phi]3[px, py, a, -1],
tpkp*[Phi]3[px, py, a, -1] - tmkp*[Phi]2[px, py, a, -1], [Mu]dw,
0, 0, 0, 0, 0, 0}
r5[px_, py_, a_] := {0,
tpkp*[Phi]3[px, py, a, -1]*W[-1] - tmkp*W[-1],
tpkp*W[1] - tmkp*[Phi]3[px, py, a, -1]*W[1], 0, [Mu]dw, 0, 0, 0,
0, 0}
r6[px_, py_, a_] := {0, tpkp*W[1] - tmkp*[Phi]2[px, py, a, -1]*W[1],
tpkp*[Phi]2[px, py, a, -1]*W[-1] - tmkp*W[-1], 0, 0, [Mu]dw, 0, 0,
0, 0}
r7[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, 0, AA[px, py, a], 0}
r8[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, 0, 0, AA[px, py, a]}
r9[px_, py_, a_] := {0, 0, 0, 0, 0, 0, AA1[px, py, a], 0, 0, 0}
r10[px_, py_, a_] := {0, 0, 0, 0, 0, 0, 0, AA1[px, py, a], 0, 0}
hh[px_, py_, a_] := {R1[px, py, a] + r1[px, py, a],
R2[px, py, a] + r2[px, py, a], R3[px, py, a] + r3[px, py, a],
R4[px, py, a] + r4[px, py, a], R5[px, py, a] + r5[px, py, a],
R6[px, py, a] + r6[px, py, a], R7[px, py, a] + r7[px, py, a],
R8[px, py, a] + r8[px, py, a], R9[px, py, a] + r9[px, py, a],
R10[px, py, a] + r10[px, py, a]}
ev4 = Plot[
Eigenvalues[
hh[Cos[[Theta]], Sin[[Theta]], 1.4]], {[Theta], -0.5 [Pi],
1.5 [Pi]}]
I can provide an answer for the fixing the visual. We'll use some dynamically scoped caching + add a function that only evaluates once Plot
is ready to go:
Block[{getEv},
getEv[θ_?NumericQ] :=
getEv[θ] = Sort@Eigenvalues[hh[Cos[θ], Sin[θ], 1.4]];
getEv[θ_?NumericQ, n_] :=
getEv[θ][[n]];
Plot[Evaluate@Table[getEv[θ, n], {n, 10}], {θ, -0.5 π, 1.5 π}]
]
The jumps happen because Mathematica returns the eigenvalues ordered by magnitude, i.e. it ignores the sign. You can plot the sorted eigenvalues to smooth that out
Block[{getEv},
getEv[θ_?NumericQ] :=
getEv[θ] = Sort@Eigenvalues[hh[Cos[θ], Sin[θ], 1.4]];
getEv[θ_?NumericQ, n_] :=
getEv[θ][[n]];
Plot[Evaluate@Table[getEv[θ, n], {n, 10}], {θ, -0.5 π, 1.5 π}]
]
Note that this leads to avoided crossings, where the bands could actually switch their energy ordering. I'll leave it up to you how you want to resolve this.
Correct answer by b3m2a1 on August 27, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP