Merge pull request #178 from zhangyan31415/support_nonorth

feat: add non-orthogonal basis support

Author: quanshengwu

src/ek_bulk.f90

@@ -21,7 +21,7 @@ subroutine ek_bulk_line
    real(Dp), allocatable :: W(:)
 
    ! Hamiltonian of bulk system
-   complex(Dp), allocatable :: Hamk_bulk(:, :)
+   complex(Dp), allocatable :: Hamk_bulk(:, :), Sk_bulk(:, :)
 
    ! eigenectors of H
    real(dp), allocatable :: eigv(:,:), eigv_mpi(:,:)
@@ -39,6 +39,10 @@ subroutine ek_bulk_line
    eigv    = 0d0; eigv_mpi= 0d0
    weight  = 0d0; weight_sum = 0d0; weight_mpi = 0d0
 
+   if (.not. Orthogonal_Basis) then
+      allocate(Sk_bulk(Num_wann, Num_wann))
+      Sk_bulk= 0d0
+   endif
 
    do ik= 1+cpuid, knv3, num_cpu
 
@@ -55,7 +59,20 @@ subroutine ek_bulk_line
          if (index(Particle,'phonon')/=0.and.LOTO_correction) then
             call ham_bulk_LOTO(k, Hamk_bulk)
          else
-            call ham_bulk_latticegauge(k, Hamk_bulk)
+            if (Is_Sparse) then
+               stop 'ERROR: Sparse mode is not supported in ek_bulk_line'
+            else
+               if (Orthogonal_Basis) then
+                  call ham_bulk_latticegauge(k, Hamk_bulk)
+               else
+                  
+                  Sk_bulk= 0d0
+                  call S_bulk_latticegauge(k, Sk_bulk)
+                  call ham_bulk_latticegauge(k, Hamk_bulk)
+                  call orthogonalize_hamiltonian(Hamk_bulk, Sk_bulk, Num_wann)
+               endif
+
+            endif
            !call ham_bulk_atomicgauge(k, Hamk_bulk)
          endif
       endif
@@ -3407,7 +3424,7 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
          endif
       enddo
       write(outfileindex, '(2a)')"# please comment the following lines to plot the fatband "
-      if (Is_Sparse_Hr) then
+      if (Is_Sparse) then
          write(outfileindex, '(4a)')"plot '", trim(adjustl(datafilename)), "' u 1:2 ",  &
             " w p pt 7  ps 0.2 lc rgb 'black', 0 w l lw 2 dt 2"
       else

src/ham_bulk.f90

@@ -402,21 +402,143 @@ subroutine ham_bulk_latticegauge(k,Hamk_bulk)
    !Hamk_bulk= (Hamk_bulk+ mat2)/2d0
 
    ! check hermitcity
-  !do i1=1, Num_wann
-  !   do i2=1, Num_wann
-  !      if(abs(Hamk_bulk(i1,i2)-conjg(Hamk_bulk(i2,i1))).ge.1e-6)then
-  !         write(stdout,*)'there is something wrong with Hamk_bulk'
-  !         write(stdout,*)'i1, i2', i1, i2
-  !         write(stdout,*)'value at (i1, i2)', Hamk_bulk(i1, i2)
-  !         write(stdout,*)'value at (i2, i1)', Hamk_bulk(i2, i1)
-  !         !stop
-  !      endif
-  !   enddo
-  !enddo
+  !打印出k点
+   ! write(*,*)'check Hamk_bulk hermitcity'
+  do i1=1, Num_wann
+    do i2=1, Num_wann
+       if(abs(Hamk_bulk(i1,i2)-conjg(Hamk_bulk(i2,i1))).ge.1e-6)then
+          write(stdout,*)'there is something wrong with Hamk_bulk'
+          write(stdout,*)'i1, i2', i1, i2
+          write(stdout,*)'value at (i1, i2)', Hamk_bulk(i1, i2)
+          write(stdout,*)'value at (i2, i1)', Hamk_bulk(i2, i1)
+          !stop
+       endif
+    enddo
+  enddo
 
    return
 end subroutine ham_bulk_latticegauge
 
