diff --git a/doc/source/features.rst b/doc/source/features.rst
index 0ab28f86..966ff77e 100644
--- a/doc/source/features.rst
+++ b/doc/source/features.rst
@@ -82,7 +82,7 @@ Consequently, there are 30 weights (for Bi2Se3), one for each projectors.
 Line mode
 ------------
 
-Calculate bulk energy band for a series k lines. This is the basic calculation after the
+Calculate bulk energy band for a series k lines and project the selected orbitals onto the band structures. This is the basic calculation after the
 construction of Wannier functions. You have to compare your Wannier interpolated bands 
 with the DFT bands. Those two bands should match well around the Fermi level.
 This comparison can help verify the accuracy of the Wannier functions, 
@@ -90,6 +90,11 @@ and ensure that they provide an appropriate description of the electronic struct
 In order to obtain the parameter **NumOccupied**, you also need to calculate the bulk energy bands and plot it with 
 the software xmgrace.
 
+In order to get the weight onto the projections, we need to set **SELECTED_WANNIERORBITALS** card. it
+The orbital index is corresponding to the **PROJECTOS** card, nothing related to spin freedom. For example
+In Bi2Se3, there are 30 Wannier functions in the hr.dat, however, there are only 15 projectors (orbitals). 
+So indices in **SELECTED_WANNIERORBITALS** are from 1 to 15.
+
 .. _bulkekin:
 
 Input
@@ -103,6 +108,12 @@ Typical flags for bulk band calculation in the wt.in. ::
   Nk1 = 101   ! Number of k points for each k line
   /
 
+  ! get projected bands onto different orbitals, here we only consider orbital and omit the spin freedom
+  SELECTED_WANNIERORBITALS
+  2
+  1-6                ! Bi
+  7-15               ! Se
+
   KPATH_BULK     ! k point path
   4              ! number of k lines only for bulk band
   G 0.00000 0.00000 0.0000 Z 0.00000 0.00000 0.5000
@@ -141,31 +152,74 @@ The data structure for **bulkek.dat** ::
 
 The subrotine for this feature is ek_bulk.f90 .
 
+Examples
+>>>>>>>>>
+
+Two examples Bi2Se3 and Graphene were prepared:
+
+Strong topological insulator Bi2Se3  (in order to get the fat-band plot, you need to modify bulkek.gnu following the instruction inside)
+
+.. code:: console
+   
+   $ cd examples/Bi2Se3
+   $ cp wt.in-bands wt.in
+   $ ../../bin/wt.x
+   $ gnuplot bulkek.gnu
+
+2D Dirac semimetal Graphene
+
+.. code:: console
+   
+   $ cd examples/Graphene
+   $ cp wt.in-bands wt.in
+   $ ../../bin/wt.x
+   $ gnuplot bulkek.gnu
+
+The band plots of Bi2Se3 and Graphene should be like this
+
+.. image:: images/bulkbands.png
+   :scale: 30 %
+
+.. note::
+
+   In order to get the latex symbol :math:`\Gamma` in gnuplot, you need to replace **"G"** to **"{/Symbol G}"** in the bulkek.gnu like this
+
+   set xtics ("M  "    0.00000,"K  "    0.84848, **"{/Symbol G}"**    2.54543,"M  "    4.01504)
 
 Plane mode
 --------------
 
 Calculate band structure in a k slice(plane) specified by KPLANE_BULK card.
-The mode is very useful to visualize the Dirac/Weyl cone. You have to set the following tags in **wt.in** ::
+The mode is very useful to visualize the Dirac/Weyl cone. You have to set the following necessary tags in **wt.in** ::
 
+
+  !> bulk band structure calculation flag
   &CONTROL
   BulkBand_plane_calc = T
   /
+  
+  &SYSTEM
+  NumOccupied = 1         ! only write out 4 bands around NumOccpuied
+  SOC = 0                 ! soc
+  E_FERMI =  -1.2533      ! e-fermi
+  /
+  
   &PARAMETERS
-  Nk1 = 51   ! Number of k points along the first vector in KPLANE_BULK
-  Nk2 = 51   ! Number of k points along the second vector in KPLANE_BULK
+  Nk1 = 201          ! number k points
+  Nk2 = 201          ! number k points
   /
-
+  
   KPLANE_BULK   ! fractional coordinates
