Extracting numerical results from hamiltonian using Quantica

Hi, I am using Quantica for finding the band structure of twisted bilayer graphene. Here is the part of my code that finds the hamiltonian:

b = hamiltonian(tbg, model) |>
bands(subdiv((0, 2, 3, 4), 20);
mapping=(0, 2, 3, 4) => (:Γ, :K, :M, :Γ),
solver=ES.ShiftInvert(ES.ArnoldiMethod(nev=10), 0.0))

I want to extract the numerical results of this operation, for example the energy values at the Γ point, and I can’t seem to figure out how to do that.

1 Like

If you want a set of eigenpairs at a specific point in the Brillouin zone, the easiest is to use spectrum directly, instead of trying to extract that through bands. In your case you would do

h = hamiltonian(tbg, model)
es, phis = spectrum(h, (0,0); solver=ES.ShiftInvert(ES.ArnoldiMethod(nev=10), 0.0))
1 Like

Thank you! And how would I find it for one of the other points, say K? Couldn’t figure out how to input the k-point character instead of the (0,0) so that it works…

Well, you need to replace (0,0) by (k*A1, k*A2), which for the K point is probably (2pi/3, -2pi/3) (or the negative of that for K’)

Sorry I’m still a bit confused by this. In the usual definition of the BZ, the K-point is defined as K = [2π/3, 2π/(3√3)] and neither this nor the value K = [2π/3, -2π/3] that you suggested gives the right result.
I don’t think I understand what convention the function input should follow. Also, could you please explain what you mean by the notation k, A1 and A2?

Thank you.

Your value of (ϕ1, ϕ2) at the K-point corresponds to a carbon-carbon distance a_CC = 1. My convention is for the bravais vectors of length 1, |A_1| = |A_2| = \sqrt{3}*a_CC = 1. Also, my notation corresponds to a certain orientation of the lattice, where the Carbon-Carbon bond inside the unit cell vectors is vertical (along the y-direction).

If your construction is different, your definition of K is different, and you should use the (ϕ1, ϕ2) syntax, where ϕ_i = K * A_i

If I am working with TBG would I use the superlattice unit vectors in this case for A_1?

Also, is there not a way to directly extract the position of the K-points from the TBG model? Because I noticed that for instance when constructing the hamiltonian the mapping argument takes in the K-points as mapping=(0, 2, 3, 4) => (:Γ, :K, :M, :Γ).

There is a relation between the two Bravais vectors A1, A2 in TBG and the monolayer Bravais vectors that depends on twist angle (or more precisely, on twist indices). The Bloch phases (ϕ1, ϕ2) of the K-point is orthogonal (relative to the Gamma point) to the Bloch phases of the untwisted monolayer. Hence, while it is common to use (2pi/3, -2pi/3) for the latter, the corresponding one is (2pi/3, 2pi/3) for the former, I believe. But do check it out.

A quick way to confirm is to build your system and actually plot the banstructure for a largish twist angle. You will quickly be able to identify the K/K´ values by the position of the Dirac points

1 Like