diff --git a/examples/Haldane_model/wt.in-BerryCurvature_kpath_sepband_calc b/examples/Haldane_model/wt.in-BerryCurvature_kpath_sepband_calc
new file mode 100644
index 0000000..b1cb7d6
--- /dev/null
+++ b/examples/Haldane_model/wt.in-BerryCurvature_kpath_sepband_calc
@@ -0,0 +1,72 @@
+&TB_FILE
+Hrfile = "Haldane_hr.dat"
+/
+
+
+!> bulk band structure calculation flag
+&CONTROL
+BerryCurvature_kpath_sepband_calc   = T  ! calculate BC sum over 1~NumOccupied bands
+/
+
+&SYSTEM
+NSLAB =60  
+NumOccupied = 1         ! NumOccupied
+SOC = 0                 ! soc
+E_FERMI = 0        ! e-fermi
+/
+
+&PARAMETERS
+Eta_Arc = 0.01     ! infinite small value, like brodening 
+E_arc = 0.0         ! energy for calculate Fermi Arc
+OmegaNum = 1000  ! omega number       
+OmegaMin = -5.0     ! energy interval
+OmegaMax =  5.0     ! energy interval
+Nk1 = 60            ! number k points 
+Nk2 = 60            ! number k points 
+NP = 1              ! number of principle layers
+/
+
+LATTICE
+Angstrom
+2.1377110  -1.2342080   0.0000000
+0.0000000   2.4684160   0.0000000
+0.0000000   0.0000000   10.000000
+
+ATOM_POSITIONS
+2                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+C 0.333333 0.666667 0.500000 
+C 0.666667 0.333333 0.500000 
+
+PROJECTORS
+1 1          ! number of projectors
+C pz
+C pz
+
+
+SURFACE            ! See doc for details
+ 0  0  1
+ 1  0  0
+ 0  1  0
+
+KPATH_BULK            ! k point path
+4              ! number of k line only for bulk band
+  M   0.50000  0.00000  0.00000   K' -.33333   -.33333  0.00000
+  K'  -.33333  -.33333  0.00000   G  0.00000   0.00000  0.00000
+  G   0.00000  0.00000  0.00000   K  0.33333   0.33333  0.00000
+  K   0.33333  0.33333  0.00000   M  0.50000   0.00000  0.00000
+
+KPATH_SLAB
+1        ! numker of k line for 2D case
+0 0.0 0.0 1 0. 1.0  ! k path for 2D case
+
+KPLANE_SLAB
+-0.5 -0.5      ! Original point for 2D k plane
+ 1.0  0.0      ! The first vector to define 2D k plane 
+ 0.0  1.0      ! The second vector to define 2D k plane  for arc plots
+
+KPLANE_BULK
+ 0.00  0.00  0.00   ! Original point for 3D k plane 
+ 1.00  0.00  0.00   ! The first vector to define 3d k space plane
+ 0.00  1.00  0.00   ! The second vector to define 3d k space plane
+
diff --git a/src/Boltz_transport_anomalous.f90 b/src/Boltz_transport_anomalous.f90
index f01ca79..0829974 100644
--- a/src/Boltz_transport_anomalous.f90
+++ b/src/Boltz_transport_anomalous.f90
@@ -425,7 +425,7 @@ subroutine sigma_ahc_vary_ChemicalPotential(NumOfmu, mulist, NumberofEta, eta_ar
 
         !> calculate Berry curvature at a single k point for all bands
         !> \Omega_n^{\gamma}(k)=i\sum_{\alpha\beta}\epsilon_{\gamma\alpha\beta}(D^{\alpha\dag}D^{\beta})_{nn}
-        call berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+        call Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
  
         do ieta= 1, NumberofEta
            eta_local = eta_array(ieta)
diff --git a/src/berrycurvature.f90 b/src/berrycurvature.f90
index 7fbe42e..4ea46a4 100644
--- a/src/berrycurvature.f90
+++ b/src/berrycurvature.f90
@@ -1,4 +1,4 @@
-  subroutine berry_curvarture_singlek_EF(k, mu, Omega_x, Omega_y, Omega_z)
+  subroutine Berry_curvature_singlek_EF(k, mu, Omega_x, Omega_y, Omega_z)
      !> Calculate Berry curvature for a sigle k point
      !> The Fermi distribution is determined by the Fermi level iso_energy
      !> ref : Physical Review B 74, 195118(2006)
@@ -78,9 +78,9 @@ subroutine berry_curvarture_singlek_EF(k, mu, Omega_x, Omega_y, Omega_z)
      enddo
 
      return
-  end subroutine berry_curvarture_singlek_EF
+  end subroutine Berry_curvature_singlek_EF
 
-  subroutine berry_curvarture_singlek_numoccupied_slab_total(k, Omega_z)
+  subroutine Berry_curvature_singlek_numoccupied_slab_total(k, Omega_z)
      !> Calculate Berry curvature for a sigle k point
      !> The Fermi distribution is determined by the NumOccupied bands, not the Fermi level
      !> ref : Physical Review B 74, 195118(2006)
@@ -158,9 +158,9 @@ subroutine berry_curvarture_singlek_numoccupied_slab_total(k, Omega_z)
      Omega_z= -aimag(Omega_z*2d0)/Angstrom2atomic**2
 
      return
-  end subroutine berry_curvarture_singlek_numoccupied_slab_total
+  end subroutine Berry_curvature_singlek_numoccupied_slab_total
 
-  subroutine berry_curvarture_singlek_numoccupied_total(k, Omega_x, Omega_y, Omega_z)
+  subroutine Berry_curvature_singlek_numoccupied_total(k, Omega_x, Omega_y, Omega_z)
      !> Calculate Berry curvature for a sigle k point using Kubo-formula
      !> The Fermi distribution is determined by the NumOccupied bands, not the Fermi level
      !> ref : Physical Review B 74, 195118(2006)
@@ -230,9 +230,9 @@ subroutine berry_curvarture_singlek_numoccupied_total(k, Omega_x, Omega_y, Omega
      deallocate(Amat, velocity_wann, velocity_Ham)
 
      return
-  end subroutine berry_curvarture_singlek_numoccupied_total
+  end subroutine Berry_curvature_singlek_numoccupied_total
 
-  subroutine berry_curvarture_singlek_numoccupied(k, Omega_x, Omega_y, Omega_z)
+  subroutine Berry_curvature_singlek_numoccupied(k, Omega_x, Omega_y, Omega_z)
      !> Calculate Berry curvature for a sigle k point
      !> The Fermi distribution is determined by the NumOccupied bands, not the Fermi level
      !> ref : Physical Review B 74, 195118(2006)
@@ -300,9 +300,9 @@ subroutine berry_curvarture_singlek_numoccupied(k, Omega_x, Omega_y, Omega_z)
      Omega_z= -Omega_z*2d0*zi
 
      return
-  end subroutine berry_curvarture_singlek_numoccupied
+  end subroutine Berry_curvature_singlek_numoccupied
 
-  subroutine berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+  subroutine Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
      !> Calculate Berry curvature for a sigle k point and all bands
      !> ref : eqn (30) Physical Review B 74, 195118(2006)
      !> \Omega_n^{\gamma}(k)=i\sum_{\alpha\beta}\epsilon_{\gamma\alpha\beta}(D^{\alpha\dag}D^{\beta})_{nn}
@@ -360,7 +360,7 @@ subroutine berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
      enddo
 
      return
-  end subroutine berry_curvarture_singlek_allbands
+  end subroutine Berry_curvature_singlek_allbands
 
 
   subroutine orbital_magnetization_singlek_allbands(Dmn_Ham, Vmn_Ham_nondiag, m_OrbMag)
@@ -462,7 +462,7 @@ subroutine get_Dmn_Ham(W, velocity_Ham, Dmn_Ham)
      return
   end subroutine get_Dmn_Ham
 
-  subroutine berry_curvarture_slab
+  subroutine Berry_curvature_slab
      !> Calculate Berry curvature 
      !
      !> ref : Physical Review B 74, 195118(2006)
@@ -521,7 +521,7 @@ subroutine berry_curvarture_slab
 
         Omega_z= 0d0
 
-        call berry_curvarture_singlek_numoccupied_slab_total(k, Omega_z(1))
+        call Berry_curvature_singlek_numoccupied_slab_total(k, Omega_z(1))
         Omega(ik) = sum(Omega_z)
 
      enddo ! ik
@@ -599,9 +599,9 @@ subroutine berry_curvarture_slab
  
      return
 
-  end subroutine berry_curvarture_slab
+  end subroutine Berry_curvature_slab
 
-  subroutine berry_curvarture_line
+  subroutine Berry_curvature_line_occupied
      !> Calculate Berry curvature for a k line defined bu KPATH_BULK
      !
      !> ref : Physical Review B 74, 195118(2006)
@@ -623,10 +623,13 @@ subroutine berry_curvarture_line
 
      real(dp), external :: norm
 
-     !> Berry curvature  (3, bands, k)
      complex(dp), allocatable :: Omega_x(:), Omega_y(:), Omega_z(:)
+     !> Berry curvature  (3, k)
      real(dp), allocatable :: Omega(:, :), Omega_mpi(:, :)
 
+     !> Berry curvature  (3, bands, k)
+     real(dp), allocatable :: Omega_sep_bk(:, :), Omega_sep_bk_mpi(:, :)
+
      !> energy bands
      real(dp), allocatable :: eigv(:,:)
      real(dp), allocatable :: eigv_mpi(:,:)
@@ -654,7 +657,7 @@ subroutine berry_curvarture_line
            ' time elapsed: ', time_end-time_start0 
 
         !> a k point in fractional coordinates
-        k= k3points(:, ik)
+        k= kpath_3d(:, ik)
 
         call now(time_start)
  
@@ -662,13 +665,13 @@ subroutine berry_curvarture_line
         Omega_y= 0d0
         Omega_z= 0d0
 
-        !call berry_curvarture_singlek_numoccupied_old(k, Omega_x, Omega_y, Omega_z)
+        !call Berry_curvature_singlek_numoccupied_old(k, Omega_x, Omega_y, Omega_z)
         if (Berrycurvature_kpath_EF_calc) then
-           call berry_curvarture_singlek_EF(k, iso_energy, Omega_x, Omega_y, Omega_z)
+           call Berry_curvature_singlek_EF(k, iso_energy, Omega_x, Omega_y, Omega_z)
         else if (BerryCurvature_kpath_Occupied_calc) then
-           call berry_curvarture_singlek_numoccupied_total(k, Omega_x(1), Omega_y(1), Omega_z(1))
+           call Berry_curvature_singlek_numoccupied_total(k, Omega_x(1), Omega_y(1), Omega_z(1))
         else
-           write(*, *) 'ERROR: In subroutine berry_curvarture_line, we only support BerryCurvature_kpath_Occupied_calc '
+           write(*, *) 'ERROR: In subroutine Berry_curvature_line, we only support BerryCurvature_kpath_Occupied_calc '
            write(*, *) ' and Berrycurvature_kpath_EF_calc '
            stop
         endif
@@ -698,7 +701,7 @@ subroutine berry_curvarture_line
            'real(Omega_x)', 'real(Omega_y)', 'real(Omega_z)'
 
         do ik= 1, nk3_band
-           k=k3points(:, ik)
+           k=kpath_3d(:, ik)
            write(outfileindex, '(20E18.8)')k3len(ik)*Angstrom2atomic, real(Omega_mpi(:, ik))/Angstrom2atomic**2
         enddo
 
@@ -757,10 +760,10 @@ subroutine berry_curvarture_line
  
      return
 
-  end subroutine berry_curvarture_line
+  end subroutine Berry_curvature_line_occupied
 
 
-  subroutine berry_curvarture_cube
+  subroutine berry_curvature_cube
      !> Calculate Berry curvature  in a cube defined in KCUBE_BULK
      !
      !> ref : Physical Review B 74, 195118(2006)
@@ -775,9 +778,9 @@ subroutine berry_curvarture_cube
      use para
      implicit none
     
-     integer :: ik, ierr, ikx, iky, ikz, knv3, i, m, n
+     integer :: ik, ierr, ikx, iky, ikz, n_kpoints, i, m, n
 
-     real(dp) :: k(3), o1(3), k_cart(3)
+     real(dp) :: k(3), o1(3), k_cart(3), emin, emax
      real(dp) :: time_start, time_end, time_start0
 
      real(dp), external :: norm
@@ -800,19 +803,23 @@ subroutine berry_curvarture_cube
      real(dp), allocatable :: W(:)
      real(dp), allocatable :: eigval_allk(:, :), eigval_allk_mpi(:, :)
 
-     knv3=Nk1*Nk2*Nk3
+     if (BerryCurvature_Cube_calc) then
+        n_kpoints=Nk1*Nk2*Nk3
+     else
+        n_kpoints= nk3_band 
+     endif
 
      allocate(Vmn_wann(Num_wann, Num_wann, 3), Vmn_Ham(Num_wann, Num_wann, 3))
      allocate(Dmn_Ham(Num_wann,Num_wann,3), Vmn_Ham_nondiag(Num_wann, Num_wann, 3))
      allocate(W(Num_wann))
      allocate(UU(Num_wann, Num_wann))
-     allocate(eigval_allk(Num_wann, knv3))
-     allocate(eigval_allk_mpi(Num_wann, knv3))
+     allocate(eigval_allk(Num_wann, n_kpoints))
+     allocate(eigval_allk_mpi(Num_wann, n_kpoints))
      allocate( Omega_BerryCurv(Num_wann, 3), m_OrbMag(Num_wann, 3))
-     allocate( Omega_allk    (Num_wann, 3, knv3))
-     allocate( Omega_allk_mpi(Num_wann, 3, knv3))
-     allocate( m_OrbMag_allk    (Num_wann, 3, knv3))
-     allocate( m_OrbMag_allk_mpi(Num_wann, 3, knv3))
+     allocate( Omega_allk    (Num_wann, 3, n_kpoints))
+     allocate( Omega_allk_mpi(Num_wann, 3, n_kpoints))
+     allocate( m_OrbMag_allk    (Num_wann, 3, n_kpoints))
+     allocate( m_OrbMag_allk_mpi(Num_wann, 3, n_kpoints))
      Omega_BerryCurv= 0d0
      m_OrbMag=0d0
      
@@ -828,19 +835,27 @@ subroutine berry_curvarture_cube
      call now(time_start0)
      time_start= time_start0
      time_end  = time_start0
-     do ik= 1+ cpuid, knv3, num_cpu
+     do ik= 1+ cpuid, n_kpoints, num_cpu
         if (cpuid==0.and. mod(ik/num_cpu, 100)==0) &
            write(stdout, '(a, i9, "  /", i10, a, f10.1, "s", a, f10.1, "s")') &
-           ' Berry curvature: ik', ik, knv3, ' time left', &
-           (knv3-ik)*(time_end- time_start)/num_cpu, &
+           ' Berry curvature: ik', ik, n_kpoints, ' time left', &
+           (n_kpoints-ik)*(time_end- time_start)/num_cpu, &
            ' time elapsed: ', time_end-time_start0 
-        ikx= (ik-1)/(nk2*nk3)+1
-        iky= ((ik-1-(ikx-1)*Nk2*Nk3)/nk3)+1
-        ikz= (ik-(iky-1)*Nk3- (ikx-1)*Nk2*Nk3)
-        k= K3D_start_cube+ K3D_vec1_cube*(ikx-1)/dble(nk1)  &
-           + K3D_vec2_cube*(iky-1)/dble(nk2)  &
-           + K3D_vec3_cube*(ikz-1)/dble(nk3) 
-           !- (K3D_vec1_cube+ K3D_vec2_cube+ K3D_vec3_cube)/2d0
+
+        !> if we calculate BC in the BZ, we generate kpoints in the BZ
+        if (BerryCurvature_Cube_calc) then
+           !> kbulk mode
+           ikx= (ik-1)/(nk2*nk3)+1
+           iky= ((ik-1-(ikx-1)*Nk2*Nk3)/nk3)+1
+           ikz= (ik-(iky-1)*Nk3- (ikx-1)*Nk2*Nk3)
+           k= K3D_start_cube+ K3D_vec1_cube*(ikx-1)/dble(nk1)  &
+              + K3D_vec2_cube*(iky-1)/dble(nk2)  &
+              + K3D_vec3_cube*(ikz-1)/dble(nk3) 
+ 
+        elseif (BerryCurvature_kpath_sepband_calc) then
+           !> kpath mode
+           k= kpath_3d(:, ik)
+        endif
 
         call now(time_start)
         
@@ -857,7 +872,7 @@ subroutine berry_curvarture_cube
         call get_Dmn_Ham(W, Vmn_Ham, Dmn_Ham)
         call get_Vmn_Ham_nondiag(Vmn_Ham, Vmn_Ham_nondiag)
 
-        call berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+        call Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
         call orbital_magnetization_singlek_allbands(Dmn_Ham, Vmn_Ham_nondiag, m_OrbMag)
         Omega_allk(:, :, ik) = Omega_BerryCurv
         m_OrbMag_allk(:, :, ik) = m_OrbMag
@@ -878,6 +893,7 @@ subroutine berry_curvarture_cube
      m_OrbMag_allk_mpi= m_OrbMag_allk
 #endif
 
+     IF (BerryCurvature_Cube_calc) THEN
      !> write out Berry curvature and orbital magnetization to a file which
      !> can be open by software Fermisurfer.
      outfileindex= outfileindex+ 1
@@ -890,12 +906,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= Omega_allk_mpi(m, :, ik)/Angstrom2atomic**2
               write(outfileindex, '(E18.10)') norm(o1)
            enddo
@@ -916,12 +932,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= Omega_allk_mpi(m, :, ik)/Angstrom2atomic**2
               write(outfileindex, '(E18.10)') o1(1)
            enddo
@@ -942,12 +958,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= Omega_allk_mpi(m, :, ik)/Angstrom2atomic**2
               write(outfileindex, '(E18.10)') o1(3)
            enddo
@@ -966,12 +982,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= m_OrbMag_allk_mpi(m, :, ik)
               write(outfileindex, '(E18.10)') norm(o1)
            enddo
@@ -990,12 +1006,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= m_OrbMag_allk_mpi(m, :, ik)
               write(outfileindex, '(E18.10)') o1(3)
            enddo
@@ -1014,12 +1030,12 @@ subroutine berry_curvarture_cube
         write(outfileindex, '(3f12.6)') Origin_cell%Kub
         write(outfileindex, '(3f12.6)') Origin_cell%Kuc
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               write(outfileindex, '(E18.10)') eigval_allk_mpi(m, ik)-iso_energy
            enddo
         enddo
         do m=1, Num_wann
-           do ik= 1, knv3
+           do ik= 1, n_kpoints
               o1= m_OrbMag_allk_mpi(m, :, ik)
               write(outfileindex, '(E18.10)') o1(1)
            enddo
@@ -1027,6 +1043,41 @@ subroutine berry_curvarture_cube
         close(outfileindex)
      endif
 
+     ELSE
+
+
+     !> output the Berry curvature to file
+     outfileindex= outfileindex+ 1
+     if (cpuid==0) then
+        open(unit=outfileindex, file='Berrycurvature_line.dat')
+        write(outfileindex, '(a)')'# Column 1 kpath with accumulated length in the kpath, Coloum 2: energy'
+        write(outfileindex, '(a)')'# Column 2-4 Berry curvature (Ang^2)'
+        write(outfileindex, "('#column', i5, 20i16)")(i, i=1, 8)
+        write(outfileindex, '(20a16)')'# k (1/A)', " eig", &
+           'Omega_x', 'Omega_y', 'Omega_z', &
+           'm_x', 'm_y', 'm_z'
+        do i=1, Num_wann
+           do ik=1, n_kpoints
+              write(outfileindex, '(200f16.9)')k3len(ik)*Angstrom2atomic, eigval_allk_mpi(i, ik), &
+                 Omega_allk_mpi(i, :, ik), m_OrbMag_allk_mpi(i, :, ik)
+           enddo
+           write(outfileindex, *)' '
+        enddo
+        close(outfileindex)
+ 
+
+        close(outfileindex)
+     endif
+
+     !> minimum and maximum value of energy bands
+     emin=  minval(eigval_allk_mpi)-0.5d0
+     emax=  maxval(eigval_allk_mpi)+0.5d0
+
+     call generate_ek_kpath_gnu('Berrycurvature_line.dat', 'Berrycurvature_line.gnu', 'Berrycurvature_line.pdf', &
+                                 emin, emax, n_kpoints, Nk3lines, &
+                                 k3line_name, k3line_stop, k3len)
+   
+     ENDIF
 
 
 #if defined (MPI)
@@ -1035,9 +1086,9 @@ subroutine berry_curvarture_cube
 
      return
 
-  end subroutine berry_curvarture_cube
+  end subroutine berry_curvature_cube
 
-  subroutine berry_curvarture_plane_full
+  subroutine Berry_curvature_plane_full
      !> Calculate Berry curvature and orbital magnetization for the selected bands
      !> ref : Physical Review B 74, 195118(2006)
      !> eqn (34)
@@ -1126,7 +1177,7 @@ subroutine berry_curvarture_plane_full
         k= kslice(:, ik)
 
         call now(time_start)
-        call berry_curvarture_orb_mag_singlek_allbands_pack(k, Omega_BerryCurv, m_OrbMag, W)
+        call Berry_curvature_orb_mag_singlek_allbands_pack(k, Omega_BerryCurv, m_OrbMag, W)
        
         do i=1, 3
            Omega_allk_Occ(1, i, ik) = sum(Omega_BerryCurv(1:NumOccupied, i))
@@ -1417,9 +1468,9 @@ subroutine berry_curvarture_plane_full
 
      return
 
-  end subroutine berry_curvarture_plane_full
+  end subroutine Berry_curvature_plane_full
 
-  subroutine berry_curvarture_orb_mag_singlek_allbands_pack(k, Omega_BerryCurv, m_OrbMag, W)
+  subroutine Berry_curvature_orb_mag_singlek_allbands_pack(k, Omega_BerryCurv, m_OrbMag, W)
      !> Calculate Berry curvature and orbital magnetization for the selected bands
      !> ref : Physical Review B 74, 195118(2006)
      !> eqn (34)
@@ -1477,19 +1528,19 @@ subroutine berry_curvarture_orb_mag_singlek_allbands_pack(k, Omega_BerryCurv, m_
      call get_Dmn_Ham(W, Vmn_Ham, Dmn_Ham)
      call get_Vmn_Ham_nondiag(Vmn_Ham, Vmn_Ham_nondiag)
 
-     call berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+     call Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
      call orbital_magnetization_singlek_allbands(Dmn_Ham, Vmn_Ham_nondiag, m_OrbMag)
 
      deallocate(Vmn_wann, Vmn_Ham, Vmn_Ham_nondiag)
      deallocate(UU, Dmn_Ham)
      return
 
-  end subroutine berry_curvarture_orb_mag_singlek_allbands_pack
+  end subroutine Berry_curvature_orb_mag_singlek_allbands_pack
 
 
 
 
-  subroutine berry_curvarture_plane_EF
+  subroutine Berry_curvature_plane_EF
      !> Calculate Berry curvature and orbital magnetization for the selected bands
      !> ref : Physical Review B 74, 195118(2006)
      !> eqn (34)
@@ -1600,7 +1651,7 @@ subroutine berry_curvarture_plane_EF
         call get_Dmn_Ham(W, Vmn_Ham, Dmn_Ham)
         call get_Vmn_Ham_nondiag(Vmn_Ham, Vmn_Ham_nondiag)
 
-        call berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+        call Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
         call orbital_magnetization_singlek_allbands(Dmn_Ham, Vmn_Ham_nondiag, m_OrbMag)
         Omega_allk(:, :, ik) = Omega_BerryCurv
         m_OrbMag_allk(:, :, ik) = m_OrbMag
@@ -1698,11 +1749,11 @@ subroutine berry_curvarture_plane_EF
  
      return
 
-  end subroutine berry_curvarture_plane_EF
+  end subroutine Berry_curvature_plane_EF
 
 
 
-  subroutine berry_curvarture_plane_selectedbands
+  subroutine Berry_curvature_plane_selectedbands
      !> Calculate Berry curvature and orbital magnetization for the selected bands
      !> ref : Physical Review B 74, 195118(2006)
      !> eqn (34)
@@ -1813,7 +1864,7 @@ subroutine berry_curvarture_plane_selectedbands
         call get_Dmn_Ham(W, Vmn_Ham, Dmn_Ham)
         call get_Vmn_Ham_nondiag(Vmn_Ham, Vmn_Ham_nondiag)
 
-        call berry_curvarture_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
+        call Berry_curvature_singlek_allbands(Dmn_Ham, Omega_BerryCurv)
         call orbital_magnetization_singlek_allbands(Dmn_Ham, Vmn_Ham_nondiag, m_OrbMag)
         Omega_allk(:, :, ik) = Omega_BerryCurv
         m_OrbMag_allk(:, :, ik) = m_OrbMag
@@ -1911,11 +1962,11 @@ subroutine berry_curvarture_plane_selectedbands
  
      return
 
-  end subroutine berry_curvarture_plane_selectedbands
+  end subroutine Berry_curvature_plane_selectedbands
 
 
 
-  subroutine berry_curvarture_plane
+  subroutine Berry_curvature_plane
      !> Calculate Berry curvature 
      !
      !> ref : Physical Review B 74, 195118(2006)
@@ -1990,11 +2041,11 @@ subroutine berry_curvarture_plane
         Omega_y= 0d0
         Omega_z= 0d0
 
-        !call berry_curvarture_singlek_numoccupied_old(k, Omega_x, Omega_y, Omega_z)
+        !call Berry_curvature_singlek_numoccupied_old(k, Omega_x, Omega_y, Omega_z)
         if (Berrycurvature_EF_calc) then
-           call berry_curvarture_singlek_EF(k, iso_energy, Omega_x, Omega_y, Omega_z)
+           call Berry_curvature_singlek_EF(k, iso_energy, Omega_x, Omega_y, Omega_z)
         else
-           call berry_curvarture_singlek_numoccupied_total(k, Omega_x(1), Omega_y(1), Omega_z(1))
+           call Berry_curvature_singlek_numoccupied_total(k, Omega_x(1), Omega_y(1), Omega_z(1))
         endif
         Omega(1, ik) = sum(Omega_x)
         Omega(2, ik) = sum(Omega_y)
@@ -2141,7 +2192,7 @@ subroutine berry_curvarture_plane
  
      return
 
-  end subroutine berry_curvarture_plane
+  end subroutine Berry_curvature_plane
 
   subroutine Fourier_R_to_k(k, ham)
      !> Fourier transform the Hamiltonian from R space to k space
@@ -2271,8 +2322,8 @@ subroutine Chern_sphere_single(k0, r0, Chern)
          ct=cos(theta)
          sp=sin(phi)
          cp=cos(phi)
-        !call berry_curvarture_singlek_numoccupied_old(k_direct, Omega_x, Omega_y, Omega_z)
-         call berry_curvarture_singlek_numoccupied    (k_direct, Omega_x, Omega_y, Omega_z)
+        !call Berry_curvature_singlek_numoccupied_old(k_direct, Omega_x, Omega_y, Omega_z)
+         call Berry_curvature_singlek_numoccupied    (k_direct, Omega_x, Omega_y, Omega_z)
          O_x= real(sum(Omega_x(1:Numoccupied)))
          O_y= real(sum(Omega_y(1:Numoccupied)))
          O_z= real(sum(Omega_z(1:Numoccupied)))
@@ -2416,7 +2467,7 @@ subroutine Chern_halftorus_single(k0, Rbig, rsmall_a, rsmall_b, Chern)
         sp=sin(phi)
         cp=cos(phi)
         rt= Rbig+ rsmall_a* ct
-        call berry_curvarture_singlek_numoccupied(k_direct, Omega_x, Omega_y, Omega_z)
+        call Berry_curvature_singlek_numoccupied(k_direct, Omega_x, Omega_y, Omega_z)
        !O_x= real(Omega_x(Numoccupied))
        !O_y= real(Omega_y(Numoccupied))
        !O_z= real(Omega_z(Numoccupied))
diff --git a/src/ek_bulk.f90 b/src/ek_bulk.f90
index 8331b98..ec9381e 100644
--- a/src/ek_bulk.f90
+++ b/src/ek_bulk.f90
@@ -42,7 +42,7 @@ subroutine ek_bulk_line
 
    do ik= 1+cpuid, knv3, num_cpu
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
       
       ! calculation bulk hamiltonian
       Hamk_bulk= 0d0
@@ -232,7 +232,7 @@ subroutine ek_bulk_line_valley
 
    do ik= 1+cpuid, knv3, num_cpu
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
       
       ! calculation bulk hamiltonian
       Hamk_bulk= 0d0
@@ -1249,7 +1249,7 @@ subroutine sparse_ekbulk
 
    do ik=1+ cpuid, nk3_band, num_cpu
       if (cpuid==0) write(stdout, '(a, 2i10)') 'BulkBand_calc in sparse mode:', ik,nk3_band
-      k3 = K3points(:, ik)
+      k3 = kpath_3d(:, ik)
       call now(time1)
       !> use atomic gauge here
       call ham_bulk_coo_sparsehr(k3,acoo,icoo,jcoo)
@@ -1509,7 +1509,7 @@ subroutine sparse_ekbulk_valley
    k3= 0
    do ik=1+ cpuid, nk3_band, num_cpu
       if (cpuid==0) write(stdout, '(a, 2i10)') 'BulkBand_calc in sparse mode:', ik,nk3_band
-      k3 = K3points(:, ik)
+      k3 = kpath_3d(:, ik)
       call now(time1)
       call ham_bulk_coo_sparsehr(k3,acoo,icoo,jcoo)
       !> change the energy unit from Hatree to eV
@@ -2042,7 +2042,7 @@ subroutine ekbulk_elpa
 
       haml_block=0.0d0
       z=0.0d0
-      k3  = k3points(:, ik)
+      k3  = kpath_3d(:, ik)
       call ham_bulk_elpa( k3, haml_block,&
          np_rows,np_cols,&
          na_rows,na_cols,&
@@ -2699,7 +2699,7 @@ subroutine ek_bulk_mengyu
       !if (cpuid==0) write(stdout, *)'BulkBand, ik, knv3 ', ik, knv3
 
       !> in fractional coordinates
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
 
       kxy_cart(:)= k(1)* Origin_cell%Kua+ k(2)* Origin_cell%Kub+ k(3)* Origin_cell%Kuc
       call rotate_k3_to_kplane(kxy_cart(:), kxy_plane(:))
@@ -2885,7 +2885,7 @@ subroutine ek_bulk_spin
    do ik= 1+cpuid, knv3, num_cpu
       if (cpuid==0) write(stdout, *)'BulkBandSpin, ik, knv3 ', ik, knv3
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
 
       ! calculation bulk hamiltonian
       Hamk_bulk= 0d0
@@ -3009,14 +3009,14 @@ subroutine ek_bulk_mirror_z
    do ik= 1+cpuid, knv3, num_cpu
       if (cpuid==0) write(stdout, *)'BulkBandmirrorz, ik, knv3 ', ik, knv3
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
 
       ! calculation bulk hamiltonian
       Hamk_bulk= 0d0
       call ham_bulk_atomicgauge    (k, Hamk_bulk)
       Hamk_bulk= Hamk_bulk/eV2Hartree
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
       !k(2)= -k(2)
       Hamk= 0d0
       call ham_bulk_latticegauge    (k, Hamk)
@@ -3177,14 +3177,14 @@ subroutine ek_bulk_mirror_x
    do ik= 1+cpuid, knv3, num_cpu
       if (cpuid==0) write(stdout, *)'EkBulk_mirror, ik, knv3 ', ik, knv3
 
-      k = k3points(:, ik)
+      k = kpath_3d(:, ik)
 
       ! calculation bulk hamiltonian
       Hamk_bulk= 0d0
       call ham_bulk_latticegauge    (k, Hamk_bulk)
       Hamk_bulk= Hamk_bulk/eV2Hartree
 
-      !k = k3points(:, ik)
+      !k = kpath_3d(:, ik)
       !k(1)= -k(1)
       !Hamk= 0d0
       !call ham_bulk_latticegauge    (k, Hamk)
@@ -3373,6 +3373,7 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
    if (cpuid==0) then
       open(unit=outfileindex, file=gnufilename)
       write(outfileindex, '(a)') '# requirement: gnuplot version>5.4'
+      write(outfileindex, '(2a)') '# Please open the data file to check the data:   ', trim(adjustl(datafilename))
       write(outfileindex, '(a)') 'set terminal pdf enhanced color font ",24"'
       write(outfileindex,'(2a)') 'set palette defined ( 0  "green", ', &
          '5 "yellow", 10 "red" )'
@@ -3407,20 +3408,20 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
       enddo
       write(outfileindex, '(2a)')"# please comment the following lines to plot the fatband "
       if (Is_Sparse_Hr) then
-         write(outfileindex, '(2a)')"plot 'bulkek.dat' u 1:2 ",  &
+         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
-         write(outfileindex, '(2a)')"plot 'bulkek.dat' u 1:2 ",  &
+         write(outfileindex, '(4a)')"plot '", trim(adjustl(datafilename)), "' u 1:2 ",  &
             " w lp lw 2 pt 7  ps 0.2 lc rgb 'black', 0 w l lw 2 dt 2"
       endif
       write(outfileindex, '(2a)')" " 
       write(outfileindex, '(2a)')"# uncomment the following lines to plot the fatband "
-      write(outfileindex, '(2a)')"#plot 'bulkek.dat' u 1:2:3 ",  &
+      write(outfileindex, '(4a)')"#plot '", trim(adjustl(datafilename)), "' u 1:2:3 ",  &
          " w lp lw 2 pt 7  ps 0.2 lc palette, 0 w l lw 2 dt 2"
-      write(outfileindex, '(2a)')"# uncomment the following lines to plot the spin if necessary"
-      write(outfileindex, '(2a)')"#plot 'bulkek.dat' u 1:2 ",  &
+      write(outfileindex, '(4a)')"# uncomment the following lines to plot the spin if necessary"
+      write(outfileindex, '(4a)')"#plot '", trim(adjustl(datafilename)), "' u 1:2 ",  &
          "w lp lw 2 pt 7  ps 0.2, \"
-      write(outfileindex, '(2a)')"     'bulkek.dat' u 1:2:($3/6):($4/6) ",  &
+      write(outfileindex, '(4a)')"#plot '", trim(adjustl(datafilename)), "' u 1:2:($3/6):($4/6) ",  &
          "w vec"
       close(outfileindex)
    endif
@@ -3430,4 +3431,4 @@ subroutine generate_ek_kpath_gnu(datafilename, gnufilename, gnuoutfilename, &
 204 format('set arrow from ',F10.5,',',A5,' to ',F10.5,',',A5, ' nohead lw 2')
 
     return
-end subroutine
+end subroutine generate_ek_kpath_gnu
diff --git a/src/lanczos_sparse.f90 b/src/lanczos_sparse.f90
index 5c7fff6..bab8f20 100644
--- a/src/lanczos_sparse.f90
+++ b/src/lanczos_sparse.f90
@@ -700,7 +700,7 @@ subroutine LandauLevel_k_dos_Lanczos
    use sparse
    use wmpi
    use para, only : Magq, Num_Wann, Bx, By, zi, pi, Fermi_broadening, Angstrom2atomic, &
-      OmegaNum, OmegaMin, OmegaMax, nk3_band, Magp, stdout, k3points, eV2Hartree, &
+      OmegaNum, OmegaMin, OmegaMax, nk3_band, Magp, stdout, kpath_3d, eV2Hartree, &
       outfileindex, K3len_mag,splen,Is_Sparse_Hr,ijmax,NumLCZVecs, MagneticSuperProjectedArea, &
       Nk3lines, k3line_mag_stop, k3line_name, NumRandomConfs
    implicit none
@@ -823,7 +823,7 @@ subroutine LandauLevel_k_dos_Lanczos
          InitialVector= InitialVector/dsqrt(dble(norm))
 
          if (cpuid==0) write(stdout, '(a, 2i10)') 'LandauLevel_k_DOS', ik, nk3_band
-         k3= k3points(:, ik)
+         k3= kpath_3d(:, ik)
          nnz= nnzmax
          if(Is_Sparse_Hr) then
             call ham_3Dlandau_sparseHR(nnz,Mdim,NQ,k3,acsr,jcsr,icsr)
@@ -930,7 +930,7 @@ subroutine bulkbandk_dos_lanczos
    use sparse
    use wmpi
    use para, only : Magq, Num_Wann, Bx, By, zi, pi, Fermi_broadening, Angstrom2atomic, &
-      OmegaNum, OmegaMin, OmegaMax, nk3_band, Magp, stdout, k3points, &
+      OmegaNum, OmegaMin, OmegaMax, nk3_band, Magp, stdout, kpath_3d, &
       outfileindex, K3len,splen,Is_Sparse_Hr,ijmax,NumLCZVecs, eV2Hartree
    implicit none
 
@@ -1016,7 +1016,7 @@ subroutine bulkbandk_dos_lanczos
       InitialVector= InitialVector/dsqrt(dble(norm))
 
       if (cpuid==0) write(stdout, '(a, 2i10)') 'LandauLevel_k_DOS', ik, nk3_band
-      k3= k3points(:, ik)
+      k3= kpath_3d(:, ik)
       nnz= nnzmax
       if(Is_Sparse_Hr) then
          call ham_bulk_coo_sparsehr(k3, acsr,jcsr,icsr)
diff --git a/src/landau_level.f90 b/src/landau_level.f90
index 559069f..48c0124 100644
--- a/src/landau_level.f90
+++ b/src/landau_level.f90
@@ -140,7 +140,7 @@ subroutine landau_level_k
 
       call now(time_start)
 
-      k3= k3points(:, ik)
+      k3= kpath_3d(:, ik)
       call ham_3Dlandau(Ndimq, Nq, k3, ham_landau)
       ham_landau= ham_landau/eV2Hartree
 
diff --git a/src/landau_level_sparse.f90 b/src/landau_level_sparse.f90
index 18e54d4..a7cd37f 100644
--- a/src/landau_level_sparse.f90
+++ b/src/landau_level_sparse.f90
@@ -1054,7 +1054,7 @@ subroutine sparse_landau_level_k
       call now(time_start)
 
       if (cpuid==0) write(stdout, '(a, 2i10)') 'LandauLevel_k_calc', ik,nk3_band
-      k3 = K3points(:, ik)
+      k3 = kpath_3d(:, ik)
       nnz= nnzmax
       call now(time1)
       acoo= 0d0
diff --git a/src/main.f90 b/src/main.f90
index 2469eb9..6cf47d8 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -482,7 +482,7 @@ program main
         if(cpuid.eq.0)write(stdout, *)' '
         if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry curvature for a slab system'
         call now(time_start)
-        call berry_curvarture_slab
+        call Berry_curvature_slab
         call now(time_end)
         call print_time_cost(time_start, time_end, 'BerryCurvature_slab')
         if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature for a slab system'
@@ -492,7 +492,7 @@ program main
         if(cpuid.eq.0)write(stdout, *)' '
         if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry curvature'
         call now(time_start)
-        call berry_curvarture_line
+        call Berry_curvature_line_occupied
         call now(time_end)
         call print_time_cost(time_start, time_end, 'BerryCurvature')
         if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature'
@@ -502,17 +502,17 @@ program main
         if(cpuid.eq.0)write(stdout, *)' '
         if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry curvature'
         call now(time_start)
-        call berry_curvarture_plane_full
+        call Berry_curvature_plane_full
         call now(time_end)
         call print_time_cost(time_start, time_end, 'BerryCurvature')
         if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature'
      endif
 
-     if (BerryCurvature_Cube_calc)then
+     if (BerryCurvature_Cube_calc.or.BerryCurvature_kpath_sepband_calc)then
         if(cpuid.eq.0)write(stdout, *)' '
         if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry curvature in a k-cube'
         call now(time_start)
-        call berry_curvarture_cube
+        call Berry_curvature_cube
         call now(time_end)
         call print_time_cost(time_start, time_end, 'BerryCurvature_Cube')
         if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature in a cube'
diff --git a/src/module.f90 b/src/module.f90
index ed0e58b..2f55d75 100644
--- a/src/module.f90
+++ b/src/module.f90
@@ -415,6 +415,7 @@ module para
      logical :: BerryCurvature_plane_selectedbands_calc ! Flag for Berry curvature calculation
      logical :: BerryCurvature_EF_calc ! Flag for Berry curvature calculation
      logical :: BerryCurvature_kpath_EF_calc ! Flag for Berry curvature calculation in kpath model at EF
+     logical :: BerryCurvature_kpath_sepband_calc ! Flag for Berry curvature calculation in kpath model for each band
      logical :: BerryCurvature_kpath_Occupied_calc ! Flag for Berry curvature calculation in kpath model sum over all occupied bands
      logical :: BerryCurvature_Cube_calc ! Flag for Berry curvature calculation
      logical :: BerryCurvature_slab_calc ! Flag for Berry curvature calculation for a slab system
@@ -469,7 +470,7 @@ module para
                           ChargeDensity_selected_bands_calc, &
                           ChargeDensity_selected_energies_calc, &
                           WireBand_calc, Wilsonloop_calc, &
-                          WannierCenter_calc,BerryPhase_calc, &
+                          WannierCenter_calc,BerryPhase_calc, BerryCurvature_kpath_sepband_calc, &
                           BerryCurvature_EF_calc, BerryCurvature_calc, &
                           Berrycurvature_kpath_EF_calc, BerryCurvature_kpath_Occupied_calc, &
                           BerryCurvature_plane_selectedbands_calc, &
@@ -818,14 +819,14 @@ module para
 
      ! k list for 3D case band
      integer :: nk3lines  ! Howmany k lines for bulk band calculation
-     integer :: nk3_band  ! Howmany k points for each k line
+     integer :: nk3_band  ! Howmany k points for whole kpath
      character(4), allocatable :: k3line_name(:) ! Name of the K points
      real(dp),allocatable :: k3line_stop(:)  ! Connet points
      real(dp),allocatable :: k3line_start(:, :) ! Start point for each k line
      real(dp),allocatable :: k3line_end(:, :) ! End point for each k line
      real(dp),allocatable :: K3list_band(:, :) ! coordinate of k points for bulk band calculation in kpath mode
      real(dp),allocatable :: K3len(:)  ! put all k points in a line in order to plot the bands 
-     real(dp),allocatable :: K3points(:, :) ! coordinate of k points for bulk band calculation in cube mode
+     real(dp),allocatable :: kpath_3d(:, :) ! coordinate of k points for bulk band calculation in cube mode
 
      !> k points in the point mode 
      integer :: Nk3_point_mode
diff --git a/src/readinput.f90 b/src/readinput.f90
index c82f2a0..15c9c78 100644
--- a/src/readinput.f90
+++ b/src/readinput.f90
@@ -125,6 +125,7 @@ subroutine readinput
    BerryCurvature_Cube_calc   = .FALSE.
    BerryCurvature_slab_calc = .FALSE.
    Berrycurvature_kpath_EF_calc = .FALSE.
+   BerryCurvature_kpath_sepband_calc = .FALSE.
    BerryCurvature_kpath_Occupied_calc = .FALSE.
    MirrorChern_calc      = .FALSE.
    Dos_calc              = .FALSE.
@@ -184,6 +185,7 @@ subroutine readinput
       write(*, *)"SlabSpintexture,wanniercenter_calc, Wilsonloop_calc, "
       write(*, *)"BerryPhase_calc,BerryCurvature_calc, BerryCurvature_EF_calc"
       write(*, *)"Berrycurvature_kpath_EF_calc, BerryCurvature_kpath_Occupied_calc"
+      write(*, *)"BerryCurvature_kpath_sepband_calc"
       write(*, *)"BerryCurvature_slab_calc, BerryCurvature_Cube_calc"
       write(*, *)"Dos_calc, JDos_calc, FindNodes_calc"
       write(*, *)"BulkFS_plane_calc"
@@ -2015,11 +2017,11 @@ subroutine readinput
    allocate(k3len(nk3_band))
    allocate(k3len_mag(nk3_band))
    allocate(k3len_unfold(nk3_band))
-   allocate(k3points(3, nk3_band))
+   allocate(kpath_3d(3, nk3_band))
    k3len=0d0
    k3len_mag=0d0
    k3len_unfold=0d0
-   k3points= 0d0
+   kpath_3d= 0d0
    t1= 0d0
    do j=1, nk3lines
       do i=1, NN
@@ -2030,7 +2032,7 @@ subroutine readinput
          !k1= kstart
          !k2= kend
 
-         k3points(:, i+ (j-1)*NN)= kstart+ (kend- kstart)*dble(i-1)/dble(NN-1)
+         kpath_3d(:, i+ (j-1)*NN)= kstart+ (kend- kstart)*dble(i-1)/dble(NN-1)
 
          temp= dsqrt((k2(1)- k1(1))**2 &
             +(k2(2)- k1(2))**2  &
diff --git a/src/sigma_OHE.f90 b/src/sigma_OHE.f90
index b47465d..f6588f8 100644
--- a/src/sigma_OHE.f90
+++ b/src/sigma_OHE.f90
@@ -1831,7 +1831,7 @@ subroutine evolve_k_ohe
         NSlice_Btau_all_mpi = 0; NSlice_Btau_all=0
         kout_all= 0d0; kout_all_mpi= 0d0
         do ik= 1+cpuid, nk3_band, num_cpu
-           k_start = k3points(:, ik)
+           k_start = kpath_3d(:, ik)
            if (cpuid.eq.0) then
               write( stdout, *)'ib, ik', ib, ik
            endif
@@ -1871,7 +1871,7 @@ subroutine evolve_k_ohe
  
 
         do ik= 1, nk3_band
-           k_start = k3points(:, ik)
+           k_start = kpath_3d(:, ik)
            call velocity_calc(Nband_Fermi_Level, bands_fermi_level, k_start, velocity_k, Ek)
            if (cpuid.eq.0) then
            endif
diff --git a/src/unfolding.f90 b/src/unfolding.f90
index a869922..c41293d 100644
--- a/src/unfolding.f90
+++ b/src/unfolding.f90
@@ -136,7 +136,7 @@ subroutine unfolding_kpath
    !> first unfold the kpoints from the kpath of the supercell
    do ik= 1+cpuid, nk3_band, num_cpu
       if (cpuid==0) write(stdout, '(a, i10," /", i10)') 'BulkBand unfolding at :', ik, nk3_band
-      k_PBZ_direct= k3points(:, ik)
+      k_PBZ_direct= kpath_3d(:, ik)
       call direct_cart_rec_unfold(k_PBZ_direct, k_cart)
       if (Landaulevel_unfold_line_calc) then 
          call cart_direct_rec_magneticcell(k_cart, k_PBZ_direct_in_SBZ)