+
+subroutine S_bulk_latticegauge(k,Sk_bulk)
+   ! This subroutine caculates Hamiltonian for
+   ! bulk system without the consideration of the atom's position
+   !
+   ! History
+   !
+   !        May/29/2011 by Quansheng Wu
+
+   use para, only : dp, pi2zi, SmnR, ndegen, nrpts, irvec, Num_wann, stdout, twopi, zi
+   implicit none
+
+   ! loop index
+   integer :: i1,i2,iR
+
+   real(dp) :: kdotr
+
+   complex(dp) :: factor
+
+   real(dp), intent(in) :: k(3)
+
+   ! Hamiltonian of bulk system
+   complex(Dp),intent(out) ::Sk_bulk(Num_wann, Num_wann)
+
+   Sk_bulk=0d0
+   do iR=1, Nrpts
+      kdotr= k(1)*irvec(1,iR) + k(2)*irvec(2,iR) + k(3)*irvec(3,iR)
+      factor= (cos(twopi*kdotr)+zi*sin(twopi*kdotr))
+
+      Sk_bulk(:, :)= Sk_bulk(:, :)+ SmnR(:, :, iR)*factor
+   enddo ! iR
+   
+   !call mat_mul(Num_wann, mirror_z, Hamk_bulk, mat1)
+   !call mat_mul(Num_wann, mat1, mirror_z, mat2)
+   !Hamk_bulk= (Hamk_bulk+ mat2)/2d0
+
+   ! check hermitcity
+   ! write(*,*)'check Sk_bulk'
+  do i1=1, Num_wann
+    do i2=1, Num_wann
+       if(abs(Sk_bulk(i1,i2)-conjg(Sk_bulk(i2,i1))).ge.1e-6)then
+          write(stdout,*)'there is something wrong with Sk_bulk'
+          write(stdout,*)'i1, i2', i1, i2
+          write(stdout,*)'value at (i1, i2)', Sk_bulk(i1, i2)
+          write(stdout,*)'value at (i2, i1)', Sk_bulk(i2, i1)
+          !stop
+       endif
+    enddo
+  enddo
+
+   return
+end subroutine S_bulk_latticegauge
+
+!---------------------------------------------------------------
+! Surroutines:change the non-orthogonal basis H and S to orthogonal basis
+!   H_new = U * H_old * U
+!   U = S_vec * M_inv * S_vec^H
+!---------------------------------------------------------------
+subroutine orthogonalize_hamiltonian(Hamk_bulk, Sk_bulk , N)
+   use para, only: Dp, stdout,E_fermi,eV2Hartree
+   implicit none
+ 
+   integer, intent(in)       :: N
+   complex(Dp), intent(in )    :: Sk_bulk(N,N)
+   complex(Dp), intent(inout)  :: Hamk_bulk(N,N)
+ 
+   !-- 局部 --
+   integer           :: i, info, j
+   real(Dp), allocatable    :: S_vals(:)
+   complex(Dp), allocatable :: S_vec(:,:), M_inv(:,:), U(:,:)
+ 
+
+   allocate(S_vals(N))
+   allocate(S_vec (N,N))
+   allocate(M_inv (N,N))
+   allocate(U     (N,N))
+   U = (0.0_dp,0.0_dp)
+   ! 1) Sk_bulk =  S_vec * diag(S_vals) * S_vec^H 
+   S_vec = Sk_bulk
+   call eigensystem_c('V','U', N, S_vec, S_vals)
+ 
+   ! 2)  M_inv = diag(1/√S_vals)
+   M_inv = (0.0_dp,0.0_dp)
+   ! do i = 1, N
+   !   M_inv(i,i) = 1.0_dp / sqrt(S_vals(i))
+   ! end do
+
+   do i = 1, N
+      if (S_vals(i) <= 0.0_dp) then
+         write(stdout,'(a,i5,1pe13.4)') ' WARNING: λ <= 0, i =', i, S_vals(i)
+      else
+         M_inv(i,i) = 1.0_dp / sqrt(S_vals(i))
+      endif
+   enddo
+ 
+   ! 3) U = S_vec * M_inv * S_vec^H
+   ! U = matmul( S_vec, matmul(M_inv, conjg(transpose(S_vec))) )
+   U = matmul( S_vec, M_inv )
+   U = matmul(U, conjg(transpose(S_vec)))
+
+   ! call mat_mul(N, S_vec, M_inv, U)
+   ! call mat_mul(N, U, conjg(transpose(S_vec)), U)
+
+ 
+   ! 4)Hamiltonian = U * Hamk_bulk * U
+   ! Hamk_bulk = matmul( U, matmul(Hamk_bulk, U) )
+   Hamk_bulk = matmul(U, Hamk_bulk)
+   Hamk_bulk = matmul(Hamk_bulk, U)
+   ! call mat_mul(N, U, Hamk_bulk, Hamk_bulk)
+   ! call mat_mul(N, Hamk_bulk, U, Hamk_bulk)
+
+   do i = 1, N
+      Hamk_bulk(i,i) = Hamk_bulk(i,i) - E_fermi*eV2Hartree
+   enddo
+
+   deallocate(S_vals, S_vec, M_inv, U)
+end subroutine orthogonalize_hamiltonian
+ 
+
+
 subroutine dHdk_latticegauge_wann(k, velocity_Wannier)
    use para, only : Nrpts, irvec, crvec, Origin_cell, &
       HmnR, ndegen, Num_wann, zi, pi2zi, dp, twopi

