Plotting a bandstructure along high-symmetry paths in the Brillouin Zone I am trying to finalize my band structure plot for twisted bilayer graphene. I have been having some problems with the plot itself.
To troubleshoot my code, I first just looked at the diagonal part of my Hamiltonian matrix which consists of block diagonal single-layer graphene Dirac Hamiltonians.
The symmetry path I chose for the Moiré BZ in its zigzag direction was $\Gamma \rightarrow K' \rightarrow K \rightarrow \Gamma \rightarrow M $. I choose around 100 wave-vector values between each symmetry point, however when I plot the bands, I get the correct dispersion along the $K'$ to $K$ path, while the others just look weird. See below...

I think this has to do with the actual scaling of the symmetry path. I literally took the values of the corners in momentum space, and that is how I constructed the path, but I am not sure if there's another convention on plotting it. Have you had this same issue, or know of a specific way to plot the bands against the symmetry paths?
Sorry if this question is too simple...
Edit: The Hamiltonian I am using for this system is one with intralayer Dirac terms and the interlayer terms "turned off". I am following the Bistritzer and Macdonald model.
$H = \hbar v_F \sigma^{\theta, \mu}_\zeta \cdot (q - \zeta K^\mu_\zeta + G_i)$
Here the $\sigma$ is the rotated Pauli operator: $e^{-i\theta/2}\sigma e^{i\theta/2}$ and the specific rotated Dirac points are the following. $K$ are the rotated Dirac points and hence the corners of the Moire BZ and the $G_i$ can be a sum of the Moire reciprocal lattice vectors.
For example:
The new Moire BZ corners which I will call $K^1$ and $K^2$ are
$K^1 = \frac{4 \pi }{3 \sqrt{3} a}(\cos(\theta/2), \sin(\theta/2))$
$K^2 = \frac{4 \pi }{3 \sqrt{3} a}(\cos(\theta/2), -\sin(\theta/2))$
When I was choosing my path I used the $K^2$ vector as reference. Using what @Franz had suggested I went from $\Gamma-K^2-M-\Gamma$ and I coded the path in the following way:
KX = list(np.linspace(0, K_2[0]*np.sqrt(3)/2, num))+ list(np.linspace(K_2[0]*np.sqrt(3)/2,K_2[0]*np.sqrt(3)/2, 50)) + list(np.linspace(K_2[0]*np.sqrt(3)/2,0,86))
KY = list(np.linspace(0, -K_2[1]/2, num)) + list(np.linspace(-K_2[1]/2,K_1[1]/2, 50)) + list(np.linspace(K_1[1]/2,0, 86))


From here I feed these KX and KY to the Hamiltonian function I wrote, and then will obtain a (len(KX), 4*n) array for the respective eigenvalues.
 A: I can give you some workarounds in order to plot your electron dispersion along high symmetry paths.
The most precise thing to do would be to calculate the absolute value of your k-vector and plot the eigenenergies against it. I had some issues with that, which is why I can give you an even simpler solution to your problem.
Let us consider your hexagonal brillouin zone in 2D. Let us furthermore think of the high symmetry path $\Gamma$-$\mathrm{K}$-$\mathrm{M}$-$\Gamma$ which I drew into the hexagon below for better imagination.

From basic geometry we know that the distance between the center of the hexagon ($\Gamma$) and the corner ($K$) is as long as the sides of the hexagon. That basically gave us all information about the relative lengths between our high symmetry points.
The only thing that is left for you to do is to choose the number of points between two high symmetry points in such a way that they fit the corresponding relative lengths of the high symmetry lines. With that you can plot your computed eigenvalues against the k-index, so basically the natural number of the point on the high symmetry line.
Say between $\Gamma$-$\mathrm{K}$ you choose $100$ points, than $\mathrm{K}$-$\mathrm{M}$ would be half of it, so $50$ points. $\mathrm{M}$-$\Gamma$ is $\frac{\sqrt{3}}{2}$ of the length of $\Gamma$-$\mathrm{K}$, so you use $\frac{\sqrt{3}}{2} \cdot 100$ points for that symmetry line.
I would also recommend you to use matplotlib.pyplot.scatter for plotting instead of the ordenary plt.plot(). Sometimes plt.plot() does not work properly for bandstructures, as you find many energies at one k-point. Another advantage of matplotlib.pyplot.scatter is that you can plot other values, such as the spin's z-expectation value, as a color on top of the energies. I therefore attached the documentation of this matplotlib function for you.
https://matplotlib.org/3.3.4/api/_as_gen/matplotlib.pyplot.scatter.html
A: This doesn't look like a problem with plotting. Look closer (ideally, zoom it in) at this part of your plot:

Notice how in the direction $K-\Gamma$ you have two periods of some oscillatory behavior of the yellow and gray curves, and then they simply diverge away from the center. This looks more like you have set up your Hamiltonian, passing $\vec k$ outside of the first Brillouin zone, and, since the Hamiltonian is a finite-dimensional matrix, rather than the full infinite-dimensional operator, you've pushed it beyond its limits of applicability.
The two oscillation periods also suggest that you may be going too fast across the Brillouin zone. In fact, this might be the sole problem with your code: if you stop at the half of the first period, it might appear to already be the $\Gamma$ point:

