diff --git a/examples/3DWeyl-model/wt.in-BerryCurvature_kpath_sepband_calc b/examples/3DWeyl-model/wt.in-BerryCurvature_kpath_sepband_calc
new file mode 100644
index 00000000..3153dd24
--- /dev/null
+++ b/examples/3DWeyl-model/wt.in-BerryCurvature_kpath_sepband_calc
@@ -0,0 +1,50 @@
+&TB_FILE
+Hrfile = "Weyl3D_hr.dat"
+/
+
+!> bulk band structure calculation flag
+&CONTROL
+BerryCurvature_kpath_sepband_calc
+/
+
+&SYSTEM
+NumOccupied = 1         ! NumOccupied
+SOC = 1                 ! soc
+E_FERMI = 0        ! e-fermi
+/
+
+&PARAMETERS
+Nk1 = 60            ! number k points 
+/
+
+LATTICE
+Angstrom
+   1.0000000   000000000   000000000    
+   000000000   1.0000000   000000000    
+   000000000   000000000   1.0000000    
+
+ATOM_POSITIONS
+1                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+ A 0     0     0. 
+
+PROJECTORS
+ 1           ! number of projectors
+A s
+
+
+SURFACE            ! See doc for details
+ 0  0  1
+ 1  0  0
+ 0  1  0
+
+KPLANE_BULK
+Direct
+ 0.00  0.00  0.00   ! Center for 3D k slice
+ 1.00  0.00  0.00   ! The first vector is along x direction
+ 0.00  0.00  1.00   ! The second vector is along z direction
+
+KPATH_BULK            ! k point path
+2              ! number of k line only for bulk band
+X 0.50000 0.00000 0.0000 G 0.00000 0.00000 0.0000
+G 0.00000 0.00000 0.0000 Z 0.00000 0.00000 0.5000  
diff --git a/examples/Fe/wt.in-BerryCurvature_kpath_sepband_calc b/examples/Fe/wt.in-BerryCurvature_kpath_sepband_calc
new file mode 100755
index 00000000..d4671617
--- /dev/null
+++ b/examples/Fe/wt.in-BerryCurvature_kpath_sepband_calc
@@ -0,0 +1,48 @@
+&TB_FILE
+Hrfile = wannier90_hr.dat
+/
+
+
+&CONTROL
+BerryCurvature_kpath_sepband_calc = T
+/
+
+&SYSTEM
+SOC = 1                 ! soc
+E_FERMI = 9.1839        ! e-fermi
+/
+
+&PARAMETERS
+Fermi_broadening = 0.01     ! infinite small value, like brodening 
+Nk1 = 61            ! number k points 
+/
+
+LATTICE
+Angstrom
+   1.4350000     1.4350000    -1.4350000
+  -1.4350000     1.4350000     1.4350000
+   1.4350000    -1.4350000     1.4350000
+
+ATOM_POSITIONS
+1                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+ Fe      0.0000000     0.0000000     0.000
+
+PROJECTORS
+ 8          ! number of projectors
+Fe pz px py dz2 dxz dyz dx2-y2 dxy
+
+SURFACE            ! See doc for details
+ 1  0  0
+ 0  1  0
+ 0  0  1
+
+KPATH_BULK            ! k point path
+6              ! number of k line only for bulk band
+G 0.0 0.0 0.0   H 0.5 -0.5 0.5   
+H 0.5 -0.5 0.5   N 0.0 0.0 0.5   
+N 0.0 0.0 0.5   G 0.0 0.0 0.0   
+G 0.0 0.0 0.0   P 0.25 0.25 0.25   
+P 0.25 0.25 0.25   H 0.5 -0.5 0.5   
+H 0.25 0.25 0.25   N 0.0 0.0 0.5  
+
diff --git a/examples/input.dat-2d b/examples/input.dat-2d
deleted file mode 100644
index e5385456..00000000
--- a/examples/input.dat-2d
+++ /dev/null
@@ -1,78 +0,0 @@
-&TB_FILE
-!Hrfile = YbMnBi2_TB_nsoc_hr.dat 
-Hrfile = new_hr.dat 
-/
-
-
-&CONTROL
-!> bulk band structure calculation flag
-BulkBand_calc         = T
-RibbonBand_calc       = F
-SlabSS_calc           = F
-Dos_calc              = T
-wanniercenter_calc    = F
-/
-
-&SYSTEM
-NSLAB = 10 
-NumOccupied = 8         ! NumOccupied
-SOC = 1                 ! soc
-E_FERMI = -0.2           ! e-fermi
-surf_onsite= 0.0        ! surf_onsite
-/
-
-&PARAMETERS
-Nk1 = 601           ! number k points 
-Nk2 = 601           ! number k points 
-NP = 1              ! number of principle layers
-OmegaMin=-2
-OmegaMax= 2
-OmegaNum= 4001
-/
-
-LATTICE
-Angstrom
- 4.460   0.000000      ! crystal lattice information
- 0.000   4.460000 
-
-ATOM_POSITIONS
-4                               ! number of atoms for projectors
-Direct                          ! Direct or Cartisen coordinate
- Bi 0.000000000         0.000000000         0.000000000  ! atom position
- Bi 0.500000000         0.500000000         0.000000000  ! atom position
- Yb 0.500000000         0.000000000         0.243000000  ! atom position
- Yb 0.000000000         0.500000000        -0.243000000  ! atom position
-
-PROJECTORS
-3 3 5 5     ! number of projectors for each atom
-Bi px py pz
-Bi px py pz
-Yb dxy dyz dxz dx2-y2 dz2
-Yb dxy dyz dxz dx2-y2 dz2
-
-EDGE            ! See doc for details
- 1  0 
- 0  1 
-
-KPATH_BULK            ! k point path
-3        ! numker of k line for 2D case
-G 0.0 0.0 M 0.5 0.5
-M 0.5 0.5 X 0.5 0.0  ! k path for 2D case
-X 0.5 0.0 G 0.0 0.0  ! k path for 2D case
-
-3        ! numker of k line for 2D case
-G 0.0 0.0 X 0.5 0.0
-X 0.5 0.0 M 0.5 0.5  ! k path for 2D case
-M 0.5 0.5 G 0.0 0.0  ! k path for 2D case
-
-KPATH_RIBBON            ! k point path
-2        ! numker of k line for 2D case
-X -0.5 G 0.0  ! k path for 2D case
-G  0.0 X 0.5
-
-
-KPLANE_BULK
--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
-
diff --git a/src/berrycurvature.f90 b/src/berrycurvature.f90
index 4ea46a48..1e2ced29 100644
--- a/src/berrycurvature.f90
+++ b/src/berrycurvature.f90
@@ -1049,7 +1049,7 @@ subroutine berry_curvature_cube
      !> output the Berry curvature to file
      outfileindex= outfileindex+ 1
      if (cpuid==0) then
