diff --git a/examples/BHZ-model/wt.in b/examples/BHZ-model/wt.in
index 8a355c34..d2196be8 100644
--- a/examples/BHZ-model/wt.in
+++ b/examples/BHZ-model/wt.in
@@ -7,6 +7,7 @@ Hrfile = "BHZ_hr.dat"
 &CONTROL
 AHC_calc              = T
 BulkBand_calc         = T
+SlabBand_calc = T
 /
 
 &SYSTEM
@@ -44,6 +45,13 @@ PROJECTORS
 2            ! number of projectors
 C s pz
 
+KPATH_SLAB
+3        ! numker of k line for 2D case
+-X 0. -0.5 G 0.0 0.0  ! k path for 2D case
+G 0.0 0.0 X 0.0 0.5
+X 0.0 0.50 M 0.5 0.5  ! k path for 2D case
+
+
 
 SURFACE            ! See doc for details
  0  0  1
diff --git a/examples/Half-BHZ-model/Half_BHZ_hr_gen-case1.py b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case1.py
new file mode 100644
index 00000000..c5e5404e
--- /dev/null
+++ b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case1.py
@@ -0,0 +1,129 @@
+#!/bin/python3
+import numpy as np
+import cmath
+
+# The Hamiltonian is 
+#     ( M-Bk^2    Delta_0+A*k+  ) 
+#     ( Delta_0+A*k_    -M+Bk^2 )
+# where k^2=kx^2+ky^2
+
+# phase II, trivial insulator
+# A=0, M*B<0
+
+# phase III, trivial insulator with band inversion
+# A=0, M*B>0
+
+# phase I, trivial insulator
+# Delta_0=0, M*B<0
+
+# phase IV, Chern insulator with band inversion
+# Delta_0=0, M*B>0
+
+
+# from the kp to TB we use sustitution
+# k->sin(k)
+# k^2->2(1-cos(k))
+
+# Constants
+dp = np.float64
+pi = np.arctan(1) * 4
+zi = 1j
+
+# Lattice constants
+M = 2.0
+B =-1.0
+A = 0.0
+Delta_0=  1.0
+             
+
+# Number of Wannier functions and R points
+num_wann = 4
+nrpts = 7
+
+# R coordinates
+Irvec = np.zeros((3, nrpts), dtype=int)
+
+# Hamiltonian m,n are band indexes
+HmnR = np.zeros((num_wann, num_wann, nrpts), dtype=complex)
+
+# No of degeneracy of R point
+ndegen = np.ones(nrpts, dtype=int)
+
+# Initialization of matrices
+Irvec[:, :] = 0
+HmnR[:, :, :] = 0.0
+
+# 0 0 0
+ir = 0
+Irvec[:, ir] = [0, 0, 0]
+HmnR[0, 0, ir] = M - 4 * B
+HmnR[1, 1, ir] = -M + 4 * B
+HmnR[2, 2, ir] = M - 4 * B
+HmnR[3, 3, ir] = -M + 4 * B
+HmnR[0, 1, ir] = Delta_0
+HmnR[1, 0, ir] = Delta_0
+HmnR[2, 3, ir] = Delta_0
+HmnR[3, 2, ir] = Delta_0
+
+# 1 0
+ir = 1
+Irvec[:, ir] = [1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir] =-0.5*zi*A
+HmnR[1, 0, ir] =-0.5*zi*A
+HmnR[2, 3, ir] = 0.5*zi*A
+HmnR[3, 2, ir] = 0.5*zi*A
+
+# 0 1
+ir = 2
+Irvec[:, ir] = [0, 1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=  -A/2
+HmnR[1, 0, ir]=   A/2
+HmnR[2, 3, ir]=  -A/2
+HmnR[3, 2, ir]=   A/2
+
+
+# -1 0
+ir = 3
+Irvec[:, ir] = [-1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]= 0.5*zi*A
+HmnR[1, 0, ir]= 0.5*zi*A
+HmnR[2, 3, ir]=-0.5*zi*A
+HmnR[3, 2, ir]=-0.5*zi*A
+
+
+# 0 -1
+ir = 4
+Irvec[:, ir] = [0, -1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=   A/2
+HmnR[1, 0, ir]=  -A/2
+HmnR[2, 3, ir]=   A/2
+HmnR[3, 2, ir]=  -A/2
+
+nrpts= ir+1
+# Writing to a file
+with open('HalfBHZ_hr.dat', 'w') as file:
+    file.write('2-band half of BHZ model\n')
+    file.write('2\n')
+    file.write(f'{nrpts}\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(2):
+            for j in range(2):
+                file.write(f"{Irvec[0, ir]:5d}{Irvec[1, ir]:5d}{Irvec[2, ir]:5d}{i+1:5d}{j+1:5d} {HmnR[i, j, ir].real:16.8f} {HmnR[i, j, ir].imag:16.8f}\n")
+
diff --git a/examples/Half-BHZ-model/Half_BHZ_hr_gen-case2.py b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case2.py
new file mode 100644
index 00000000..13ff2afc
--- /dev/null
+++ b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case2.py
@@ -0,0 +1,129 @@
+#!/bin/python3
+import numpy as np
+import cmath
+
+# The Hamiltonian is 
+#     ( M-Bk^2    Delta_0+A*k+  ) 
+#     ( Delta_0+A*k_    -M+Bk^2 )
+# where k^2=kx^2+ky^2
+
+# phase II, trivial insulator
+# A=0, M*B<0
+
+# phase III, trivial insulator with band inversion
+# A=0, M*B>0
+
+# phase I, trivial insulator
+# Delta_0=0, M*B<0
+
+# phase IV, Chern insulator with band inversion
+# Delta_0=0, M*B>0
+
+
+# from the kp to TB we use sustitution
+# k->sin(k)
+# k^2->2(1-cos(k))
+
+# Constants
+dp = np.float64
+pi = np.arctan(1) * 4
+zi = 1j
+
+# Lattice constants
+M = 2.0
+B = 1.0
+A = 0.0
+Delta_0=  1.0
+             
+
+# Number of Wannier functions and R points
+num_wann = 4
+nrpts = 7
+
+# R coordinates
+Irvec = np.zeros((3, nrpts), dtype=int)
+
+# Hamiltonian m,n are band indexes
+HmnR = np.zeros((num_wann, num_wann, nrpts), dtype=complex)
+
+# No of degeneracy of R point
+ndegen = np.ones(nrpts, dtype=int)
+
+# Initialization of matrices
+Irvec[:, :] = 0
+HmnR[:, :, :] = 0.0
+
+# 0 0 0
+ir = 0
+Irvec[:, ir] = [0, 0, 0]
+HmnR[0, 0, ir] = M - 4 * B
+HmnR[1, 1, ir] = -M + 4 * B
+HmnR[2, 2, ir] = M - 4 * B
+HmnR[3, 3, ir] = -M + 4 * B
+HmnR[0, 1, ir] = Delta_0
+HmnR[1, 0, ir] = Delta_0
+HmnR[2, 3, ir] = Delta_0
+HmnR[3, 2, ir] = Delta_0
+
+# 1 0
+ir = 1
+Irvec[:, ir] = [1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir] =-0.5*zi*A
+HmnR[1, 0, ir] =-0.5*zi*A
+HmnR[2, 3, ir] = 0.5*zi*A
+HmnR[3, 2, ir] = 0.5*zi*A
+
+# 0 1
+ir = 2
+Irvec[:, ir] = [0, 1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=  -A/2
+HmnR[1, 0, ir]=   A/2
+HmnR[2, 3, ir]=  -A/2
+HmnR[3, 2, ir]=   A/2
+
+
+# -1 0
+ir = 3
+Irvec[:, ir] = [-1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]= 0.5*zi*A
+HmnR[1, 0, ir]= 0.5*zi*A
+HmnR[2, 3, ir]=-0.5*zi*A
+HmnR[3, 2, ir]=-0.5*zi*A
+
+
+# 0 -1
+ir = 4
+Irvec[:, ir] = [0, -1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=   A/2
+HmnR[1, 0, ir]=  -A/2
+HmnR[2, 3, ir]=   A/2
+HmnR[3, 2, ir]=  -A/2
+
+nrpts= ir+1
+# Writing to a file
+with open('HalfBHZ_hr.dat', 'w') as file:
+    file.write('2-band half of BHZ model\n')
+    file.write('2\n')
+    file.write(f'{nrpts}\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(2):
+            for j in range(2):
+                file.write(f"{Irvec[0, ir]:5d}{Irvec[1, ir]:5d}{Irvec[2, ir]:5d}{i+1:5d}{j+1:5d} {HmnR[i, j, ir].real:16.8f} {HmnR[i, j, ir].imag:16.8f}\n")
+
diff --git a/examples/Half-BHZ-model/Half_BHZ_hr_gen-case3.py b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case3.py
new file mode 100644
index 00000000..b544f547
--- /dev/null
+++ b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case3.py
@@ -0,0 +1,129 @@
+#!/bin/python3
+import numpy as np
+import cmath
+
+# The Hamiltonian is 
+#     ( M-Bk^2    Delta_0+A*k+  ) 
+#     ( Delta_0+A*k_    -M+Bk^2 )
+# where k^2=kx^2+ky^2
+
+# phase II, trivial insulator
+# A=0, M*B<0
+
+# phase III, trivial insulator with band inversion
+# A=0, M*B>0
+
+# phase I, trivial insulator
+# Delta_0=0, M*B<0
+
+# phase IV, Chern insulator with band inversion
+# Delta_0=0, M*B>0
+
+
+# from the kp to TB we use sustitution
+# k->sin(k)
+# k^2->2(1-cos(k))
+
+# Constants
+dp = np.float64
+pi = np.arctan(1) * 4
+zi = 1j
+
+# Lattice constants
+M = 2.0
+B =-1.0
+A = 1.0
+Delta_0=  0.0
+             
+
+# Number of Wannier functions and R points
+num_wann = 4
+nrpts = 7
+
+# R coordinates
+Irvec = np.zeros((3, nrpts), dtype=int)
+
+# Hamiltonian m,n are band indexes
+HmnR = np.zeros((num_wann, num_wann, nrpts), dtype=complex)
+
+# No of degeneracy of R point
+ndegen = np.ones(nrpts, dtype=int)
+
+# Initialization of matrices
+Irvec[:, :] = 0
+HmnR[:, :, :] = 0.0
+
+# 0 0 0
+ir = 0
+Irvec[:, ir] = [0, 0, 0]
+HmnR[0, 0, ir] = M - 4 * B
+HmnR[1, 1, ir] = -M + 4 * B
+HmnR[2, 2, ir] = M - 4 * B
+HmnR[3, 3, ir] = -M + 4 * B
+HmnR[0, 1, ir] = Delta_0
+HmnR[1, 0, ir] = Delta_0
+HmnR[2, 3, ir] = Delta_0
+HmnR[3, 2, ir] = Delta_0
+
+# 1 0
+ir = 1
+Irvec[:, ir] = [1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir] =-0.5*zi*A
+HmnR[1, 0, ir] =-0.5*zi*A
+HmnR[2, 3, ir] = 0.5*zi*A
+HmnR[3, 2, ir] = 0.5*zi*A
+
+# 0 1
+ir = 2
+Irvec[:, ir] = [0, 1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=  -A/2
+HmnR[1, 0, ir]=   A/2
+HmnR[2, 3, ir]=  -A/2
+HmnR[3, 2, ir]=   A/2
+
+
+# -1 0
+ir = 3
+Irvec[:, ir] = [-1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]= 0.5*zi*A
+HmnR[1, 0, ir]= 0.5*zi*A
+HmnR[2, 3, ir]=-0.5*zi*A
+HmnR[3, 2, ir]=-0.5*zi*A
+
+
+# 0 -1
+ir = 4
+Irvec[:, ir] = [0, -1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=   A/2
+HmnR[1, 0, ir]=  -A/2
+HmnR[2, 3, ir]=   A/2
+HmnR[3, 2, ir]=  -A/2
+
+nrpts= ir+1
+# Writing to a file
+with open('HalfBHZ_hr.dat', 'w') as file:
+    file.write('2-band half of BHZ model\n')
+    file.write('2\n')
+    file.write(f'{nrpts}\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(2):
+            for j in range(2):
+                file.write(f"{Irvec[0, ir]:5d}{Irvec[1, ir]:5d}{Irvec[2, ir]:5d}{i+1:5d}{j+1:5d} {HmnR[i, j, ir].real:16.8f} {HmnR[i, j, ir].imag:16.8f}\n")
+
diff --git a/examples/Half-BHZ-model/Half_BHZ_hr_gen-case4.py b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case4.py
new file mode 100644
index 00000000..661cec8a
--- /dev/null
+++ b/examples/Half-BHZ-model/Half_BHZ_hr_gen-case4.py
@@ -0,0 +1,129 @@
+#!/bin/python3
+import numpy as np
+import cmath
+
+# The Hamiltonian is 
+#     ( M-Bk^2    Delta_0+A*k+  ) 
+#     ( Delta_0+A*k_    -M+Bk^2 )
+# where k^2=kx^2+ky^2
+
+# phase II, trivial insulator
+# A=0, M*B<0
+
+# phase III, trivial insulator with band inversion
+# A=0, M*B>0
+
+# phase I, trivial insulator
+# Delta_0=0, M*B<0
+
+# phase IV, Chern insulator with band inversion
+# Delta_0=0, M*B>0
+
+
+# from the kp to TB we use sustitution
+# k->sin(k)
+# k^2->2(1-cos(k))
+
+# Constants
+dp = np.float64
+pi = np.arctan(1) * 4
+zi = 1j
+
+# Lattice constants
+M = 2.0
+B = 1.0
+A = 1.0
+Delta_0=  0.0
+             
+
+# Number of Wannier functions and R points
+num_wann = 4
+nrpts = 7
+
+# R coordinates
+Irvec = np.zeros((3, nrpts), dtype=int)
+
+# Hamiltonian m,n are band indexes
+HmnR = np.zeros((num_wann, num_wann, nrpts), dtype=complex)
+
+# No of degeneracy of R point
+ndegen = np.ones(nrpts, dtype=int)
+
+# Initialization of matrices
+Irvec[:, :] = 0
+HmnR[:, :, :] = 0.0
+
+# 0 0 0
+ir = 0
+Irvec[:, ir] = [0, 0, 0]
+HmnR[0, 0, ir] = M - 4 * B
+HmnR[1, 1, ir] = -M + 4 * B
+HmnR[2, 2, ir] = M - 4 * B
+HmnR[3, 3, ir] = -M + 4 * B
+HmnR[0, 1, ir] = Delta_0
+HmnR[1, 0, ir] = Delta_0
+HmnR[2, 3, ir] = Delta_0
+HmnR[3, 2, ir] = Delta_0
+
+# 1 0
+ir = 1
+Irvec[:, ir] = [1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir] =-0.5*zi*A
+HmnR[1, 0, ir] =-0.5*zi*A
+HmnR[2, 3, ir] = 0.5*zi*A
+HmnR[3, 2, ir] = 0.5*zi*A
+
+# 0 1
+ir = 2
+Irvec[:, ir] = [0, 1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=  -A/2
+HmnR[1, 0, ir]=   A/2
+HmnR[2, 3, ir]=  -A/2
+HmnR[3, 2, ir]=   A/2
+
+
+# -1 0
+ir = 3
+Irvec[:, ir] = [-1, 0, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]= 0.5*zi*A
+HmnR[1, 0, ir]= 0.5*zi*A
+HmnR[2, 3, ir]=-0.5*zi*A
+HmnR[3, 2, ir]=-0.5*zi*A
+
+
+# 0 -1
+ir = 4
+Irvec[:, ir] = [0, -1, 0]
+HmnR[0, 0, ir] = B
+HmnR[1, 1, ir] = -B
+HmnR[2, 2, ir] = B
+HmnR[3, 3, ir] = -B
+HmnR[0, 1, ir]=   A/2
+HmnR[1, 0, ir]=  -A/2
+HmnR[2, 3, ir]=   A/2
+HmnR[3, 2, ir]=  -A/2
+
+nrpts= ir+1
+# Writing to a file
+with open('HalfBHZ_hr.dat', 'w') as file:
+    file.write('2-band half of BHZ model\n')
+    file.write('2\n')
+    file.write(f'{nrpts}\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(2):
+            for j in range(2):
+                file.write(f"{Irvec[0, ir]:5d}{Irvec[1, ir]:5d}{Irvec[2, ir]:5d}{i+1:5d}{j+1:5d} {HmnR[i, j, ir].real:16.8f} {HmnR[i, j, ir].imag:16.8f}\n")
+
diff --git a/examples/Half-BHZ-model/readme.txt b/examples/Half-BHZ-model/readme.txt
new file mode 100644
index 00000000..ed6b2c3a
--- /dev/null
+++ b/examples/Half-BHZ-model/readme.txt
@@ -0,0 +1,21 @@
+# An example to show how does the band inversion and hybridization affect the band topology
+# we prepared four python scripts and one wt.in file to calculate the band structure, slab bands
+# AHC, Wilson loop, Berry curvature
+
+# to run it, first, you need to generate the wannier TB model HalfBHZ_hr.dat
+python3 Half_BHZ_hr_gen-case1.py
+
+# 2nd run wanniertools wt.x
+wt.x 
+
+# 3rd get the plots
+
+gnuplot bulkek.gnu
+gnuplot slabek.gnu
+gnuplot wcc.gnu
+gnuplot sigma_ahc.gnu
+gnuplot Berrycurvature.gnu
+
+# use ll -tr to check the latest files
+ll -tr
+
diff --git a/examples/Half-BHZ-model/wt.in b/examples/Half-BHZ-model/wt.in
new file mode 100644
index 00000000..ae25fa3f
--- /dev/null
+++ b/examples/Half-BHZ-model/wt.in
@@ -0,0 +1,81 @@
+&TB_FILE
+Hrfile = "HalfBHZ_hr.dat"
+/
+
+!> bulk band structure calculation flag
+&CONTROL
+AHC_calc              = T
+BulkBand_calc         = T
+BulkBand_plane_calc   = F
+SlabBand_calc         = T
+FindNodes_calc        = F
+SlabSS_calc           = F
+SlabArc_calc          = F
+Wanniercenter_calc    = T
+LandauLevel_B_calc    = F
+LandauLevel_B_dos_calc = F
+BerryCurvature_calc   = T
+/
+
+&SYSTEM
+NSLAB =40
+NumOccupied = 1         ! NumOccupied
+SOC = 0                 ! soc
+E_FERMI = 0        ! e-fermi
+/
+
+&PARAMETERS
+Fermi_broadening = 0.01     ! infinite small value, like brodening 
+iso_energy = 0.0         ! energy for calculate Fermi Arc
+OmegaNum = 400  ! omega number       
+OmegaMin = -4.0     ! energy interval
+OmegaMax =  4.0     ! energy interval
+Nk1 = 60            ! number k points 
+Nk2 = 60            ! number k points 
+Nk3 = 1             ! number k points 
+NP = 2              ! number of principle layers
+/
+
+LATTICE
+Angstrom
+3 0 0 
+0 3 0
+0 0 10
+
+ATOM_POSITIONS
+1                               ! number of atoms for projectors
+Direct                          ! Direct or Cartisen coordinate
+C 0 0 0
+
+PROJECTORS
+2            ! number of projectors
+C s pz
+
+
+SURFACE            ! See doc for details
+ 0  0  1
+ 1  0  0
+ 0  1  0
+
+KPATH_BULK            ! k point path
+2              ! number of k line only for bulk band
+  X   0.50000  0.00000  0.00000   G  0.00000   0.00000  0.00000
+  G   0.00000  0.00000  0.00000   Y  0.00000   0.50000  0.00000
+
+KPATH_SLAB
+2        ! numker of k line for 2D case
+-X 0. -0.5 G 0.0 0.0  ! k path for 2D case
+G 0.0 0.0 X 0.0 0.5
+
+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
+
+
+KCUBE_BULK
+-0.50 -0.50 -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  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
+
diff --git a/examples/MoS2-1Tp/wt.in b/examples/MoS2-1Tp/wt.in
index fa5031e1..e5725823 100644
--- a/examples/MoS2-1Tp/wt.in
+++ b/examples/MoS2-1Tp/wt.in
@@ -6,15 +6,15 @@ Package = 'VASP'
 
 !> bulk band structure calculation flag
 &CONTROL
-BulkBand_calc         = T
+BulkBand_calc         = F
 DOS_calc              = F
 SlabBand_calc         = F
 SlabSS_calc           = T
-wanniercenter_calc    = T
+wanniercenter_calc    = F
 /
 
 &SYSTEM
-NSLAB = 20 
+NSLAB = 10 
 NumOccupied = 28        ! NumOccupied
 SOC = 1                 ! soc
 E_FERMI =  -3.60        ! e-fermi
@@ -56,8 +56,8 @@ S pz px py
 
 SURFACE            ! MoS2 conventional (010) surface
  0  1  0
- 0  0  1
  1  0  0
+ 0  0  1
 
 KPATH_BULK            ! k point path
 4              ! number of k line only for bulk band
diff --git a/examples/WSe2-1Tp/wt.in b/examples/WSe2-1Tp/wt.in
index e0c8977f..5fb66ac2 100644
--- a/examples/WSe2-1Tp/wt.in
+++ b/examples/WSe2-1Tp/wt.in
@@ -6,9 +6,9 @@ Package = 'VASP'
 
 !> bulk band structure calculation flag
 &CONTROL
-BulkBand_calc         = T
+BulkBand_calc         = F
 SlabBand_calc         = F
-SlabSS_calc           = T
+SlabSS_calc           = F
 wanniercenter_calc    = T
 /
 
@@ -25,7 +25,7 @@ OmegaMin = -1.0     ! energy interval
 OmegaMax =  0.5     ! energy interval
 Nk1 = 61          ! number k points 
 Nk2 = 101          ! number k points 
-NP = 1              ! number of principle layers
+NP = 2              ! number of principle layers
 /
 
 LATTICE
diff --git a/src/wanniercenter.f90 b/src/wanniercenter.f90
index 067aa6c6..8da38c6d 100644
--- a/src/wanniercenter.f90
+++ b/src/wanniercenter.f90
@@ -2389,7 +2389,7 @@ subroutine  wannier_center3D_plane0
       write(outfileindex, '(a)')'set ytics format "%4.1f" nomirror out'
       write(outfileindex, '(a)')'set ylabel "WCC"'
       write(outfileindex, '(a)')'set ylabel offset 2, 0.0 '
-      write(outfileindex, '(a)')'set xrange [0: 0.5]'
+      write(outfileindex, '(a)')'set xrange [0: 1.0]'
       write(outfileindex, '(a)')'set yrange [0:1]'
       write(outfileindex, '(a, i5, a)')"plot for [i=4: ", NumberofSelectedOccupiedBands+3, &
          "] 'wcc.dat' u 1:i w p  pt 7  ps 0.5 lc 'black'"
diff --git a/src/wanniercenter_adaptive.f90 b/src/wanniercenter_adaptive.f90
index c5f86823..c23169ca 100644
--- a/src/wanniercenter_adaptive.f90
+++ b/src/wanniercenter_adaptive.f90
@@ -332,7 +332,7 @@ subroutine  wannier_center3D_plane_adaptive
          write(outfileindex, '(a)')'set ytics format "%4.1f" nomirror out'
          write(outfileindex, '(a)')'set ylabel "WCC"'
          write(outfileindex, '(a)')'set ylabel offset 2, 0.0 '
-         write(outfileindex, '(a)')'set xrange [0: 0.5]'
+         write(outfileindex, '(a)')'set xrange [0: 1.0]'
          write(outfileindex, '(a)')'set yrange [0:1]'
          write(outfileindex, '(a, i5, a)')"plot for [i=4: ", NumberofSelectedOccupiedBands+3, &
             "] 'wcc.dat' u 1:i w p  pt 7  ps 0.5 lc 'black'"