-   0.00  0.00  0.30   ! Middle point for a k slice(plane) in 3D BZ. Usually, the position of Dirac points.
-   0.50  0.00  0.00   ! The first vector to define k plane(slice) in 3D BZ
-   0.00  0.50  0.00   ! The second vector to define k plane(slice) in 3D BZ
-
+     0.333333  0.333333  0.000000   ! Middle point for a k slice(plane) in 3D BZ. Usually, the position of Dirac points.
+     0.100000  0.000000  0.000000   ! The first vector to define k plane(slice) in 3D BZ
+     0.000000  0.100000  0.000000   ! The second vector to define k plane(slice) in 3D BZ
 
 The output file is **bulkek_plane.dat**, **bulkek_plane-matlab.dat** and **bulkek_plane.gnu**. You can get 
-**bulkek_plane.png** with ::
+**bulkek_plane.png** with 
+
+.. code:: console
 
-   gnuplot bulkek_plane.gnu
+   $ gnuplot bulkek_plane.gnu
 
 The **bulkek_plane-matlab.dat** is in MATLAB data format. You can plot the Dirac cone with matlab.  
 
@@ -181,8 +235,23 @@ where the x and y direction are line in the k plane and the z direction is perpe
 Column 7-10th are energies at each k point. Here we only print out 4 energy bands around the fermilevel. It depends on **NumOccupied**.
 Usually, I choose column 4th and 5th as k coordinates and choose 8 and 9 as energy bands to show the Dirac cone shown below.
 
-.. image:: images/bulk_plane.jpeg
-   :scale: 30 %
+
+Example
+--------------
+
+There is one example in the examples/Graphene
+
+.. code:: console
+
+   $ cd examples/Graphene
+   $ cp wt.in-bands-plane wt.in
+   $ ../../bin/wt.x
+   $ gnuplot bulkek_plane.gnu
+
+The output is bulkek_plane.png and should like this
+
+.. image:: images/graphene_dirac_cone.png
+   :scale: 50 %
 
 
 .. _bulkfscalculation:
@@ -192,9 +261,9 @@ Usually, I choose column 4th and 5th as k coordinates and choose 8 and 9 as ener
 Fermi surface calculation of the bulk system.
 
 When calculating systems with SOC=0, meaning no spin-orbit coupling, 
-WannierTools computes 12 bands from the *NumOccupied-5* to NumOccupied+6 energy levels. 
+WannierTools computes 12 bands from the **NumOccupied-5** to **NumOccupied+6** energy levels. 
 However, when calculating systems with SOC=1, meaning spin-orbit coupling is present, 
-WannierTools computes 16 bands from the *NumOccupied-7* to NumOccupied+8 energy levels.
+WannierTools computes 16 bands from the **NumOccupied-7** to **NumOccupied+8** energy levels.
 
 .. NOTE::
    1. In order to reduce the storage of the Fermi surface, we only write out few energy bands around NumOccupied'th bands.
diff --git a/doc/source/images/bulkbands.png b/doc/source/images/bulkbands.png
new file mode 100644
index 00000000..ca1ae3ee
Binary files /dev/null and b/doc/source/images/bulkbands.png differ
diff --git a/doc/source/images/graphene_dirac_cone.png b/doc/source/images/graphene_dirac_cone.png
new file mode 100644
index 00000000..cf57eebf
Binary files /dev/null and b/doc/source/images/graphene_dirac_cone.png differ
diff --git a/examples/Bi2Se3/wt.in b/examples/Bi2Se3/wt.in
index 8c3fdc51..ed4ade9e 100644
--- a/examples/Bi2Se3/wt.in
+++ b/examples/Bi2Se3/wt.in
@@ -1,6 +1,5 @@
 &TB_FILE
 Hrfile = 'wannier90_hr.dat'      
-!Hrfile = 'wannier90_hr.dat_nsymm12'
 Package = 'VASP'             ! obtained from VASP, it could be 'VASP', 'QE', 'Wien2k', 'OpenMx'
 /
 
@@ -27,45 +26,27 @@ Se pz px py
 Se pz px py
 Se pz px py
 
-SURFACE            ! See doc for details
- 1  0  0
- 0  1  0
-
-
 !> bulk band structure calculation flag
 &CONTROL
-SlabBandWaveFunc_calc = F
 BulkBand_calc         = T
-BulkBand_points_calc  = F
-SlabBand_calc         = T
-WireBand_calc         = F
-SlabSS_calc           = F
-SlabArc_calc          = F
-Z2_3D_calc            = F
-SlabSpintexture_calc  = F
-Wanniercenter_calc    = F
-SHC_calc              = F
+BulkBand_points_calc  = T
 /
 
 &SYSTEM
