diff --git a/src/ek_bulk.f90 b/src/ek_bulk.f90 index 61887189..e12e2994 100644 --- a/src/ek_bulk.f90 +++ b/src/ek_bulk.f90 @@ -210,13 +210,13 @@ subroutine ek_bulk_line_valley knv3= nk3_band allocate(W(Num_wann)) allocate(Hamk_bulk(Num_wann, Num_wann)) - allocate(valley_k(Num_wann, Num_wann)) allocate( eigv (Num_wann, knv3)) allocate( eigv_mpi(Num_wann, knv3)) allocate( weight (Num_wann, knv3)) allocate( weight_mpi(Num_wann, knv3)) allocate( psi(Num_wann), vpsi(Num_wann)) allocate( psi1(Num_wann), psi2(Num_wann)) + allocate(valley_k(Num_wann, Num_wann)) allocate(valley_plus(Num_wann, knv3)) allocate(valley_plus_mpi(Num_wann, knv3)) allocate( valley_k_nd(NDMAX, NDMAX), valley_eig(NDMAX)) @@ -251,9 +251,12 @@ subroutine ek_bulk_line_valley IE=1 do while (ie.le.Num_wann) ND=1 - do while ((ie+ND).le.Num_wann .and. (W(ie+ND)- W(ie)).lt.tolde) - ND= ND+1 - enddo + if (ie+ND.le.Num_wann) then + do while ((ie+ND).le.Num_wann .and. (W(ie+ND)- W(ie)).lt.tolde) + if (ie+ND .ge. Num_wann) exit + ND= ND+1 + enddo + endif !> if the degeneracy is larger than 1, we need to calculate the matrix valley_k_nd= 0d0 @@ -302,13 +305,14 @@ subroutine ek_bulk_line_valley enddo ! ik + #if defined (MPI) call mpi_allreduce(eigv,eigv_mpi,size(eigv),& mpi_dp,mpi_sum,mpi_cmw,ierr) call mpi_allreduce(weight, weight_mpi,size(weight),& mpi_dp,mpi_sum,mpi_cmw,ierr) call mpi_allreduce(valley_plus, valley_plus_mpi,size(valley_plus),& - mpi_logical,mpi_land,mpi_cmw,ierr) + mpi_logical,mpi_lor,mpi_cmw,ierr) #else eigv_mpi= eigv weight_mpi= weight @@ -1319,6 +1323,7 @@ subroutine sparse_ekbulk_valley !> some temporary integers integer :: ik, i, j, ierr, ib, ig + integer :: ie, ND, NDMAX, ie1, ie2 ! wave vector real(dp) :: k3(3) @@ -1339,10 +1344,16 @@ subroutine sparse_ekbulk_valley integer, allocatable :: icoo(:), icoo_valley(:) !> eigenvector of the sparse matrix acoo. Dim=(Num_wann, neval) - complex(dp), allocatable :: psi(:), vpsi(:) complex(dp), allocatable :: psi_project(:) complex(dp), allocatable :: zeigv(:, :) + ! Hamiltonian of bulk system + complex(Dp), allocatable :: psi(:), vpsi(:), valley_k_nd(:, :), psi1(:), psi2(:) + complex(Dp), allocatable :: valley_k(:, :) + complex(Dp), allocatable :: VL(:, :), VR(:, :), valley_eig(:) + logical, allocatable :: valley_plus(:, :), valley_plus_mpi(:, :) + + !> print the weight for the Selected_WannierOrbitals real(dp), allocatable :: weight_valley(:, :), weight_valley_mpi(:, :) @@ -1362,10 +1373,10 @@ subroutine sparse_ekbulk_valley real(dp) :: time1, time2, time3 !> - character :: matdescra(6) ! vector multiply vector v1*v2 complex(dp), external :: zdotc + real(dp) :: tolde= 1E-4 !if (OmegaNum==0) OmegaNum= Num_wann @@ -1393,6 +1404,7 @@ subroutine sparse_ekbulk_valley sigma=(1d0,0d0)*E_arc nnzmax=splen+Num_wann nnz=splen + NDMAX= 12 allocate( acoo(nnzmax)) allocate( jcoo(nnzmax)) allocate( icoo(nnzmax)) @@ -1405,10 +1417,18 @@ subroutine sparse_ekbulk_valley allocate( W( neval)) allocate( eigv( neval, nk3_band)) allocate( eigv_mpi( neval, nk3_band)) - allocate( psi(Num_wann), vpsi(Num_wann)) allocate( psi_project(Num_wann)) allocate( zeigv(Num_wann,nvecs)) allocate( weight_valley(neval, nk3_band), weight_valley_mpi(neval, nk3_band)) + allocate( psi(Num_wann), vpsi(Num_wann)) + allocate( psi1(Num_wann), psi2(Num_wann)) + allocate(valley_k(Num_wann, Num_wann)) + allocate(valley_plus(Num_wann, nk3_band)) + allocate(valley_plus_mpi(Num_wann, nk3_band)) + allocate( valley_k_nd(NDMAX, NDMAX), valley_eig(NDMAX)) + allocate( VL(NDMAX, NDMAX), VR(NDMAX, NDMAX)) + valley_k_nd= 0d0 + valley_plus= .False.; valley_plus_mpi= .False. psi=0d0; vpsi=0d0; psi_project= 0d0; zeigv= 0d0 weight_valley_mpi=0d0; weight_valley=0d0 @@ -1417,11 +1437,6 @@ subroutine sparse_ekbulk_valley ritzvec= .True. - matdescra(1)= 'H' - matdescra(2)= 'U' - matdescra(3)= 'U' - matdescra(4)= 'F' - !> calculate the energy bands along special k line k3= 0 do ik=1+ cpuid, nk3_band, num_cpu @@ -1444,16 +1459,75 @@ subroutine sparse_ekbulk_valley call valley_k_coo_sparsehr(nnzmax_valley, k3, acoo_valley, icoo_valley, jcoo_valley) !> get the valley projection - do ib= 1, neval - psi(:)= zeigv(:, ib) !> the eigenvector of ib'th band +! do ib= 1, neval +! psi(:)= zeigv(:, ib) !> the eigenvector of ib'th band - !> weight_valley= - call mkl_zcoogemv('N', Num_wann, acoo_valley, icoo_valley, jcoo_valley, nnzmax_valley, psi, vpsi) +! !> weight_valley= +! call mkl_zcoogemv('N', Num_wann, acoo_valley, icoo_valley, jcoo_valley, nnzmax_valley, psi, vpsi) + +! weight_valley(ib, ik)= real(zdotc(Num_wann, psi, 1, vpsi, 1)) + +! enddo + + !> check the degeneracy of each band + IE=1 + do while (ie.le.neval) + ND=1 + if (ie+ND.le.neval) then + do while ((ie+ND).le.neval .and. (W(ie+ND)- W(ie)).lt.tolde) + if (ie+ND .ge. neval) exit + ND= ND+1 + enddo + endif + + + !> if the degeneracy is larger than 1, we need to calculate the matrix + valley_k_nd= 0d0 + do ie1= 1, ND + psi1= zeigv(:, IE+ie1-1) + + do ie2= 1, ND + psi2= zeigv(:, IE+ie2-1) + vpsi=0d0 + call mkl_zcoogemv('N', Num_wann, acoo_valley, icoo_valley, jcoo_valley, nnzmax_valley, psi2, vpsi) + valley_k_nd(ie1, ie2)= zdotc(Num_wann, psi1, 1, vpsi, 1) + enddo + !if (abs(W(IE+ie1-1)/eV2Hartree)<2d0) & + !write( *, '(a, 2i5, 200f8.3)') "ND, ie1: ", ND, ie1, valley_k_nd(ie1, 1:ND) + enddo + + !> diagonalize the matrix + VL = 0d0; VR= 0d0 + if (ND>1) then + call zgeev_sys(ND, valley_k_nd(1:ND, 1:ND), valley_eig(1:ND),'N',VL(1:ND, 1:ND),"V",VR(1:ND, 1:ND) ) + do ie1= 1, ND + weight_valley(IE+ie1-1, ik) = real(valley_eig(ie1)) + enddo + else + valley_eig(1)= valley_k_nd(1, 1) + weight_valley(IE, ik) = real(valley_k_nd(1, 1)) + endif + ! if (abs(W(IE+ie1-1)/eV2Hartree)<2d0) & + ! write(*, '(a, 2i5, a, 200f8.3)'), 'ND, ie', ND, ie, ' valley_eig', real(valley_eig(1:ND)) + ! if (abs(W(IE+ie1-1)/eV2Hartree)<2d0) & + ! print *, ' ' - weight_valley(ib, ik)= real(zdotc(Num_wann, psi, 1, vpsi, 1)) + !> get new eigenvector + !do ie1= 1, ND + ! psi= 0d0 + ! do ie2= 1, ND + ! psi= psi+ VR(ie2, ie1)* zeigv(:, IE+ie2-1) + ! enddo + !enddo + !> if the degeneracy is larger than 1, we need to calculate the matrix + IE= IE+ ND + enddo + do ie=1, neval + if (weight_valley(ie, ik)>0) valley_plus(ie, ik)=.true. enddo + if (cpuid==0)write(stdout, '(a, f20.2, a)')' >> Time cost for constructing H: ', time2-time1, ' s' if (cpuid==0)write(stdout, '(a, f20.2, a)')' >> Time cost for diagonalize H: ', time3-time2, ' s' enddo !ik @@ -1463,9 +1537,12 @@ subroutine sparse_ekbulk_valley mpi_dp,mpi_sum,mpi_cmw,ierr) call mpi_allreduce(weight_valley, weight_valley_mpi,size(weight_valley),& mpi_dp,mpi_sum,mpi_cmw,ierr) + call mpi_allreduce(valley_plus, valley_plus_mpi,size(valley_plus),& + mpi_logical,mpi_lor,mpi_cmw,ierr) #else eigv_mpi= eigv weight_valley_mpi= weight_valley + valley_plus_mpi= valley_plus #endif !> minimum and maximum value of energy bands @@ -1487,6 +1564,37 @@ subroutine sparse_ekbulk_valley close(outfileindex) endif + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek_valley_plus.dat') + do i=1, neval + do ik=1, nk3_band + if (valley_plus_mpi(i, ik)) then + write(outfileindex, '(200F16.8)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), & + weight_valley_mpi(i, ik) + endif + enddo + write(outfileindex, *)' ' + enddo + close(outfileindex) + endif + + outfileindex= outfileindex+ 1 + if (cpuid==0)then + open(unit=outfileindex, file='bulkek_valley_minus.dat') + do i=1, neval + do ik=1, nk3_band + if (.not.valley_plus_mpi(i, ik)) then + write(outfileindex, '(200F16.8)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), & + weight_valley_mpi(i, ik) + endif + enddo + write(outfileindex, *)' ' + enddo + close(outfileindex) + endif + + call generate_ek_kpath_gnu('bulkek.dat', 'bulkek.gnu', 'bulkek.pdf', & emin, emax, nk3_band, Nk3lines, & k3line_name, k3line_stop, k3len)