try to add valley operator, not tested yet

src/ek_bulk.f90

@@ -236,7 +236,8 @@ subroutine ek_bulk_line_valley
       eigv(:, ik)= W
       do j=1, Num_wann  !> band
          psi= Hamk_bulk(:, j)
-         call zgemv('N', Num_wann, Num_wann, One_complex, psi, Num_wann, valley_k, 1, zzero, vpsi, 1)
+         vpsi=0d0
+         call zgemv('N', Num_wann, Num_wann, One_complex,  valley_k, Num_wann, psi, 1, zzero, vpsi, 1)
          weight(j, ik)= real(zdotc(Num_wann, psi, 1, vpsi, 1))
       enddo ! i
    enddo ! ik
@@ -258,7 +259,7 @@ subroutine ek_bulk_line_valley
 
       do i=1, Num_wann
          do ik=1, knv3
-            write(outfileindex, '(200E16.5)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), &
+            write(outfileindex, '(200F16.8)')k3len(ik)*Angstrom2atomic,eigv_mpi(i, ik), &
                weight_mpi(i, ik)
          enddo
          write(outfileindex, *)' '
@@ -1354,8 +1355,7 @@ subroutine sparse_ekbulk_valley
          psi(:)= zeigv(:, ib)  !> the eigenvector of ib'th band
 
          !> weight_valley= <psi|vz|psi>
-         call mkl_zcoomv('N', Num_wann, Num_wann, One_complex, matdescra, acoo_valley, &
-            icoo_valley, jcoo_valley, nnzmax_valley, psi, zzero, vpsi)
+         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))
 

src/ham_bulk.f90

@@ -94,16 +94,9 @@ end subroutine ham_bulk_atomicgauge
 
 
 subroutine valley_k_atomicgauge(k,valley_k)
-   ! This subroutine caculates Hamiltonian for
-   ! bulk system with the consideration of the atom's position
-   !
+   ! This subroutine performs the Fourier transform of avalley operator
    ! History
-   !
-   !        May/29/2011 by Quansheng Wu
-   !  Atomic gauge Guan Yifei 2019
-   !  Lattice gauge Hl
-   !  Atomic gauge Ha= U* Hl U 
-   !  where U = e^ik.wc(i) on diagonal
+   !        Nov/5/2023 by Quansheng Wu
 
    use para
    implicit none
@@ -146,6 +139,20 @@ subroutine valley_k_atomicgauge(k,valley_k)
       enddo ! i2
    enddo ! iR
 
+   ! check hermitcity
+   do i1=1, Num_wann
+      do i2=1, Num_wann
+         if(abs(valley_k(i1,i2)-conjg(valley_k(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)', valley_k(i1, i2)
+            write(stdout,*)'value at (i2, i1)', valley_k(i2, i1)
+            !stop
+         endif
+      enddo
+   enddo
+
+
    return
 end subroutine valley_k_atomicgauge
 
@@ -819,14 +826,14 @@ subroutine valley_k_coo_sparsehr(nnz, k,acoo,icoo,jcoo)
    integer :: i,j,ir
 
    do i=1,nnz
-      ir=valley_operator_irv(i)
+      ir= valley_operator_irv(i)
       icoo(i)=valley_operator_icoo(i)
       jcoo(i)=valley_operator_jcoo(i)
       posij=irvec_valley(:, ir)+ Origin_cell%wannier_centers_direct(:, jcoo(i))- Origin_cell%wannier_centers_direct(:, icoo(i))
       kdotr=posij(1)*k(1)+posij(2)*k(2)+posij(3)*k(3)
       ratio= (cos(twopi*kdotr)+zi*sin(twopi*kdotr))
       acoo(i)=ratio*valley_operator_acoo(i)
-   end do
+   enddo
 
    return
 end subroutine valley_k_coo_sparsehr

src/readHmnR.f90

@@ -237,16 +237,11 @@ subroutine readNormalHmnR()
 end subroutine readNormalHmnR
 
 subroutine read_valley_operator()
-   !>> Read in the tight-binding model from wannier90_hr.dat
-   !>  The format is defined by the wannier90 software
-   ! Constructed by quansheng wu 4/2/2010
-   !
-   ! Yifei Guan added the sparse hr file parsing June/2018
+   !>> Read in the valley operator from valley_operator.dat
+   ! Constructed by quansheng wu 04 Nov. 2023
    ! License: GPL V3
 
    use para
-   !> in: N of wann
-   !> out : nth atom
 
    implicit none
 
@@ -257,14 +252,20 @@ subroutine read_valley_operator()
    integer :: n, m, ir0
    integer :: add_electric_field
    integer :: nwann, nwann_nsoc
+   logical :: exists
 
    real(dp) :: static_potential
    real(dp) :: tot, rh, ih
    real(dp) :: pos(Origin_cell%Num_atoms)
 
 
    if(cpuid.eq.0)write(stdout,*)' '
-   open(12, file="valley_operator.dat", status='OLD')
+   inquire (file ="valley_operator.dat", EXIST = exists)
+   if (exists)then
+      open(12, file="valley_operator.dat", status='OLD')
+   else
+      STOP ">> for valley projection , you have to prepare a file valley_operator.dat"
+   endif
 
    !> skip a comment line
    read(12, *)
@@ -727,13 +728,21 @@ subroutine readsparse_valley_operator
    implicit none
    integer:: i,j,nwann,nwann_nsoc,i1,i2,i3,i4,i5,ir,n,m 
    real(dp) :: r1, r2
+   logical :: exists
 
    !> the direction which adding electric field which is also the stacking direction
    integer :: add_electric_field
    real(dp) :: Bx_in_au, By_in_au, Bz_in_au
    complex(dp) :: h_value
 
-   open(13, file='valley_operator.dat')
+   if(cpuid.eq.0)write(stdout,*)' '
+   inquire (file ="valley_operator.dat", EXIST = exists)
+   if (exists)then
+      open(13, file="valley_operator.dat", status='OLD')
+   else
+      STOP ">> for valley projection , you have to prepare a file valley_operator.dat"
+   endif
+
 
    !> skip a comment line
    read(13, *)
@@ -797,10 +806,12 @@ subroutine readsparse_valley_operator
       enddo
    enddo
 1002 continue
+   close(13)
    !> correct nrpts
    Nrpts_valley=ir
 
    if (cpuid.eq.0) write(stdout, '(a, i6)')' >> Nrpts_valley is ', Nrpts_valley
+   if (cpuid.eq.0) write(stdout, '(a)')' >> splen_valley_input', splen_valley_input, j
    if (cpuid.eq.0) write(stdout, '(a)')' >> Valley operator reading finished '
 
    return