-NSLAB = 20              ! for thin film system
-NSLAB1= 4               ! nanowire system 
-NSLAB2= 4               ! nanowire system 
 NumOccupied = 18        ! NumOccupied
 SOC = 1                 ! soc
 E_FERMI = 4.4195        ! e-fermi
 /
 
+! get projected bands onto different orbitals, here we only consider orbital and omit the spin freedom
+SELECTED_WANNIERORBITALS
+2 
+1-6                ! Bi
+7-15               ! Se
+
+
 &PARAMETERS
-Eta_Arc = 0.050     ! infinite small value, like brodening 
-E_arc = 0.0         ! energy for calculate Fermi Arc
-OmegaNum = 300      ! omega number       
-OmegaMin = -5.0     ! energy interval
-OmegaMax =  5.0     ! energy interval
 Nk1 = 41            ! number k points  odd number would be better
-Nk2 = 41            ! number k points  odd number would be better
-Nk3 = 41            ! number k points  odd number would be better
-NP = 2              ! number of principle layers
 /
 
 KPATH_BULK            ! k point path
@@ -75,34 +56,6 @@ Z 0.00000 0.00000 0.5000 F 0.50000 0.50000 0.0000
 F 0.50000 0.50000 0.0000 G 0.00000 0.00000 0.0000
 G 0.00000 0.00000 0.0000 L 0.50000 0.00000 0.0000  
 
-KPATH_SLAB
-2        ! numker of k line for 2D case
-K 0.33 0.67 G 0.0 0.0  ! k path for 2D case
-G 0.0 0.0 M 0.5 0.5
-
-KPLANE_SLAB
--0.1 -0.1      ! Original point for 2D k plane
- 0.2  0.0      ! The first vector to define 2D k plane 
- 0.0  0.2      ! The second vector to define 2D k plane  for arc plots
-
-KPLANE_BULK
- 0.00  0.00  0.50   ! Original point for 3D k plane 
- 1.00  0.00  0.00   ! The first vector to define 3d k space plane
- 0.00  0.50  0.00   ! The second vector to define 3d k space plane
-
-
-KCUBE_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
- 0.00  0.00  1.00   ! The third vector to define 3d k cube
-
-
-EFFECTIVE_MASS      ! optional
-2                   ! The i'th band to be calculated
-0.01                ! k step in unit of (1/Angstrom)
-0.0 0.0 0.0         ! k point where the effective mass calculated.
-
 KPOINTS_3D
 4 
 Direct
@@ -111,11 +64,6 @@ Direct
 0.0  0.5  0.0
 0.0  0.0  0.5
 
-SINGLEKPOINT_2D
-Direct 
-0.0 0.0
-
-
 WANNIER_CENTRES     ! copy from wannier90.wout
 Cartesian
  -0.000040  -1.194745   6.638646 
