fix bugs in the surface state calculation

src/Makefile

@@ -0,0 +1,47 @@
+OBJ =  module.o sparse.o wt_aux.o math_lib.o symmetry.o readHmnR.o inverse.o proteus.o \
+       eigen.o ham_qlayer2qlayer.o psi.o unfolding.o rand.o \
+  		 ham_slab.o ham_bulk.o ek_slab.o ek_bulk_polar.o ek_bulk.o \
+       readinput.o fermisurface.o surfgreen.o surfstat.o \
+  		 mat_mul.o ham_ribbon.o ek_ribbon.o \
+       fermiarc.o berrycurvature.o \
+  		 wanniercenter.o dos.o  orbital_momenta.o \
+  		 landau_level_sparse.o landau_level.o lanczos_sparse.o \
+		 berry.o wanniercenter_adaptive.o \
+  		 effective_mass.o findnodes.o \
+		 sigma_OHE.o sigma.o Boltz_transport_anomalous.o \
+	main.o
+
+# compiler
+F90  = mpiifort -fpp -DMPI -DINTELMKL -fpe3
+
+INCLUDE = -I${MKLROOT}/include
+WFLAG = -nogen-interface
+OFLAG = -O3 -static-intel
+FFLAG = $(OFLAG) $(WFLAG)
+LFLAG = $(OFLAG)
+
+# ARPACK LIBRARY
+ARPACK=/home/yanyi2024/0-tools/0-Basic/ARPACK/libarpack_MAC.a
+
+# blas and lapack libraries
+# static linking
+#LIBS = -Wl,--start-group ${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \
+        ${MKLROOT}/lib/intel64/libmkl_sequential.a \
+        ${MKLROOT}/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread -lm -ldl \
+		  ${ARPACK}
+
+# dynamic linking
+LIBS = ${ARPACK} -L/${MKLROOT}/lib/intel64 -lmkl_core -lmkl_sequential -lmkl_intel_lp64 -lpthread
+
+main : $(OBJ)
+	$(F90) $(LFLAG) $(OBJ) -o wt2.x $(LIBS)
+	cp -f wt2.x ../bin
+
+.SUFFIXES: .o .f90
+
+.f90.o :
+	$(F90) $(FFLAG) $(INCLUDE) -c $*.f90
+
+clean :
+	rm -f *.o *.mod *~ wt2.x
+

src/ham_qlayer2qlayer.f90

@@ -17,6 +17,7 @@ subroutine ham_qlayer2qlayer(k,H00new,H01new)
      ! new index used to sign irvec     
      real(dp) :: new_ia,new_ib,new_ic
      integer :: inew_ic
+     integer :: int_ic
 
      ! wave vector k times lattice vector R  
      real(Dp) :: kdotr  
@@ -44,21 +45,19 @@ subroutine ham_qlayer2qlayer(k,H00new,H01new)
      allocate(Hij(Num_wann,Num_wann,-ijmax:ijmax))
 
      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
+     do iR=1,nrpts_surfacecell
+        ia=irvec_surfacecell(1,iR)
+        ib=irvec_surfacecell(2,iR)
+        ic=irvec_surfacecell(3,iR)
+
+        int_ic= int(ic)
+        if (abs(ic).le.ijmax)then
+           kdotr=k(1)*ia+ k(2)*ib
            ratio=cos(2d0*pi*kdotr)+zi*sin(2d0*pi*kdotr)
 
-           Hij(1:Num_wann, 1:Num_wann, inew_ic )&
-           =Hij(1:Num_wann, 1:Num_wann, inew_ic)&
-           +HmnR(:,:,iR)*ratio/ndegen(iR)
+           Hij(1:Num_wann, 1:Num_wann, int_ic )&
+           =Hij(1:Num_wann, 1:Num_wann, int_ic)&
+           +HmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR)
         endif
 
      enddo
@@ -443,7 +442,7 @@ subroutine ham_qlayer2qlayer2(k,Hij)
      implicit none
 
      ! loop index
-     integer :: iR, inew_ic
+     integer :: iR, inew_ic, int_ic
      real(dp) :: ia, ib, ic
 
      ! new index used to sign irvec     
@@ -460,22 +459,21 @@ subroutine ham_qlayer2qlayer2(k,Hij)
      complex(Dp), intent(out) :: Hij(-ijmax:ijmax,Num_wann,Num_wann)
 
      Hij=0.0d0
-     do iR=1,Nrpts
-        ia=irvec(1,iR)
-        ib=irvec(2,iR)
-        ic=irvec(3,iR)
 
-        !> new lattice
-        call latticetransform(ia, ib, ic, new_ia, new_ib, new_ic)
+     do iR=1,nrpts_surfacecell
+        ia=irvec_surfacecell(1,iR)
+        ib=irvec_surfacecell(2,iR)
+        ic=irvec_surfacecell(3,iR)
 