src/ham_qlayer2qlayer.f90

@@ -101,6 +101,110 @@ subroutine ham_qlayer2qlayer(k,H00new,H01new)
   return
   end subroutine ham_qlayer2qlayer
 
+  subroutine s_qlayer2qlayer(k,H00new,H01new)
+   ! This subroutine caculates Hamiltonian between
+   ! slabs For surface state calculation
+   ! 4/23/2010 by QS Wu
+   ! Copyright (c) 2010 QuanSheng Wu. All rights reserved.
+
+   use para
+
+   implicit none
+
+   ! loop index
+   integer :: i,j,iR
+
+   ! index used to sign irvec     
+   real(dp) :: ia,ib,ic
+
+   ! 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  
+
+   ! input wave vector k's cooridinates
+   real(Dp),intent(in) :: k(2)
+
+   ! H00 Hamiltonian between nearest neighbour-quintuple-layers
+   ! the factor 2 is induced by spin
+
+   complex(dp) :: ratio
+
+   complex(Dp), allocatable :: Hij(:, :, :)
+
+   ! H00 Hamiltonian between nearest neighbour-quintuple-layers
+   ! the factor 2 is induced by spin
+
+   !     complex(Dp),allocatable,intent(out) :: H00new(:,:)
+   complex(Dp),intent(out) :: H00new(Ndim,Ndim)
+
+   ! H01 Hamiltonian between next-nearest neighbour-quintuple-layers
+   !     complex(Dp),allocatable,intent(out) :: H01new(:,:)
+   complex(Dp),intent(out) :: H01new(Ndim,Ndim)
+
+   allocate(Hij(Num_wann,Num_wann,-ijmax:ijmax))
+
+   Hij=0.0d0
+   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, int_ic )&
+         =Hij(1:Num_wann, 1:Num_wann, int_ic)&
+         +SmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(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
+
+   enddo
+   enddo
+
+   deallocate(Hij)
+
+ return
+ end subroutine s_qlayer2qlayer
+
+
   subroutine ham_qlayer2qlayer_bak(k,H00new,H01new)
      ! This subroutine caculates Hamiltonian between
      ! slabs For surface state calculation
@@ -481,7 +585,55 @@ subroutine ham_qlayer2qlayer2(k,Hij)
   return
   end subroutine ham_qlayer2qlayer2
 
+  subroutine s_qlayer2qlayer2(k,Hij)
+   ! This subroutine caculates Hamiltonian between
+   ! slabs  
+   ! 4/23/2010 by QS Wu
+   ! Copyright (c) 2010 QuanSheng Wu. All rights reserved.
+
+   use para
+
+   implicit none
+
+   ! loop index
+   integer :: iR, inew_ic, int_ic
+   real(dp) :: ia, ib, ic
+
+   ! new index used to sign irvec     
+   real(dp) :: new_ia,new_ib,new_ic
+
+   ! wave vector k times lattice vector R  
+   real(Dp) :: kdotr  
+
+   ! input wave vector k's cooridinates
+   real(Dp),intent(in) :: k(2)
+
+   complex(dp) :: ratio
+
+   complex(Dp), intent(out) :: Hij(-ijmax:ijmax,Num_wann,Num_wann)
+
+   Hij=0.0d0
+
+   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(int_ic, 1:Num_wann, 1:Num_wann )&
+         =Hij(int_ic, 1:Num_wann, 1:Num_wann )&
+         +SmnR_surfacecell(:,:,iR)*ratio/ndegen_surfacecell(iR)
+
+      endif
+
+   enddo
 
+return
+end subroutine s_qlayer2qlayer2
 
   subroutine ham_qlayer2qlayer2_LOTO(k,Hij)
      ! This subroutine caculates Hamiltonian between