# How to interpret WannierTools Python code that generates tight-binding model in the wannier90_hr.dat format?

I am getting myself acquainted with Wannier Tools. Wannier Tools requires two inputs, a wt.in file, and a .dat file. This .dat file should have the structure explained in the manual here. This file may be parsed from other DFT software (such as Quantum ESPRESSO via Wannier90). However, it may also be generated manually, as the developers do in an example using a Python file. I am going to focus on their implementation of the 2D Haldane model for the QAHE (link to implementation; the original paper). My issue is that I do not understand how they went from the tight-binding model to getting the matrix elements of the Hamiltonian. In particular, I do not understand why they use 7 points, and why each point gives each matrix element.

Here is the tight-binding model in first quantization, from this source (I used another source because the original work only gives a Hamiltonian in k-space): Here is part of an example Python file that they use to generate the .dat file. They use $$m$$ instead of $$M$$, $$phi$$ instead of $$\phi$$, and $$t1$$ instead of $$t$$:

# Chern = 1
m=0.2; phi= np.pi/2.0; t1=1.0; t2=m/3.0/np.sqrt(3)*2.0;

# maximum dimension for hr matrix
ndim = 2
nrpts = 7
num_patom=2

# hr matrix
norbs = num_patom*1
hmnr= np.zeros((norbs,norbs,nrpts),dtype = np.complex128)

# WS points
irvec   = np.zeros((3,nrpts),dtype = np.int32)

# degeneracy
dege   = np.zeros((nrpts),dtype = np.int32)+1

# complex unit
zi=1j

ir= 0
irvec[0, ir]= 0
irvec[1, ir]= 0
hmnr[0, 0, ir]=  m
hmnr[1, 1, ir]= -m
hmnr[0, 1, ir]=  t1
hmnr[1, 0, ir]=  t1

# 1 0
ir= ir+1
irvec[0, ir]= 1
irvec[1, ir]= 0
hmnr[0, 0, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2
hmnr[0, 1, ir]= t1

# 0 1
ir= ir+1
irvec[0, ir]= 0
irvec[1, ir]= 1
hmnr[0, 0, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2
hmnr[1, 0, ir]= t1

# 1 1
ir= ir+1
irvec[0, ir]= 1
irvec[1, ir]= 1
hmnr[0, 0, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2

#-1 0
ir= ir+1
irvec[0, ir]=-1
irvec[1, ir]= 0
hmnr[0, 0, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2
hmnr[1, 0, ir]= t1

# 0-1
ir= ir+1
irvec[0, ir]= 0
irvec[1, ir]=-1
hmnr[0, 0, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2
hmnr[0, 1, ir]= t1

#-1-1
ir= ir+1
irvec[0, ir]=-1
irvec[1, ir]=-1
hmnr[0, 0, ir]= (np.cos(phi)-zi*np.sin(phi)) *t2
hmnr[1, 1, ir]= (np.cos(phi)+zi*np.sin(phi)) *t2

#print "dump hr.dat..."
with open('Haldane_hr.dat','w') as f:
line="Haldane model with m="+str(m)+", phi="+str(phi/np.pi)+"pi, t1="+str(t1)+", t2="+str(t2)+"Ref:Physical Review Letters 61, 18(1988)"+'\n'
f.write(line)
nl = np.int32(np.ceil(nrpts/15.0))
f.write(str(norbs)+'\n')
f.write(str(nrpts)+'\n')
for l in range(nl):
line="    "+'    '.join([str(np.int32(i)) for i in dege[l*15:(l+1)*15]])
f.write(line)
f.write('\n')
for irpt in range(nrpts):
rx = irvec[0,irpt];ry = irvec[1,irpt];rz = irvec[2,irpt]
for jatomorb in range(norbs):
for iatomorb in range(norbs):
rp =hmnr[iatomorb,jatomorb,irpt].real
ip =hmnr[iatomorb,jatomorb,irpt].imag
line="{:8d}{:8d}{:8d}{:8d}{:8d}{:20.10f}{:20.10f}\n".format(rx,ry,rz,jatomorb+1,iatomorb+1,rp,ip)
f.write(line)


Here is my attempt at understanding the above implementation:

I am not sure what the seven cases indicated by comments "# 1 0", "#-1 0", ... are. However, I figured they corresponded to real space (x,y) points because there are 7 such code blocks, and there is "nrpts = 7". Assuming that each hmnr value corresponds to a matrix element of the Hamiltonian (i.e. it can be a hopping, on-site energy, etc), I do not see why we need 7 points when we could have used 4 to capture all on-site, NN and NNN potentials/interactions. Indeed, it appears as if 3 of these code blocks have values identical to another block. So, why 7? It makes the resulting matrix larger...

Next, I figured that the $$\cos\phi \pm i \sin\phi$$ came from the Peierls substitution for the NNN hopping parameter $$t_2 \rightarrow t_2 e^{\pm i \phi} = t_2 (\cos\phi \pm i \sin\phi)$$. This checks out.

But, what still doesn't make sense is how each of the 7 points makes them assign only some matrix elements. To exemplify this, consider the first case above

ir= 0
irvec[0, ir]= 0
irvec[1, ir]= 0
hmnr[0, 0, ir]=  m
hmnr[1, 1, ir]= -m
hmnr[0, 1, ir]=  t1
hmnr[1, 0, ir]=  t1


First, I assumed the above to have the following meanings:

ir= a number that iterates through all 7 points
irvec[0, ir]= set x coordinate of point (0,0) to 0
irvec[1, ir]= set y coordinate of point (0,0) to 0
hmnr[0, 0, ir]= set (0,0) Hamiltonian matrix element value to +m
hmnr[1, 1, ir]= set (1,1) Hamiltonian matrix element value to -m
hmnr[0, 1, ir]= set (0,1) Hamiltonian matrix element value to t1
hmnr[1, 0, ir]= set (1,0) Hamiltonian matrix element value to t1


I took this to mean that the Hamiltonian at this point is $$\begin{bmatrix} \langle m | H_{00} | m \rangle & \langle m | H_{01} | n \rangle \\ \langle n | H_{10} | m \rangle & \langle n | H_{11} | n \rangle \end{bmatrix} = \begin{bmatrix} m & t_1 \\ t_1 & -m \end{bmatrix}.$$

But, how exactly do they determine the values assigned to each matrix element? In addition to the above, I don't understand why they point $$(0,1)$$ gives $$\begin{bmatrix} t_2 (\cos\phi + i \sin\phi) & 0 \\ t_1 & t_2 (\cos\phi - i \sin\phi) \end{bmatrix}.$$ What justifies setting the $$\langle m | H_{01} | n \rangle$$ element to 0? Clearly, I have either misinterpreted what the code stands for, or just don't understand what's going on. I also considered interpreting aforementioned seven points as pairs of values of $$(i,j)$$ or some other vector. But it's very hard to decrypt this without knowing too well the context the authors wrote this in. So, any help would be appreciated! I tried comparing this syntax with other tools such as PythTB and TBmodels, but they don't match too well with this syntax.