diff --git a/examples/Bi2Se3/wt.in-bands b/examples/Bi2Se3/wt.in-bands
new file mode 100644
index 00000000..ed4ade9e
--- /dev/null
+++ b/examples/Bi2Se3/wt.in-bands
@@ -0,0 +1,99 @@
+&TB_FILE
+Hrfile = 'wannier90_hr.dat'      
+Package = 'VASP'             ! obtained from VASP, it could be 'VASP', 'QE', 'Wien2k', 'OpenMx'
+/
+
+LATTICE
+Angstrom
+-2.069  -3.583614  0.000000     ! crystal lattice information
+ 2.069  -3.583614  0.000000
+ 0.000   2.389075  9.546667
+
+ATOM_POSITIONS
+5                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+ Bi 0.3990    0.3990    0.6970
+ Bi 0.6010    0.6010    0.3030
+ Se 0.0000    0.0000    0.5000
+ Se 0.2060    0.2060    0.1180
+ Se 0.7940    0.7940    0.8820
+
+PROJECTORS
+ 3 3 3 3 3          ! number of projectors
+Bi pz px py         ! projectors
+Bi pz px py
+Se pz px py
+Se pz px py
+Se pz px py
+
+!> bulk band structure calculation flag
+&CONTROL
+BulkBand_calc         = T
+BulkBand_points_calc  = T
+/
+
+&SYSTEM
+NumOccupied = 18        ! NumOccupied
+SOC = 1                 ! soc
+E_FERMI = 4.4195        ! e-fermi
+/
+
+! get projected bands onto different orbitals, here we only consider orbital and omit the spin freedom
+SELECTED_WANNIERORBITALS
+2 
+1-6                ! Bi
+7-15               ! Se
+
+
+&PARAMETERS
+Nk1 = 41            ! number k points  odd number would be better
+/
+
+KPATH_BULK            ! k point path
+4              ! number of k line only for bulk band
+G 0.00000 0.00000 0.0000 Z 0.00000 0.00000 0.5000
+Z 0.00000 0.00000 0.5000 F 0.50000 0.50000 0.0000
+F 0.50000 0.50000 0.0000 G 0.00000 0.00000 0.0000
+G 0.00000 0.00000 0.0000 L 0.50000 0.00000 0.0000  
+
+KPOINTS_3D
+4 
+Direct
+0.0  0.0  0.0
+0.5  0.0  0.0
+0.0  0.5  0.0
+0.0  0.0  0.5
+
+WANNIER_CENTRES     ! copy from wannier90.wout
+Cartesian
+ -0.000040  -1.194745   6.638646 
+  0.000038  -1.196699   6.640059 
+ -0.000032  -1.192363   6.640243 
+ -0.000086  -3.583414   2.908040 
+  0.000047  -3.581457   2.906587 
+ -0.000033  -3.585864   2.906443 
+ -0.000001   1.194527   4.773338 
+  0.000003   1.194538   4.773336 
+ -0.000037   1.194536   4.773327 
+  0.000006  -1.194384   1.130261 
+ -0.000018  -1.216986   1.140267 
+  0.000007  -1.172216   1.140684 
+  0.000011  -3.583770   8.416406 
+ -0.000002  -3.561169   8.406398 
+ -0.000007  -3.605960   8.405979 
+  0.000086  -1.194737   6.638626 
+ -0.000047  -1.196693   6.640080 
+  0.000033  -1.192286   6.640223 
+  0.000040  -3.583406   2.908021 
+ -0.000038  -3.581452   2.906608 
+  0.000032  -3.585788   2.906424 
+  0.000001   1.194548   4.773330 
+ -0.000003   1.194537   4.773332 
+  0.000037   1.194539   4.773340 
+ -0.000011  -1.194381   1.130260 
+  0.000002  -1.216981   1.140268 
+  0.000007  -1.172191   1.140687 
+ -0.000006  -3.583766   8.416405 
+  0.000018  -3.561165   8.406400 
+ -0.000007  -3.605935   8.405982 
+
diff --git a/examples/Graphene/wt.in-bands b/examples/Graphene/wt.in-bands
new file mode 100644
index 00000000..41dd63a4
--- /dev/null
+++ b/examples/Graphene/wt.in-bands
@@ -0,0 +1,42 @@
+&TB_FILE
+Hrfile = 'Graphene_hr.dat'
+/
+
+
+!> bulk band structure calculation flag
+&CONTROL
+BulkBand_calc         = T
+/
+
+&SYSTEM
+SOC = 0                 ! soc
+E_FERMI =  -1.2533        ! e-fermi
+/
+
+&PARAMETERS
+Nk1 = 101          ! number k points 
+/
+
+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
+C 0.666667 0.333333 0.500000 C
+
+PROJECTORS
+1 1        ! number of projectors
+C  pz
+C  pz
+
+KPATH_BULK            ! k point path
+3              ! number of k line only for bulk band
+  M   0.50000  0.00000  0.00000   K   0.33333333  0.33333333  0.00000 
+  K   0.33333333  0.33333333  0.00000   G   0.00000  0.00000  0.00000 
+  G   0.00000  0.00000  0.00000   M   0.50000  0.00000  0.00000 
+
diff --git a/examples/Graphene/wt.in-bands-plane b/examples/Graphene/wt.in-bands-plane
new file mode 100644
index 00000000..d898efb3
--- /dev/null
+++ b/examples/Graphene/wt.in-bands-plane
@@ -0,0 +1,43 @@
+&TB_FILE
+Hrfile = 'Graphene_hr.dat'
+/
+
+!> bulk band structure calculation flag
+&CONTROL
+BulkBand_plane_calc = T
+/
+
+&SYSTEM
+NumOccupied = 1
+SOC = 0                 ! soc
+E_FERMI =  -1.2533      ! e-fermi
+/
+
+
+&PARAMETERS
+Nk1 = 201          ! number k points 
+Nk2 = 201          ! number k points 
+/
+
+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
+C 0.666667 0.333333 0.500000 C
+
+PROJECTORS
+1 1        ! number of projectors
+C  pz
+C  pz
+
+KPLANE_BULK   ! fractional coordinates
+   0.333333  0.333333  0.000000   ! Middle point for a k slice(plane) in 3D BZ. Usually, the position of Dirac points.
+   0.100000  0.000000  0.000000   ! The first vector to define k plane(slice) in 3D BZ
+   0.000000  0.100000  0.000000   ! The second vector to define k plane(slice) in 3D BZ
+
diff --git a/examples/Graphene/wt.in-band b/examples/Graphene/wt.in-fsplane
similarity index 61%
rename from examples/Graphene/wt.in-band
rename to examples/Graphene/wt.in-fsplane
index 82865ca3..2c5597dd 100644
--- a/examples/Graphene/wt.in-band
+++ b/examples/Graphene/wt.in-fsplane
@@ -5,13 +5,10 @@ Hrfile = 'Graphene_hr.dat'
 
 !> bulk band structure calculation flag
 &CONTROL