-        inew_ic= int(new_ic)
-        if (abs(new_ic).le.ijmax)then
-           kdotr=k(1)*new_ia+k(2)*new_ib
+        int_ic= int(ic)
+        if (abs(ic).le.ijmax)then
+           kdotr=k(1)*ia+k(2)*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)
+           Hij(int_ic, 1:Num_wann, 1:Num_wann )&
+           =Hij(int_ic, 1:Num_wann, 1:Num_wann )&
+           +HmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR)
+
         endif
 
      enddo

src/module.f90

@@ -897,11 +897,18 @@ module para
 
      integer, allocatable     :: ndegen(:)  ! degree of degeneracy of R point
 
-     complex(dp), allocatable :: HmnR_newcell(:,:,:)   ! Hamiltonian m,n are band indexes
+     complex(dp), allocatable :: HmnR_surfacecell(:,:,:)   ! Hamiltonian m,n are band indexes
      real(dp), allocatable :: Atom_position_cart_newcell(:,:)   ! Hamiltonian m,n are band indexes
      real(dp), allocatable :: Atom_position_direct_newcell(:,:)   ! Hamiltonian m,n are band indexes
-     integer, allocatable     :: irvec_newcell(:,:)   ! R coordinates
-     integer, allocatable     :: ndegen_newcell(:)  ! degree of degeneracy of R point
+     integer, allocatable     :: irvec_surfacecell(:,:)   ! R coordinates
+     integer, allocatable     :: ndegen_surfacecell(:)  ! degree of degeneracy of R point
+
+     real(dp), allocatable     :: Rmn_old(:)
+     real(dp), allocatable     :: Rmn_new(:)
+     real(dp), allocatable     :: irvec_new(:)
+     integer, allocatable     :: irvec_new_int(:)
+     integer, allocatable     :: nrpts_surfacecell
+
      real(dp),public, save :: Rua_newcell(3) !> three rotated primitive vectors in old coordinate system
      real(dp),public, save :: Rub_newcell(3) !> three rotated primitive vectors in old coordinate system
      real(dp),public, save :: Ruc_newcell(3) !> three rotated primitive vectors in old coordinate system

src/readHmnR.f90

@@ -337,7 +337,7 @@ subroutine get_hmnr_cell(cell)
    real(dp) :: shift_vec_direct(3)
 
    !> for newcell
-   real(dp) :: apos1d(3),apos2d(3)
+   real(dp) :: apos1d(3),apos2d(3),apos1dprime(3),apos2dprime(3)
    !>count newcell nrpts
    integer :: max_ir
    integer :: nir1,nir2,nir3, ir_cell
@@ -349,116 +349,152 @@ subroutine get_hmnr_cell(cell)
    integer, allocatable :: allirs(:, :, :, :)
    integer :: irn1(3),irn2(3)
 
+   !> The Allocate Process
+   integer :: ia1,ia2,ia1prime,ia2prime
+
    max_ir=8
    nrpts_max=(2*max_ir+1)**3
    allocate( rpts_array(-max_ir:max_ir,-max_ir:max_ir,-max_ir:max_ir))
    allocate( rpts_map(-max_ir:max_ir,-max_ir:max_ir,-max_ir:max_ir))
 
    allocate(allirs(nrpts, num_wann, num_wann, 3))
 
-   call cart_direct_real(shift_to_topsurface_cart, shift_vec_direct, cell%lattice)
+   allocate(Rmn_old(3))
+   allocate(Rmn_new(3))
+   allocate(irvec_new_int(3))
+
+   ! call cart_direct_real(shift_to_topsurface_cart, shift_vec_direct, cell%lattice)
 
    call date_and_time(DATE=date_now,ZONE=zone_now, TIME=time_now)
    !> Get new Hrs
    rpts_array=0
-   nrpts_new=0
+   nrpts_surfacecell=0
    rpts_map= 0
 
-   !> get number of R points for the new cell first 
+   !> 1. Get number of R points for the new cell
    do ir=1, nrpts
