change center_atom_for_electric_field from one atom to two atoms

src/ek_bulk.f90

@@ -187,6 +187,7 @@ subroutine ek_bulk_line_valley
    implicit none
 
    integer :: ik, il, ig, io, i, j, knv3, ierr
+   integer :: ie, ND, NDMAX
    real(dp) :: emin,  emax,  k(3)
    character*40 :: filename
    complex(dp), external :: zdotc
@@ -195,12 +196,13 @@ subroutine ek_bulk_line_valley
    real(Dp), allocatable :: W(:)
 
    ! Hamiltonian of bulk system
-   complex(Dp), allocatable :: psi(:), vpsi(:)
+   complex(Dp), allocatable :: psi(:), vpsi(:), valley_k_nd(:, :)
    complex(Dp), allocatable :: Hamk_bulk(:, :), valley_k(:, :)
 
    ! eigenectors of H
    real(dp), allocatable :: eigv(:,:), eigv_mpi(:,:)
    real(dp), allocatable :: weight(:,:), weight_mpi(:,:)
+   real(dp) :: tolde= 1E-6
 
    knv3= nk3_band
    allocate(W(Num_wann))
@@ -240,6 +242,20 @@ subroutine ek_bulk_line_valley
          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
+
+      !> check the degeneracy of each band
+      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 the degeneracy is larger than 1, we need to calculate the matrix
+         IE= IE+ ND
+
+      enddo
+
    enddo ! ik
 
 #if defined (MPI)
@@ -256,6 +272,10 @@ subroutine ek_bulk_line_valley
    outfileindex= outfileindex+ 1
    if (cpuid==0)then
       open(unit=outfileindex, file='bulkek.dat')
+      open(unit=outfileindex, file='bulkek_valley_plus.dat')
+      open(unit=outfileindex, file='bulkek_valley_minus.dat')
+
+
 
       do i=1, Num_wann
          do ik=1, knv3

src/module.f90

@@ -611,7 +611,7 @@ module para
      !> Electric field along the stacking direction of a 2D system in eV/Angstrom
      real(dp) :: Electric_field_in_eVpA
      real(dp) :: Symmetrical_Electric_field_in_eVpA
-     integer :: center_atom_for_electric_field  ! At this atom, the electric potential is zero
+     integer :: center_atom_for_electric_field(2)  ! At this atom, the electric potential is zero
      logical :: Inner_symmetrical_Electric_Field
 
      !> a parameter to control the Vacumm thickness for the slab system 

src/readHmnR.f90

@@ -895,8 +895,9 @@ subroutine get_stacking_direction_and_pos(add_electric_field, pos)
    if (add_electric_field>0) then
       pos=Origin_cell%Atom_position_direct(add_electric_field, :)
       pos= mod(pos, 1d0)-0.5d0
-      if (center_atom_for_electric_field>0) then
-         center=  Origin_cell%Atom_position_direct(add_electric_field, center_atom_for_electric_field)
+      if (sum(center_atom_for_electric_field)>0) then
+         center= (Origin_cell%Atom_position_direct(add_electric_field, center_atom_for_electric_field(1)) &
+                + Origin_cell%Atom_position_direct(add_electric_field, center_atom_for_electric_field(2)))/2d0
          center= mod(center, 1d0)-0.5d0
       else
          center= (maxval(pos)+minval(pos))/2d0