diff --git a/examples/BHZ-model/BHZ_hr_gen-case1.py b/examples/BHZ-model/BHZ_hr_gen-case1.py
new file mode 100644
index 00000000..e4455315
--- /dev/null
+++ b/examples/BHZ-model/BHZ_hr_gen-case1.py
@@ -0,0 +1,123 @@
+#!/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
+
+# Case I, QSHE with band inversion and no trivial hybridization
+# Delta_0=0, M*B>0, |A|>0
+
+# Case I, QSHE with band inversion and with trivial and nontrivial hybridization
+# Delta_0=0.5, M*B>0, |A|>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('BHZ_hr.dat', 'w') as file:
+    file.write('4-band of BHZ model\n')
+    file.write('4 !num_wann \n')
+    file.write(f'{nrpts} ! nrpts\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(4):
+            for j in range(4):
+                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/BHZ-model/BHZ_hr_gen-case2.py b/examples/BHZ-model/BHZ_hr_gen-case2.py
new file mode 100644
index 00000000..f0a68739
--- /dev/null
+++ b/examples/BHZ-model/BHZ_hr_gen-case2.py
@@ -0,0 +1,123 @@
+#!/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
+
+# Case I, QSHE with band inversion and no trivial hybridization
+# Delta_0=0, M*B>0, |A|>0
+
+# Case I, QSHE with band inversion and with trivial and nontrivial hybridization
+# Delta_0=0.5, M*B>0, |A|>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.5
+             
+
+# 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('BHZ_hr.dat', 'w') as file:
+    file.write('4-band of BHZ model\n')
+    file.write('4 !num_wann \n')
+    file.write(f'{nrpts} ! nrpts\n')
+    file.write(' '.join(f'{x:5d}' for x in ndegen) + '\n')
+    for ir in range(nrpts):
+        for i in range(4):
+            for j in range(4):
+                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/BHZ-model/wt.in b/examples/BHZ-model/wt.in
index a3304969..a4bd8c3a 100644
--- a/examples/BHZ-model/wt.in
+++ b/examples/BHZ-model/wt.in
@@ -7,6 +7,7 @@ Hrfile = "BHZ_hr.dat"
 BulkBand_calc   = T   ! bulk band structure
 SlabBand_calc   = T   ! slab band structure
 SHC_calc        = T   ! spin Hall conductivity
+AHC_calc        = T   ! anomalous Hall conductivity
 Wilsonloop_calc = T   ! Wilson loop
 /
 
@@ -14,6 +15,9 @@ Wilsonloop_calc = T   ! Wilson loop
 SOC = 1                 ! soc
 E_FERMI = 0        ! e-fermi
 NumOccupied = 2
+Bx= 0, By= 0, Bz= 1     ! Bx By Bz
+Add_Zeeman_Field = T
+Zeeman_energy_in_eV = 0.5 ! in eV
 /
 
 &PARAMETERS
@@ -58,6 +62,11 @@ KPATH_BULK            ! k point path
   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.50000
 
+KPLANE_BULK
+0 0 0
+1 0 0
+0 0.5 0
+
 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
diff --git a/examples/BHZ-model/wt.in-normal b/examples/BHZ-model/wt.in-normal
new file mode 100644
index 00000000..a3304969
--- /dev/null
+++ b/examples/BHZ-model/wt.in-normal
@@ -0,0 +1,66 @@
+&TB_FILE
+Hrfile = "BHZ_hr.dat"
+/
+
+!>  flags to control different functionalities
+&CONTROL
+BulkBand_calc   = T   ! bulk band structure
+SlabBand_calc   = T   ! slab band structure
+SHC_calc        = T   ! spin Hall conductivity
+Wilsonloop_calc = T   ! Wilson loop
+/
+
+&SYSTEM
+SOC = 1                 ! soc
+E_FERMI = 0        ! e-fermi
+NumOccupied = 2
+/
+
+&PARAMETERS
+Fermi_broadening = 0.001     ! infinite small value, like brodening 
+iso_energy = 0.0         ! energy for calculate Fermi Arc
+OmegaNum = 400  ! omega number       
+OmegaMin = -8.0     ! energy interval
+OmegaMax =  8.0     ! energy interval
+Nk1 = 60            ! number k points 
+Nk2 = 60            ! number k points 
+Nk3 = 2             ! number k points 
+/
+
+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
+
+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
+ 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.50000
+
+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/BHZ-model/wt.in-ahe-withzeeman b/examples/BHZ-model/wt.in-with-zeeman
similarity index 58%
rename from examples/BHZ-model/wt.in-ahe-withzeeman
rename to examples/BHZ-model/wt.in-with-zeeman
index a9336d9f..367283b2 100644
--- a/examples/BHZ-model/wt.in-ahe-withzeeman
+++ b/examples/BHZ-model/wt.in-with-zeeman
@@ -2,47 +2,55 @@
 Hrfile = "BHZ_hr.dat"
 /
 
-
-!> bulk band structure calculation flag
+!>  flags to control different functionalities
 &CONTROL
-AHC_calc              = T
-BulkBand_calc         = T
+BulkBand_calc   = T   ! bulk band structure
+SlabBand_calc   = T   ! slab band structure
+SHC_calc        = T   ! spin Hall conductivity
+AHC_calc        = T   ! anomalous Hall conductivity
+Wilsonloop_calc = T   ! Wilson loop
 /
 
 &SYSTEM
 SOC = 1                 ! soc
 E_FERMI = 0        ! e-fermi
+NumOccupied = 2
 Bx= 0, By= 0, Bz= 1     ! Bx By Bz
 Add_Zeeman_Field = T
 Zeeman_energy_in_eV = 0.5 ! in eV
 /
 
 &PARAMETERS
-Eta_Arc = 0.001     ! infinite small value, like brodening 
-E_arc = 0.0         ! energy for calculate Fermi Arc
+Fermi_broadening = 0.001     ! infinite small value, like brodening 
+iso_energy = 0.0         ! energy for calculate Fermi Arc
 OmegaNum = 400  ! omega number       
 OmegaMin = -8.0     ! energy interval
 OmegaMax =  8.0     ! energy interval
 Nk1 = 60            ! number k points 
 Nk2 = 60            ! number k points 
-Nk3 = 1             ! number k points 
+Nk3 = 2             ! number k points 
 /
 
 LATTICE
 Angstrom
-1 0 0 
-0 1 0
-0 0 1
+3 0 0 
+0 3 0
+0 0 10
 
 ATOM_POSITIONS
-1                               ! number of atoms for projectors
-Direct                          ! Direct or Cartisen coordinate
+1                       ! number of atoms for projectors
+Direct                  ! Direct or Cartisen coordinate
 C 0 0 0
 
 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
@@ -54,7 +62,6 @@ KPATH_BULK            ! k point path
   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.50000
 
-
 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
diff --git a/src/readinput.f90 b/src/readinput.f90
index 5218460a..c82f2a07 100644
--- a/src/readinput.f90
+++ b/src/readinput.f90
@@ -609,8 +609,6 @@ subroutine readinput
    polarization_alpha_arpes= (45d0/180d0)*3.14159265358979d0
    polarization_delta_arpes= (0d0/180d0)*3.14159265358979d0
 
-
-
    !> by default, we only project on atoms for a given wave function
    projection_weight_mode = "NORMAL"
 
@@ -627,7 +625,6 @@ subroutine readinput
       if (OmegaNum_unfold==0) OmegaNum_unfold= 200
    endif
 
-   if (NumLCZVecs> Num_wann) NumLCZVecs= Num_wann
 
    if (stat>0) then
 
@@ -1136,7 +1133,7 @@ subroutine readinput
       stop
    endif
 
-
+   if (NumLCZVecs> NumberOfspinorbitals) NumLCZVecs= NumberOfspinorbitals
 
 110 continue
 
diff --git a/src/wanniercenter_adaptive.f90 b/src/wanniercenter_adaptive.f90
index c23169ca..cf4d28b5 100644
--- a/src/wanniercenter_adaptive.f90
+++ b/src/wanniercenter_adaptive.f90
@@ -743,12 +743,26 @@ subroutine  Wcc_integrate_func(k2, kvec1, wcc, &
       VT= 0d0
 
       if (Is_Sparse) then
-         neval=OmegaNum
-         if (neval>Num_wann-2) neval= Num_wann- 2
+
+         if (OmegaNum==0) OmegaNum= Num_wann
+         if (NumSelectedEigenVals==0) NumSelectedEigenVals=OmegaNum
+      
+         neval= NumSelectedEigenVals
+      
+         if (neval>Num_wann-2) then
+            neval= Num_wann- 2
+            nvecs= Num_wann
+         endif
+      
+         !> ncv= NumLCZVecs if specfied in wt.in
+         if (NumLCZVecs.ne.0) then
+            nvecs= NumLCZVecs
+         else
+            nvecs=int(2*neval)
+            if (nvecs<50) nvecs= int(6*neval)
+         endif
       
-         !> ncv
-         nvecs=int(2*neval)
-         if (nvecs<50) nvecs= int(6*neval)
+         if (neval+2>=nvecs) neval= nvecs-2
       
          if (nvecs>Num_wann) nvecs= Num_wann