-        open(unit=outfileindex, file='Berrycurvature_line.dat')
+        open(unit=outfileindex, file='Berrycurvature_line_sepband.txt')
         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)
@@ -1058,7 +1058,7 @@ subroutine berry_curvature_cube
            '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), &
+              write(outfileindex, '(200E16.6)')k3len(ik)*Angstrom2atomic, eigval_allk_mpi(i, ik)/eV2Hartree, &
                  Omega_allk_mpi(i, :, ik), m_OrbMag_allk_mpi(i, :, ik)
            enddo
            write(outfileindex, *)' '
@@ -1070,13 +1070,60 @@ subroutine berry_curvature_cube
      endif
 
      !> minimum and maximum value of energy bands
-     emin=  minval(eigval_allk_mpi)-0.5d0
-     emax=  maxval(eigval_allk_mpi)+0.5d0
+     emin=  minval(eigval_allk_mpi)/eV2Hartree-0.5d0
+     emax=  maxval(eigval_allk_mpi)/eV2Hartree+0.5d0
+
+      !> write script for gnuplot
+     outfileindex= outfileindex+ 1
+     if (cpuid==0) then
+        open(unit=outfileindex, file='Berrycurvature_line_sepband.gnu')
+        write(outfileindex, '(a)') '# requirement: gnuplot version>5.4'
+        write(outfileindex, '(2a)') '# Please open the data file to check the data: Berrycurvature_line_sepband.txt  '
+        write(outfileindex, '(a)') 'set terminal pdf enhanced color font ",24"'
+        write(outfileindex,'(2a)') 'set palette defined ( 0  "green", ', &
+           '5 "yellow", 10 "red" )'
+        write(outfileindex, '(3a)')"set output 'Berrycurvature_line_sepband.pdf'"
+        write(outfileindex, '(a)')'set style data linespoints'
+        write(outfileindex, '(a)')'unset key'
+        write(outfileindex, '(a)')'set pointsize 0.8'
+        write(outfileindex, '(a)')'#set xtics font ",24"'
+        write(outfileindex, '(a)')'#set ytics font ",24"'
+        write(outfileindex, '(a)')'#set ylabel font ",24"'
+        write(outfileindex, '(a)')'set ylabel offset 0.5,0'
+        write(outfileindex, '(a)')'set border lw 2'
+        write(outfileindex, '(a, f10.5, a)')'set xrange [0: ', maxval(k3len*Angstrom2atomic), ']'
+        write(outfileindex, '(a,f12.6)')'emin=', emin
+        write(outfileindex, '(a,f12.6)')'emax=', emax
+        write(outfileindex, '(a)')' '
+        write(outfileindex, '(a)')'#set cbrange if the number is too large meeting some band crossings'
+        write(outfileindex, '(a)')'set cbrange [-100:100]'
+        write(outfileindex, '(a)')' '
+        if (index(Particle,'phonon')/=0) then
+           write(outfileindex, '(a)')'set yrange [0: emax]'
+           write(outfileindex, '(a)')'set ylabel "Frequency (THz)"'
+        else
+           write(outfileindex, '(a)')'set ylabel "Energy (eV)"'
+           write(outfileindex, '(a)')'set yrange [ emin : emax ]'
+        endif
+        write(outfileindex, 202, advance="no") (k3line_name(i), k3line_stop(i)*Angstrom2atomic, i=1, Nk3lines)
+        write(outfileindex, 203)k3line_name(Nk3lines+1), k3line_stop(Nk3lines+1)*Angstrom2atomic
+  
+        do i=1, Nk3lines-1
+           if (index(Particle,'phonon')/=0) then
+              write(outfileindex, 204)k3line_stop(i+1)*Angstrom2atomic, '0.0', k3line_stop(i+1)*Angstrom2atomic, 'emax'
+           else
+              write(outfileindex, 204)k3line_stop(i+1)*Angstrom2atomic, 'emin', k3line_stop(i+1)*Angstrom2atomic, 'emax'
+           endif
+        enddo
+        write(outfileindex, '(4a)')"plot 'Berrycurvature_line_sepband.txt' u 1:2:5 ",  &
+           " w lp lw 2 pt 7  ps 0.2 lc palette, 0 w l lw 2 dt 2"
+        close(outfileindex)
+     endif
+
+202 format('set xtics (',20('"',A3,'" ',F10.5,','))
+203 format(A3,'" ',F10.5,')')
+204 format('set arrow from ',F10.5,',',A5,' to ',F10.5,',',A5, ' nohead lw 2')
 
-     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
 
 
diff --git a/src/main.f90 b/src/main.f90
index 6cf47d85..e54c32ee 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -508,7 +508,7 @@ program main
         if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature'
      endif
 
-     if (BerryCurvature_Cube_calc.or.BerryCurvature_kpath_sepband_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)