# Data format of wannier90_hr.dat from wannier90

I did a calculation on a bulk system with DFT and fit it with the tight binding matrix with wannier90. There is question on this. I try to write my own code to construct a slab tight binding matrix from this bulk matrix. But I can't find any documentation mention how the orbital basis was arrange in this matrix. (for example, I don't know which column correspond to particular orbital and particular atom.) I can only obtain the energy spectrum from it. Is there anyone who did this before ?

The following code is a python function to read the wannier90_hr.dat, from which you should figure out its data structure.

def read_hamiltonian(path):

"""
Read hopping matrix element from the wannier90 output file:wannier90_hr.dat
wan_num: number of wannier functions
wsc_num: number of wigner-sitez cells
wsc_count: variables related to degeneracy
"""

with open(path,"r") as f:

wan_num=int(lines[1]); wsc_num=int(lines[2])
ski_row_num=int(np.ceil(wsc_num/15.0))      # skip row numbers
wsc_count=[]
for i in range(ski_row_num):
wsc_count.extend(list(map(int,lines[i+3].split())))

wsc_tot=np.zeros((wan_num**2*wsc_num,3))
tem_tot=np.zeros((wan_num**2*wsc_num,2))
for i in range(wan_num**2*wsc_num):
wsc_tot[i,:]=list(map(int,lines[3+ski_row_num+i].split()[:3]))
tem_tot[i,:]=list(map(float,lines[3+ski_row_num+i].split()[5:]))

wsc_idx=wsc_tot[0:-1:wan_num**2,:]   # the translational vector between wigner-sitez cells
hop_mat=np.reshape(tem_tot[:,0]+1j*tem_tot[:,1],[wan_num,wan_num,wsc_num],order='F')

return hop_mat,wsc_idx,wsc_count,wan_num


To run this function, just input the path to wannier90_hr.dat.

PS: this function is tested only for the Wannier1.2.

Hope it helps.

• Thank you for your reply. I did the same things before. But Im wondering how the orbital basis was arranged in the matrix. Since Im doing SOC, I don't know the matrix element correspond to which orbital, spin and atom. For example the element of the second row and third column, which hopping does it correspond to ? Apr 10, 2021 at 11:07
• You may take a look at wannier90_centres.xyz, which will tell you how to define the projection operator.
– Jack
Apr 10, 2021 at 13:14

You can read the wanniertools code. In wanniertools, to calculate surface state, they write a surfstate subroutine in surfstate.f90 file. The slab Hamiltonian is restructured from a bulk Hamiltonian in ham_qlayer2qlayer.f90 file.

enter code here
! This is a fortran code.
! H00 Hamiltonian between nearest neighbour-quintuple-layers
! the factor 2 is induced by spin

! to read the matrix from hr file
Hij=0.0d0
do iR=1,Nrpts
ia=irvec(1,iR)
ib=irvec(2,iR)
ic=irvec(3,iR)

call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic)

inew_ic= int(new_ic)
if (abs(new_ic).le.ijmax)then
kdotr=k(1)*new_ia+ k(2)*new_ib
ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr)

Hij(inew_ic, 1:Num_wann, 1:Num_wann )&
=Hij(inew_ic, 1:Num_wann, 1:Num_wann )&
+HmnR(:,:,iR)*ratio/ndegen(iR)
endif

enddo

H00new=0.0d0
H01new=0.0d0

! nslab's principle layer
! H00new
do i=1,Np
do j=1,Np
if (abs(i-j).le.(ijmax)) then
H00new(Num_wann*(i-1)+1:Num_wann*i,Num_wann*(j-1)+1:Num_wann*j)&
=Hij(j-i,:,:)
endif
enddo
enddo

! H01new
do i=1,Np
do j=Np+1,Np*2
if (j-i.le.ijmax) then
H01new(Num_wann*(i-1)+1:Num_wann*i,&
Num_wann*(j-1-Np)+1:Num_wann*(j-Np))=Hij(j-i,:,:)
endif
enddo
enddo

do i=1,Ndim
do j=1,Ndim
if(abs(H00new(i,j)-conjg(H00new(j,i))).ge.1e-4)then
!  write(stdout,*)'there are something wrong with ham_qlayer2qlayer'
!stop
endif


• This extremely short answer, which talks about what you "think", can easily fit in just a comment, and with 381 more characters left to spare. I would suggest you delete the answer and write a comment instead, as much as I'd love to have more answers here. Unless of course you can expand this answer significantly and write a couple paragraphs or give a chunk of code block like Jack did. Apr 11, 2021 at 2:04