-      do i=1, Num_wann
-         do j=1, Num_wann
-            apos1d= Origin_cell%wannier_centers_direct(:, i)- shift_vec_direct
-            apos2d= Origin_cell%wannier_centers_direct(:, j)- shift_vec_direct+irvec(:, ir)
-            call latticetransform(apos1d(1),apos1d(2),apos1d(3),new_ia,new_ib,new_ic)
-            irn1=floor([new_ia,new_ib,new_ic])
-
-            call latticetransform(apos2d(1),apos2d(2),apos2d(3),new_ia,new_ib,new_ic)
-            irn2=floor([new_ia,new_ib,new_ic])
-
-            nir1=irn2(1)-irn1(1)
-            nir2=irn2(2)-irn1(2)
-            nir3=irn2(3)-irn1(3)
-            if (abs(nir1)>max_ir .or. abs(nir2)>max_ir .or. abs(nir3)>max_ir) cycle
-            rpts_array(nir1,nir2,nir3)=1
+      do j=1, Num_wann
+         do i=1, Num_wann
+            ia1 = Origin_cell%spinorbital_to_atom_index(i)
+            ia2 = Origin_cell%spinorbital_to_atom_index(j)
+            apos1d = Origin_cell%Atom_position_direct(:, ia1)
+            apos2d = Origin_cell%Atom_position_direct(:, ia2)
+
+            ia1prime = Cell_defined_by_surface%spinorbital_to_atom_index(i)
+            ia2prime = Cell_defined_by_surface%spinorbital_to_atom_index(j)
+            apos1dprime = Cell_defined_by_surface%Atom_position_direct(:,ia1prime)
+            apos2dprime = Cell_defined_by_surface%Atom_position_direct(:,ia2prime)
+
+            Rmn_old = irvec(:, ir) + apos2d - apos1d
+            call latticetransform(Rmn_old(1),Rmn_old(2),Rmn_old(3),Rmn_new(1),Rmn_new(2),Rmn_new(3))
+            irvec_new = Rmn_new - (apos2dprime - apos1dprime)
+
+            !> Due to the accuracy of computing, need to rounding  (but always tiny)
+            irvec_new_int(1) = ANINT(irvec_new(1))
+            irvec_new_int(2) = ANINT(irvec_new(2))
+            irvec_new_int(3) = ANINT(irvec_new(3))
+
+            if (abs(irvec_new_int(1))>max_ir .or. abs(irvec_new_int(2))>max_ir .or. abs(irvec_new_int(3))>max_ir) cycle
+            rpts_array(irvec_new_int(1),irvec_new_int(2),irvec_new_int(3))=1
+
          enddo
       enddo
    enddo
 
-   !> find all irvec
-   nrpts_new= sum(rpts_array)
-   allocate(irvec_newcell(3, Nrpts_new))
-   iter= 0
+   !> The total number of lattice points searched above
+   nrpts_surfacecell= sum(rpts_array)
+
+   allocate(irvec_surfacecell(3, nrpts_surfacecell))
+   irvec_surfacecell = 0
+
+   !> 2. Create an order map to sign the R points generated in step1
+
+   iter = 0
    do nir3=-max_ir,max_ir
       do nir2=-max_ir,max_ir
          do nir1=-max_ir,max_ir
             if (rpts_array(nir1, nir2, nir3)==1) then
                iter=iter+1
-               irvec_newcell(:, iter)=[nir1, nir2, nir3]
+               irvec_surfacecell(:, iter)=[nir1, nir2, nir3]
                rpts_map(nir1, nir2, nir3)=iter
             endif
          enddo
       enddo
    enddo
 
-   !>> hamiltonian for the new cell
-   allocate(HmnR_newcell(Num_wann, Num_wann, Nrpts_new))
-   allocate(ndegen_newcell(Nrpts_new))
-   HmnR_newcell= 0d0
-   ndegen_newcell= 1
 
-   !> get number of R points for the new cell first 
+   allocate(HmnR_surfacecell(Num_wann, Num_wann, nrpts_surfacecell))
+   allocate(ndegen_surfacecell(nrpts_surfacecell))
+   HmnR_surfacecell= 0d0
+   ndegen_surfacecell= 1
+
+   !> 3. Allocate the old HmnR to the new HmnR.
+   !>  Note: The cell can't be treated as a mass point. We need to consider the relative coordinates of atoms.
    do ir=1, nrpts
       do j=1, Num_wann
          do i=1, Num_wann
