Skip to content

Commit

Permalink
update the calculation of valley chern number version 202311101224
Browse files Browse the repository at this point in the history
  • Loading branch information
@quanshengwu
quanshengwu committed Nov 9, 2023
1 parent e228387 commit 5722d75
Showing 1 changed file with 126 additions and 18 deletions.
144 changes: 126 additions & 18 deletions src/ek_bulk.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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(:, :)

Expand All @@ -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
Expand Down Expand Up @@ -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))
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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= <psi|vz|psi>
call mkl_zcoogemv('N', Num_wann, acoo_valley, icoo_valley, jcoo_valley, nnzmax_valley, psi, vpsi)
! !> weight_valley= <psi|vz|psi>
! 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
Expand All @@ -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
Expand All @@ -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)
Expand Down

0 comments on commit 5722d75

Please sign in to comment.