-BulkBand_calc         = T
 BulkFS_Plane_calc     = T
 /
 
 &SYSTEM
-NSLAB = 100 
-NumOccupied = 1         ! NumOccupied
 SOC = 0                 ! soc
 E_FERMI =  -1.2533        ! e-fermi
 /
@@ -19,9 +16,6 @@ E_FERMI =  -1.2533        ! e-fermi
 &PARAMETERS
 Eta_Arc = 0.05     ! infinite small value, like brodening 
 E_arc = -1.00         ! energy for calculate Fermi Arc
-OmegaNum = 2001     ! omega number       
-OmegaMin = -1.      ! energy interval
-OmegaMax =  1.0     ! energy interval
 Nk1 = 101          ! number k points 
 Nk2 = 101          ! number k points 
 /
@@ -43,19 +37,8 @@ PROJECTORS
 C  pz
 C  pz
 
-SURFACE          
- 0  0  1   ! magnetic field direction
- 1  0  0
- 0  1  0
-
-KPATH_BULK            ! k point path
-3              ! number of k line only for bulk band
-  M   0.50000  0.00000  0.00000   K   0.33333333  0.33333333  0.00000 
-  K   0.33333333  0.33333333  0.00000   G   0.00000  0.00000  0.00000 
-  G   0.00000  0.00000  0.00000   M   0.50000  0.00000  0.00000 
-
 KPLANE_BULK  ! in fractional coordinates
-    0.00  0.00  0.00   ! Original point for 3D k plane 
+    0.00  0.00  0.00   ! center point for 3D k plane 
     2.00  0.00  0.00   ! The first vector to define 3d k space plane
     0.00  2.00  0.00   ! The second vector to define 3d k space plane
  
diff --git a/src/ek_bulk.f90 b/src/ek_bulk.f90
index 45a313e9..074ca2df 100644
--- a/src/ek_bulk.f90
+++ b/src/ek_bulk.f90
@@ -750,8 +750,9 @@ subroutine ek_bulk_plane
       write(outfileindex, '(a)')'#set xtics font ",24"'
       write(outfileindex, '(a)')'#set ytics font ",24"'
       write(outfileindex, '(a)')'set size ratio -1'
-      write(outfileindex, '(a)')'set xtics'
-      write(outfileindex, '(a)')'set ytics'
+      write(outfileindex, '(a)')'set ticslevel 0'
+      write(outfileindex, '(a)')'unset xtics'
+      write(outfileindex, '(a)')'unset ytics'
       write(outfileindex, '(a)')'set view 80,60'
       write(outfileindex, '(a)')'set xlabel "k_1"'
       write(outfileindex, '(a)')'set ylabel "k_2"'
diff --git a/src/main.f90 b/src/main.f90
index 9dec669e..82633ec5 100644
--- a/src/main.f90
+++ b/src/main.f90
@@ -222,6 +222,7 @@ program main
 #endif
       else
          call ek_bulk_line
+        !call ek_bulk_spin
         !call ek_bulk_mirror_z
       end if
       call now(time_end)
@@ -464,7 +465,6 @@ program main
        !call psik_bulk
        !call ek_bulk_polar
        !call ek_bulk_fortomas
-       !call ek_bulk_spin
        !call ek_bulk2D
        !call ek_bulk2D_spin
         call gapshape