-           !ia1=Origin_cell%spinorbital_to_atom_index(i)
-           !ia2=Origin_cell%spinorbital_to_atom_index(j)
-           !apos1d= Origin_cell%Atom_position_direct(:, ia1)- shift_vec_direct
-           !apos2d= Origin_cell%Atom_position_direct(:, ia2)- shift_vec_direct+irvec(:, ir)
-            apos1d= Origin_cell%wannier_centers_direct(:, i)- shift_vec_direct
-            apos2d= Origin_cell%wannier_centers_direct(:, j)- shift_vec_direct+irvec(:, ir)
-            call latticetransform(apos1d(1),apos1d(2),apos1d(3),new_ia,new_ib,new_ic)
-            irn1=floor([new_ia,new_ib,new_ic])
-
-            call latticetransform(apos2d(1),apos2d(2),apos2d(3),new_ia,new_ib,new_ic)
-            irn2=floor([new_ia,new_ib,new_ic])
-
-            nir1=irn2(1)-irn1(1)
-            nir2=irn2(2)-irn1(2)
-            nir3=irn2(3)-irn1(3)
-            if (abs(nir1)>max_ir .or. abs(nir2)>max_ir .or. abs(nir3)>max_ir) cycle
-            ir_cell= rpts_map(nir1, nir2, nir3)
-            HmnR_newcell(i,j,ir_cell)=HmnR(i,j,ir)/ndegen(ir)
+
+            ia1 = Origin_cell%spinorbital_to_atom_index(i)
+            ia2 = Origin_cell%spinorbital_to_atom_index(j)
+            apos1d = Origin_cell%Atom_position_direct(:, ia1)
+            apos2d = Origin_cell%Atom_position_direct(:, ia2) 
+
+            ia1prime = Cell_defined_by_surface%spinorbital_to_atom_index(i)
+            ia2prime = Cell_defined_by_surface%spinorbital_to_atom_index(j)
+            apos1dprime = Cell_defined_by_surface%Atom_position_direct(:,ia1prime)
+            apos2dprime = Cell_defined_by_surface%Atom_position_direct(:,ia2prime)
+
+            !> R'_mn = R' + tau'_2 - tau'_1
+            !> R_mn = R + tau_2 - tau_1
+            !> R'_mn = Pinv * R_mn
+            !> R' = (Pinv * R_mn) - (tau'_2 - tau'_1)
+
+            Rmn_old = irvec(:, ir) + apos2d - apos1d
+            call latticetransform(Rmn_old(1),Rmn_old(2),Rmn_old(3),Rmn_new(1),Rmn_new(2),Rmn_new(3))
+            irvec_new = Rmn_new - (apos2dprime - apos1dprime)
+
+            !> For safety, we perform a rounding operation due to the accuracy of computing
+            irvec_new_int(1) = ANINT(irvec_new(1))
+            irvec_new_int(2) = ANINT(irvec_new(2))
+            irvec_new_int(3) = ANINT(irvec_new(3))
+
+            if (abs(irvec_new_int(1))>max_ir .or. abs(irvec_new_int(2))>max_ir .or. abs(irvec_new_int(3))>max_ir) cycle
+            ir_cell = rpts_map(irvec_new_int(1), irvec_new_int(2), irvec_new_int(3))
+            HmnR_surfacecell(i,j,ir_cell) = HmnR(i,j,ir)/ndegen(ir)
+
          enddo
       enddo
    enddo
 
-   !> do cut-off according to the hopping value
-   iter= 0 
-   do ir=1, nrpts_new
-      max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree
-      if (max_val<eps4) then
-         iter= iter+ 1
-      endif
-   enddo
+   ! !> do cut-off according to the hopping value
+   ! iter= 0 
+   ! do ir=1, nrpts_new
+   !    max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree
+   !    if (max_val<eps4) then
+   !       iter= iter+ 1
+   !    endif
+   ! enddo
 
    if (cpuid==0.and.export_newhr) then
+      !> Problem: Just cut the nrpts_surfacecell, but the irvec is not well ordered. The order (iter) can be anywhere.
       !> write to new_hr.dat
-      nrpts_new=sum(rpts_array)-iter
+      ! nrpts_surfacecell=sum(rpts_array)-iter
+
+      nrpts_surfacecell=sum(rpts_array)
       outfileindex= outfileindex+ 1
       open(unit=outfileindex, file='wannier90_hr_newcell.dat')
       write(outfileindex, '(a,1X,a,1X,a,1X, a, a)')  &
          'HmnR for new cell generated at ', time_now, date_now, 'UTC', zone_now
       write(outfileindex, *)Num_wann
-      write(outfileindex, *)nrpts_new
-      write(outfileindex, '(15I5)')(1 , i=1, nrpts_new)
-      do ir=1, nrpts_new
-         max_val=maxval(abs(HmnR_newcell(:, :, ir)))/eV2Hartree
+      write(outfileindex, *)nrpts_surfacecell
+      write(outfileindex, '(15I5)')(1 , i=1, nrpts_surfacecell)
+      do ir=1, nrpts_surfacecell
+         max_val=maxval(abs(HmnR_surfacecell(:, :, ir)))/eV2Hartree
         !if (max_val<eps4) cycle
          do j=1, Num_wann
             do i=1, Num_wann
                write( outfileindex, '(5I5, 2f16.6)') &
-                  irvec_newcell(:, ir), i, j, HmnR_newcell(i, j, ir)/eV2Hartree
+                  irvec_surfacecell(:, ir), i, j, HmnR_surfacecell(i, j, ir)/eV2Hartree
             end do
          end do
       end do