diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
index 422e2f9fef..d85e6e0f27 100644
--- a/source/source_esolver/esolver_ks_lcao.cpp
+++ b/source/source_esolver/esolver_ks_lcao.cpp
@@ -450,7 +450,9 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
                                                   PARAM.inp.ks_solver,
                                                   PARAM.globalv.kpar_lcao,
                                                   PARAM.globalv.nlocal,
-                                                  PARAM.inp.nelec);
+                                                  PARAM.inp.nbands,
+                                                  PARAM.inp.nelec,
+                                                  PARAM.inp.device == "gpu");
         hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, 
           this->chr, PARAM.inp.nspin, skip_charge);
     }
diff --git a/source/source_esolver/esolver_ks_lcao_tddft.cpp b/source/source_esolver/esolver_ks_lcao_tddft.cpp
index 154a3a65e2..f0b0e6de78 100644
--- a/source/source_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/source_esolver/esolver_ks_lcao_tddft.cpp
@@ -351,7 +351,9 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
                                                                         PARAM.inp.ks_solver,
                                                                         PARAM.globalv.kpar_lcao,
                                                                         PARAM.globalv.nlocal,
-                                                                        PARAM.inp.nelec);
+                                                                        PARAM.inp.nbands,
+                                                                        PARAM.inp.nelec,
+                                                                        PARAM.inp.device == "gpu");
             hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                    this->psi[0],
                                    this->pelec,
diff --git a/source/source_hsolver/diago_cusolver.cpp b/source/source_hsolver/diago_cusolver.cpp
index b02c413257..34ced44cb8 100644
--- a/source/source_hsolver/diago_cusolver.cpp
+++ b/source/source_hsolver/diago_cusolver.cpp
@@ -6,7 +6,6 @@
 #include "source_base/module_external/scalapack_connector.h"
 #include "source_base/tool_title.h"
 #include "source_base/timer.h"
-#include "source_io/module_parameter/parameter.h"
 
 #include <memory>
 #include <type_traits>
@@ -22,7 +21,7 @@ template <typename T>
 int DiagoCusolver<T>::DecomposedState = 0;
 
 template <typename T>
-DiagoCusolver<T>::DiagoCusolver()
+DiagoCusolver<T>::DiagoCusolver(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in)
 {
 }
 
@@ -42,13 +41,13 @@ void DiagoCusolver<T>::diag(
     ModuleBase::TITLE("DiagoCusolver", "diag");
     ModuleBase::timer::start("DiagoCusolver", "cusolver");
     // Allocate memory for eigenvalues
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(h_mat.row * h_mat.col);
     this->dc.Dngvd(h_mat.row, h_mat.col, h_mat.p, s_mat.p, eigen.data(), eigenvectors.data());
     const int size = psi.get_nbands() * psi.get_nbasis();
     BlasConnector::copy(size, eigenvectors.data(), 1, psi.get_pointer(), 1);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
     ModuleBase::timer::end("DiagoCusolver", "cusolver");
 }
 
diff --git a/source/source_hsolver/diago_cusolver.h b/source/source_hsolver/diago_cusolver.h
index 28ebdc852c..50d9a393c6 100644
--- a/source/source_hsolver/diago_cusolver.h
+++ b/source/source_hsolver/diago_cusolver.h
@@ -19,9 +19,11 @@ class DiagoCusolver
 
   public:
 
-    DiagoCusolver();
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    DiagoCusolver(const int nlocal_in, const int nbands_in);
     ~DiagoCusolver();
-    
+
     // Override the diag function for CUSOLVER diagonalization
     void diag(
       hamilt::MatrixBlock<T>& h_mat,
@@ -40,6 +42,9 @@ class DiagoCusolver
     // Function to check if ELPA handle needs to be created or reused in MPI settings
     bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const;
 #endif
+
+    const int nlocal;
+    const int nbands;
 };
 
 } // namespace hsolver
diff --git a/source/source_hsolver/diago_cusolvermp.cpp b/source/source_hsolver/diago_cusolvermp.cpp
index b5bb80eb5f..395add83d2 100644
--- a/source/source_hsolver/diago_cusolvermp.cpp
+++ b/source/source_hsolver/diago_cusolvermp.cpp
@@ -1,6 +1,5 @@
 #ifdef __CUSOLVERMP
 
-#include "source_io/module_parameter/parameter.h"
 #include "diago_cusolvermp.h"
 
 #include "source_base/module_external/blas_connector.h"
@@ -18,7 +17,7 @@ void DiagoCusolverMP<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real*
     hamilt::MatrixBlock<T> h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<Real> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(h_mat.row * h_mat.col);
 
     MPI_Comm COMM_DIAG = MPI_COMM_WORLD; // use all processes
@@ -30,7 +29,7 @@ void DiagoCusolverMP<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real*
         ModuleBase::timer::end("DiagoCusolverMP", "Diag_CusolverMP_gvd");
     }
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
     const int size = psi.get_nbands() * psi.get_nbasis();
     BlasConnector::copy(size, eigenvectors.data(), inc, psi.get_pointer(), inc);
 }
diff --git a/source/source_hsolver/diago_cusolvermp.h b/source/source_hsolver/diago_cusolvermp.h
index d173315e31..155a4bf94d 100644
--- a/source/source_hsolver/diago_cusolvermp.h
+++ b/source/source_hsolver/diago_cusolvermp.h
@@ -16,11 +16,17 @@ class DiagoCusolverMP
     using Real = typename GetTypeReal<T>::type;
 
   public:
-    DiagoCusolverMP()
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    DiagoCusolverMP(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in)
     {
     }
     // the diag function for CUSOLVERMP diagonalization
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
+
+  private:
+    const int nlocal;
+    const int nbands;
 };
 } // namespace hsolver
 #endif // __CUSOLVERMP
diff --git a/source/source_hsolver/diago_elpa.cpp b/source/source_hsolver/diago_elpa.cpp
index a264a41b9c..55ca982539 100644
--- a/source/source_hsolver/diago_elpa.cpp
+++ b/source/source_hsolver/diago_elpa.cpp
@@ -2,7 +2,6 @@
 #include "source_base/global_function.h"
 #include "source_base/module_external/blas_connector.h"
 
-#include "source_io/module_parameter/parameter.h"
 #include "module_genelpa/elpa_solver.h"
 #include "source_base/module_external/blacs_connector.h"
 #include "source_base/global_variable.h"
@@ -75,13 +74,13 @@ void DiagoElpa<std::complex<double>>::diag(
     matcd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
 
     bool isReal = false;
     MPI_Comm COMM_DIAG = setmpicomm(); // set mpi_comm needed
     ELPA_Solver es((const bool)isReal,
                    COMM_DIAG,
-                   (const int)PARAM.inp.nbands,
+                   (const int)this->nbands,
                    (const int)h_mat.row,
                    (const int)h_mat.col,
                    (const int*)h_mat.desc);
@@ -97,7 +96,7 @@ void DiagoElpa<std::complex<double>>::diag(
     es.exit();
 
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 #else
     ModuleBase::WARNING_QUIT("DiagoElpa",
                              "DiagoElpa only can be used with macro __MPI");
@@ -113,15 +112,13 @@ void DiagoElpa<double>::diag(hamilt::Hamilt<double>* phm_in,
     matd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
 
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
 
     bool isReal = true;
     MPI_Comm COMM_DIAG = setmpicomm(); // set mpi_comm needed
-    // ELPA_Solver es(isReal, COMM_DIAG, PARAM.inp.nbands, h_mat.row, h_mat.col,
-    // h_mat.desc);
     ELPA_Solver es((const bool)isReal,
                    COMM_DIAG,
-                   (const int)PARAM.inp.nbands,
+                   (const int)this->nbands,
                    (const int)h_mat.row,
                    (const int)h_mat.col,
                    (const int*)h_mat.desc);
@@ -135,7 +132,7 @@ void DiagoElpa<double>::diag(hamilt::Hamilt<double>* phm_in,
     es.exit();
 
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 #else
     ModuleBase::WARNING_QUIT("DiagoElpa",
                              "DiagoElpa only can be used with macro __MPI");
@@ -151,11 +148,11 @@ void DiagoElpa<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex
     Real* eigenvalue_in,
     MPI_Comm& comm)
 {
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     bool isReal = false;
     ELPA_Solver es((const bool)isReal,
                    comm,
-                   (const int)PARAM.inp.nbands,
+                   (const int)this->nbands,
                    (const int)h_mat.row,
                    (const int)h_mat.col,
                    (const int*)h_mat.desc);
@@ -170,7 +167,7 @@ void DiagoElpa<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex
     ModuleBase::timer::end("DiagoElpa", "elpa_solve");
     es.exit();
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 
 template <>
@@ -180,14 +177,12 @@ void DiagoElpa<double>::diag_pool(hamilt::MatrixBlock<double>& h_mat,
     Real* eigenvalue_in,
     MPI_Comm& comm)
 {
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
 
     bool isReal = true;
-    // ELPA_Solver es(isReal, COMM_DIAG, PARAM.inp.nbands, h_mat.row, h_mat.col,
-    // h_mat.desc);
     ELPA_Solver es((const bool)isReal,
                    comm,
-                   (const int)PARAM.inp.nbands,
+                   (const int)this->nbands,
                    (const int)h_mat.row,
                    (const int)h_mat.col,
                    (const int*)h_mat.desc);
@@ -203,7 +198,7 @@ void DiagoElpa<double>::diag_pool(hamilt::MatrixBlock<double>& h_mat,
     const int inc = 1;
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
                                 "K-S equation was solved by genelpa2");
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,
                                 "eigenvalues were copied to ekb");
 }
diff --git a/source/source_hsolver/diago_elpa.h b/source/source_hsolver/diago_elpa.h
index c6050af224..2bee7a14fa 100644
--- a/source/source_hsolver/diago_elpa.h
+++ b/source/source_hsolver/diago_elpa.h
@@ -15,6 +15,10 @@ class DiagoElpa
     using Real = typename GetTypeReal<T>::type;
 
   public:
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    DiagoElpa(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {};
+
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
@@ -30,6 +34,9 @@ class DiagoElpa
     bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const;
     static int lastmpinum; // last using mpi;
 #endif
+
+    const int nlocal;
+    const int nbands;
 };
 
 template <typename T>
diff --git a/source/source_hsolver/diago_elpa_native.cpp b/source/source_hsolver/diago_elpa_native.cpp
index ee003a08e0..fb4cfbc448 100644
--- a/source/source_hsolver/diago_elpa_native.cpp
+++ b/source/source_hsolver/diago_elpa_native.cpp
@@ -4,7 +4,6 @@
 #include "source_base/module_external/blas_connector.h"
 #include "source_base/module_external/blacs_connector.h"
 #include "source_base/global_variable.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_base/timer.h"
 #include "source_base/tool_quit.h"
 #include "source_hsolver/module_genelpa/elpa_new.h"
@@ -59,7 +58,7 @@ void DiagoElpaNative<T>::diag_pool(hamilt::MatrixBlock<T>& h_mat,
 
     ModuleBase::timer::start("DiagoElpaNative", "elpa_solve");
 
-    int nev = PARAM.inp.nbands;
+    int nev = this->nbands;
     int narows = h_mat.row;
     int nacols = h_mat.col;
 
@@ -70,7 +69,7 @@ void DiagoElpaNative<T>::diag_pool(hamilt::MatrixBlock<T>& h_mat,
     int nprows, npcols, myprow, mypcol;
 
     Cblacs_gridinfo(cblacs_ctxt, &nprows, &npcols, &myprow, &mypcol);
-    std::vector<Real> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(narows * nacols);
 
     if (elpa_init(20210430) != ELPA_OK)
@@ -107,7 +106,7 @@ void DiagoElpaNative<T>::diag_pool(hamilt::MatrixBlock<T>& h_mat,
     #define ELPA_WITH_SYCL_GPU_VERSION 0
  */
 #if ELPA_WITH_NVIDIA_GPU_VERSION
-    if (PARAM.inp.device == "gpu")
+    if (this->use_gpu)
     {
         elpa_set(handle, "nvidia-gpu", 1, &success);
         elpa_set(handle, "real_kernel", ELPA_2STAGE_REAL_NVIDIA_GPU, &success);
@@ -138,7 +137,7 @@ void DiagoElpaNative<T>::diag_pool(hamilt::MatrixBlock<T>& h_mat,
     }
 
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
     const int size = psi.get_nbands() * psi.get_nbasis();
     BlasConnector::copy(size, eigenvectors.data(), inc, psi.get_pointer(), inc);
 }
diff --git a/source/source_hsolver/diago_elpa_native.h b/source/source_hsolver/diago_elpa_native.h
index 2c556254e8..2cbd2c27d6 100644
--- a/source/source_hsolver/diago_elpa_native.h
+++ b/source/source_hsolver/diago_elpa_native.h
@@ -15,6 +15,12 @@ class DiagoElpaNative
     using Real = typename GetTypeReal<T>::type;
 
   public:
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    /// @param use_gpu_in offload to the NVIDIA-GPU ELPA kernels when ELPA was built with GPU support
+    DiagoElpaNative(const int nlocal_in, const int nbands_in, const bool use_gpu_in)
+        : nlocal(nlocal_in), nbands(nbands_in), use_gpu(use_gpu_in) {};
+
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
@@ -27,6 +33,10 @@ class DiagoElpaNative
 
     static int DecomposedState;
 
+  private:
+    const int nlocal;
+    const int nbands;
+    const bool use_gpu;
 };
 
 template <typename T>
diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp
index 6bdba41ae9..0020efc1e8 100644
--- a/source/source_hsolver/diago_lapack.cpp
+++ b/source/source_hsolver/diago_lapack.cpp
@@ -1,7 +1,5 @@
 // Refactored according to diago_scalapack
 // This code will be futher refactored to remove the dependency of psi and hamilt
-#include "source_io/module_parameter/parameter.h"
-
 #include "diago_lapack.h"
 
 #include "source_base/global_variable.h"
@@ -42,14 +40,14 @@ void DiagoLapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>&
     phm_in->matrix(h_mat, s_mat);
 
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
 
     // Diag
     this->dsygvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     // Copy result
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 
 template <>
@@ -61,11 +59,11 @@ void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>
     matcd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
     this->zhegvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 
 #ifdef __MPI
@@ -78,11 +76,11 @@ void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>
 {
     ModuleBase::TITLE("DiagoLapack", "diag_pool");
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
     this->dsygvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
     template<>
     void DiagoLapack<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex<double>>& h_mat,
@@ -93,11 +91,11 @@ void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>
 {
     ModuleBase::TITLE("DiagoLapack", "diag_pool");
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
     this->zhegvx_diag(h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 #endif
 
@@ -115,20 +113,20 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::dsygvx_once(const int ncol,
     memcpy(s_tmp.c, s_mat, sizeof(double) * ncol * nrow);
 
     const char jobz = 'V', range = 'I', uplo = 'U';
-    const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
+    const int itype = 1, il = 1, iu = this->nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0;
     double vl = 0, vu = 0;
     const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC;
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
-    std::vector<int> ifail(PARAM.globalv.nlocal, 0);
+    std::vector<int> ifail(this->nlocal, 0);
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> gap(GlobalV::DSIZE);
 
     // LAPACK dsygvx signature:
     // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU,
     //  ABSTOL, M, W, Z, LDZ, WORK, LWORK, IWORK, IFAIL, INFO)
-    int n = PARAM.globalv.nlocal;
+    int n = this->nlocal;
     int lda = n, ldb = n, ldz = n;
     dsygvx_(&itype,
         &jobz,
@@ -224,7 +222,7 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     memcpy(s_tmp.c, s_mat, sizeof(std::complex<double>) * ncol * nrow);
 
     const char jobz = 'V', range = 'I', uplo = 'U';
-    const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
+    const int itype = 1, il = 1, iu = this->nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0;
     const double abstol = LAPACK_ABSTOL, orfac = LAPACK_ORFAC;
     
@@ -232,14 +230,14 @@ std::pair<int, std::vector<int>> DiagoLapack<T>::zhegvx_once(const int ncol,
     std::vector<std::complex<double>> work(1, 0);
     std::vector<double> rwork(3, 0);
     std::vector<int> iwork(1, 0);
-    std::vector<int> ifail(PARAM.globalv.nlocal, 0);
+    std::vector<int> ifail(this->nlocal, 0);
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> gap(GlobalV::DSIZE);
 
     // LAPACK zhegvx signature:
     // (ITYPE, JOBZ, RANGE, UPLO, N, A, LDA, B, LDB, VL, VU, IL, IU,
     //  ABSTOL, M, W, Z, LDZ, WORK, LWORK, RWORK, IWORK, IFAIL, INFO)
-    int n = PARAM.globalv.nlocal;
+    int n = this->nlocal;
     int lda = n, ldb = n, ldz = n;
     zhegvx_(&itype,
         &jobz,
@@ -416,8 +414,7 @@ void DiagoLapack<T>::post_processing(const int info, const std::vector<int>& vec
     {
         const std::string str_M = "M = " + ModuleBase::GlobalFunc::TO_STRING(vec[0]) + ".\n";
         const std::string str_NZ = "NZ = " + ModuleBase::GlobalFunc::TO_STRING(vec[1]) + ".\n";
-        const std::string str_NBANDS
-            = "PARAM.inp.nbands = " + ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.nbands) + ".\n";
+        const std::string str_NBANDS = "nbands = " + ModuleBase::GlobalFunc::TO_STRING(this->nbands) + ".\n";
         throw std::runtime_error(str_info_FILE + str_M + str_NZ + str_NBANDS);
     }
     else if (info / 16 % 2)
diff --git a/source/source_hsolver/diago_lapack.h b/source/source_hsolver/diago_lapack.h
index bfdf78ac34..9f674fc20f 100644
--- a/source/source_hsolver/diago_lapack.h
+++ b/source/source_hsolver/diago_lapack.h
@@ -26,6 +26,10 @@ class DiagoLapack
     using Real = typename GetTypeReal<T>::type;
 
   public:
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    DiagoLapack(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {};
+
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
   #ifdef __MPI
     // diagnolization used in parallel-k case
@@ -61,6 +65,10 @@ class DiagoLapack
     int degeneracy_max = 12; // For reorthogonalized memory. 12 followes siesta.
 
     void post_processing(const int info, const std::vector<int>& vec);
+
+  private:
+    const int nlocal;
+    const int nbands;
 };
 
 } // namespace hsolver
diff --git a/source/source_hsolver/diago_scalapack.cpp b/source/source_hsolver/diago_scalapack.cpp
index 366478fdc2..f2ddb8289b 100644
--- a/source/source_hsolver/diago_scalapack.cpp
+++ b/source/source_hsolver/diago_scalapack.cpp
@@ -1,6 +1,5 @@
 //=====================
 // AUTHOR : Peize Lin
-#include "source_io/module_parameter/parameter.h"
 // DATE : 2021-11-02
 // REFACTORING AUTHOR : Daye Zheng
 // DATE : 2022-04-14
@@ -28,10 +27,10 @@ namespace hsolver
     matd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     this->pdsygvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
     template<>
     void DiagoScalapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>>* phm_in, psi::Psi<std::complex<double>>& psi, Real* eigenvalue_in)
@@ -40,10 +39,10 @@ namespace hsolver
     matcd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     this->pzhegvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 
 #ifdef __MPI
@@ -56,10 +55,10 @@ namespace hsolver
 {
     ModuleBase::TITLE("DiagoScalapack", "diag_pool");
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     this->pdsygvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
     template<>
     void DiagoScalapack<std::complex<double>>::diag_pool(hamilt::MatrixBlock<std::complex<double>>& h_mat,
@@ -70,10 +69,10 @@ namespace hsolver
 {
     ModuleBase::TITLE("DiagoScalapack", "diag_pool");
     assert(h_mat.col == s_mat.col && h_mat.row == s_mat.row && h_mat.desc == s_mat.desc);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
+    std::vector<double> eigen(this->nlocal, 0.0);
     this->pzhegvx_diag(h_mat.desc, h_mat.col, h_mat.row, h_mat.p, s_mat.p, eigen.data(), psi);
     const int inc = 1;
-    BlasConnector::copy(PARAM.inp.nbands, eigen.data(), inc, eigenvalue_in, inc);
+    BlasConnector::copy(this->nbands, eigen.data(), inc, eigenvalue_in, inc);
 }
 #endif
 
@@ -92,13 +91,13 @@ namespace hsolver
     memcpy(s_tmp.c, s_mat, sizeof(double) * ncol * nrow);
 
     const char jobz = 'V', range = 'I', uplo = 'U';
-    const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
+    const int itype = 1, il = 1, iu = this->nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, liwork = -1, info = 0;
     double vl = 0, vu = 0;
     const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC;
     std::vector<double> work(3, 0);
     std::vector<int> iwork(1, 0);
-    std::vector<int> ifail(PARAM.globalv.nlocal, 0);
+    std::vector<int> ifail(this->nlocal, 0);
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> gap(GlobalV::DSIZE);
 
@@ -106,7 +105,7 @@ namespace hsolver
              &jobz,
              &range,
              &uplo,
-             &PARAM.globalv.nlocal,
+             &this->nlocal,
              h_tmp.c,
              &one,
              &one,
@@ -152,7 +151,7 @@ namespace hsolver
              &jobz,
              &range,
              &uplo,
-             &PARAM.globalv.nlocal,
+             &this->nlocal,
              h_tmp.c,
              &one,
              &one,
@@ -217,7 +216,7 @@ namespace hsolver
     memcpy(s_tmp.c, s_mat, sizeof(std::complex<double>) * ncol * nrow);
 
     const char jobz = 'V', range = 'I', uplo = 'U';
-    const int itype = 1, il = 1, iu = PARAM.inp.nbands, one = 1;
+    const int itype = 1, il = 1, iu = this->nbands, one = 1;
     int M = 0, NZ = 0, lwork = -1, lrwork = -1, liwork = -1, info = 0;
     const double abstol = SCALAPACK_ABSTOL, orfac = SCALAPACK_ORFAC;
     //Note: pzhegvx_ has a bug
@@ -227,7 +226,7 @@ namespace hsolver
     std::vector<std::complex<double>> work(1, 0);
     std::vector<double> rwork(3, 0);
     std::vector<int> iwork(1, 0);
-    std::vector<int> ifail(PARAM.globalv.nlocal, 0);
+    std::vector<int> ifail(this->nlocal, 0);
     std::vector<int> iclustr(2 * GlobalV::DSIZE);
     std::vector<double> gap(GlobalV::DSIZE);
 
@@ -235,7 +234,7 @@ namespace hsolver
              &jobz,
              &range,
              &uplo,
-             &PARAM.globalv.nlocal,
+             &this->nlocal,
              h_tmp.c,
              &one,
              &one,
@@ -276,7 +275,7 @@ namespace hsolver
     //	GlobalV::ofs_running<<"lwork="<<work[0]<<"\t"<<"lrwork="<<rwork[0]<<"\t"<<"liwork="<<iwork[0]<<std::endl;
     lwork = work[0].real();
     work.resize(lwork, 0);
-    lrwork = rwork[0] + this->degeneracy_max * PARAM.globalv.nlocal;
+    lrwork = rwork[0] + this->degeneracy_max * this->nlocal;
     int maxlrwork = std::max(lrwork,3);
     rwork.resize(maxlrwork, 0);
     liwork = iwork[0];
@@ -286,7 +285,7 @@ namespace hsolver
              &jobz,
              &range,
              &uplo,
-             &PARAM.globalv.nlocal,
+             &this->nlocal,
              h_tmp.c,
              &one,
              &one,
@@ -431,7 +430,7 @@ namespace hsolver
         const std::string str_M = "M = " + ModuleBase::GlobalFunc::TO_STRING(vec[0]) + ".\n";
         const std::string str_NZ = "NZ = " + ModuleBase::GlobalFunc::TO_STRING(vec[1]) + ".\n";
         const std::string str_NBANDS
-            = "PARAM.inp.nbands = " + ModuleBase::GlobalFunc::TO_STRING(PARAM.inp.nbands) + ".\n";
+            = "nbands = " + ModuleBase::GlobalFunc::TO_STRING(this->nbands) + ".\n";
         throw std::runtime_error(str_info_FILE + str_M + str_NZ + str_NBANDS);
     }
     else if (info / 16 % 2)
diff --git a/source/source_hsolver/diago_scalapack.h b/source/source_hsolver/diago_scalapack.h
index 05b1144836..7dce030616 100644
--- a/source/source_hsolver/diago_scalapack.h
+++ b/source/source_hsolver/diago_scalapack.h
@@ -27,6 +27,10 @@ namespace hsolver
 private:
     using Real = typename GetTypeReal<T>::type;
   public:
+    /// @param nlocal_in global dimension of the NAO Hamiltonian
+    /// @param nbands_in number of lowest eigenpairs to compute
+    DiagoScalapack(const int nlocal_in, const int nbands_in) : nlocal(nlocal_in), nbands(nbands_in) {};
+
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
@@ -67,6 +71,9 @@ namespace hsolver
     int degeneracy_max = 12; // For reorthogonalized memory. 12 followes siesta.
 
     void post_processing(const int info, const std::vector<int> &vec);
+
+    const int nlocal;
+    const int nbands;
 };
 
 } // namespace hsolver
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
index d31658d8ad..3c8bcc3ff0 100644
--- a/source/source_hsolver/hsolver_lcao.cpp
+++ b/source/source_hsolver/hsolver_lcao.cpp
@@ -141,19 +141,19 @@ void HSolverLCAO<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>&
     if (this->method == "scalapack_gvx")
     {
 #ifdef __MPI
-        DiagoScalapack<T> sa;
+        DiagoScalapack<T> sa(this->nlocal, this->nbands);
         sa.diag(hm, psi, eigenvalue);
 #endif
     }
 #ifdef __ELPA
     else if (this->method == "genelpa")
     {
-        DiagoElpa<T> el;
+        DiagoElpa<T> el(this->nlocal, this->nbands);
         el.diag(hm, psi, eigenvalue);
     }
     else if (this->method == "elpa")
     {
-        DiagoElpaNative<T> el;
+        DiagoElpaNative<T> el(this->nlocal, this->nbands, this->use_gpu);
         el.diag(hm, psi, eigenvalue);
     }
 #endif
@@ -161,7 +161,7 @@ void HSolverLCAO<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>&
     else if (this->method == "cusolver")
     {
         // Note: This branch will only be executed in the single-process case
-        DiagoCusolver<T> cu;
+        DiagoCusolver<T> cu(this->nlocal, this->nbands);
         hamilt::MatrixBlock<T> hk, sk;
         hm->matrix(hk, sk);
         cu.diag(hk, sk, psi, eigenvalue);
@@ -169,14 +169,14 @@ void HSolverLCAO<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>&
 #ifdef __CUSOLVERMP
     else if (this->method == "cusolvermp")
     {
-        DiagoCusolverMP<T> cm;
+        DiagoCusolverMP<T> cm(this->nlocal, this->nbands);
         cm.diag(hm, psi, eigenvalue);
     }
 #endif
 #endif
     else if (this->method == "lapack") // only for single core
     {
-        DiagoLapack<T> la;
+        DiagoLapack<T> la(this->nlocal, this->nbands);
         la.diag(hm, psi, eigenvalue);
     }
     else
@@ -253,23 +253,23 @@ void HSolverLCAO<T, Device>::parakSolve(hamilt::Hamilt<T>* pHamilt,
             /// solve eigenvector and eigenvalue for H(k)
             if (this->method == "scalapack_gvx")
             {
-                DiagoScalapack<T> sa;
+                DiagoScalapack<T> sa(this->nlocal, this->nbands);
                 sa.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D);
             }
             else if (this->method == "lapack")
             {
-                DiagoLapack<T> la;
+                DiagoLapack<T> la(this->nlocal, this->nbands);
                 la.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D);
             }
 #ifdef __ELPA
             else if (this->method == "genelpa")
             {
-                DiagoElpa<T> el;
+                DiagoElpa<T> el(this->nlocal, this->nbands);
                 el.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D);
             }
             else if (this->method == "elpa")
             {
-                DiagoElpaNative<T> el;
+                DiagoElpaNative<T> el(this->nlocal, this->nbands, this->use_gpu);
                 el.diag_pool(hk_pool, sk_pool, psi_pool, &(pes->ekb(ik_global, 0)), k2d.POOL_WORLD_K2D);
             }
 #endif
@@ -427,7 +427,7 @@ void HSolverLCAO<T, Device>::parakSolve_cusolver(hamilt::Hamilt<T>* pHamilt,
         if(kpt_assigned != -1)
         {
             psi_local.resize(1, ncol, nrow);
-            DiagoCusolver<T> cu{};
+            DiagoCusolver<T> cu(this->nlocal, this->nbands);
             hamilt::MatrixBlock<T> hk_local = hamilt::MatrixBlock<T>{
                     hk_mat.data(), (size_t)nrow, (size_t)ncol,
                     mat_para_local.desc};
diff --git a/source/source_hsolver/hsolver_lcao.h b/source/source_hsolver/hsolver_lcao.h
index 60a777e3a4..3b54d92d1f 100644
--- a/source/source_hsolver/hsolver_lcao.h
+++ b/source/source_hsolver/hsolver_lcao.h
@@ -19,8 +19,11 @@ class HSolverLCAO
                 const std::string method_in,
                 const int kpar_lcao_in,
                 const int nlocal_in,
-                const double nelec_in)
-        : ParaV(ParaV_in), method(method_in), kpar_lcao(kpar_lcao_in), nlocal(nlocal_in), nelec(nelec_in) {};
+                const int nbands_in,
+                const double nelec_in,
+                const bool use_gpu_in)
+        : ParaV(ParaV_in), method(method_in), kpar_lcao(kpar_lcao_in), nlocal(nlocal_in), nbands(nbands_in),
+          nelec(nelec_in), use_gpu(use_gpu_in) {};
 
     void solve(hamilt::Hamilt<TK>* pHamilt,
                psi::Psi<TK>& psi,
@@ -49,8 +52,10 @@ class HSolverLCAO
     const std::string method;
 
     const int kpar_lcao; // number of pools for LCAO diagonalization
-    const int nlocal;    // global dimension of the NAO Hamiltonian, only used by the pexsi branch
+    const int nlocal;    // global dimension of the NAO Hamiltonian
+    const int nbands;    // number of bands to be solved for
     const double nelec;  // total number of electrons, only used by the pexsi branch
+    const bool use_gpu;  // true if running on GPU, only used by the native-ELPA branch
 };
 
 } // namespace hsolver
diff --git a/source/source_hsolver/test/diago_lapack_test.cpp b/source/source_hsolver/test/diago_lapack_test.cpp
index 8620eaafa4..e3d358052e 100644
--- a/source/source_hsolver/test/diago_lapack_test.cpp
+++ b/source/source_hsolver/test/diago_lapack_test.cpp
@@ -1,8 +1,5 @@
 // Author: Zhang Xiaoyang
 // A modified version of diago_lcao_test.cpp
-#define private public
-#include "source_io/module_parameter/parameter.h"
-#undef private
 // Remove some useless functions and dependencies. Serialized the full code
 // and refactored some function.
 
@@ -184,21 +181,14 @@ class DiagoLapackPrepare
         hmtest.ncol = nlocal;
     }
 
-    void set_env()
-    {
-        PARAM.sys.nlocal = nlocal;
-        PARAM.input.nbands = nbands;
-    }
-
     void diago()
     {
         this->pb2d();
         this->print_hs();
-        this->set_env();
 
         for (int i = 0; i < REPEATRUN; i++)
         {
-            hsolver::DiagoLapack<T> dh;
+            hsolver::DiagoLapack<T> dh(nlocal, nbands);
             dh.diag(&hmtest, psi, e_solver.data());
             // dh->diag(&hmtest, psi, e_solver.data());
         }
diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
index 9f340419b7..3a5154af10 100644
--- a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
@@ -1,8 +1,5 @@
 #include "source_hsolver/diago_scalapack.h"
 #include "source_hsolver/test/diago_elpa_utils.h"
-#define private public
-#include "source_io/module_parameter/parameter.h"
-#undef private
 #include "mpi.h"
 #include "string.h"
 
@@ -204,8 +201,6 @@ class DiagoPrepare
 
     void set_env()
     {
-        PARAM.sys.nlocal = nlocal;
-        PARAM.input.nbands = nbands;
         GlobalV::DSIZE = dsize;
     }
 
@@ -225,13 +220,13 @@ class DiagoPrepare
             hmtest.s_local = this->s_local;
             if (ks_solver == "scalapack_gvx")
             {
-                hsolver::DiagoScalapack<T> dh;
+                hsolver::DiagoScalapack<T> dh(nlocal, nbands);
                 dh.diag(&hmtest, psi, e_solver.data());
             }
     #ifdef __CUDA
             else if (ks_solver == "cusolver")
                 {
-                    hsolver::DiagoCusolver<T> dh;
+                    hsolver::DiagoCusolver<T> dh(nlocal, nbands);
                     hamilt::MatrixBlock<T> h_mat, s_mat;
                     hmtest.matrix(h_mat, s_mat);
                     dh.diag(h_mat, s_mat, psi, e_solver.data());
diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
index d0ed0fdae6..6d66269415 100644
--- a/source/source_hsolver/test/diago_lcao_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -1,9 +1,6 @@
 #include "source_hsolver/diago_scalapack.h"
 #include "source_hsolver/diago_lapack.h"
 #include "source_hsolver/test/diago_elpa_utils.h"
-#define private public
-#include "source_io/module_parameter/parameter.h"
-#undef private
 #include "mpi.h"
 #include "string.h"
 
@@ -205,8 +202,6 @@ class DiagoPrepare
 
     void set_env()
     {
-        PARAM.sys.nlocal = nlocal;
-        PARAM.input.nbands = nbands;
         GlobalV::DSIZE = dsize;
     }
 
@@ -226,18 +221,18 @@ class DiagoPrepare
             hmtest.s_local = this->s_local;
             if (ks_solver == "scalapack_gvx")
             {
-                hsolver::DiagoScalapack<T> dh;
+                hsolver::DiagoScalapack<T> dh(nlocal, nbands);
                 dh.diag(&hmtest, psi, e_solver.data());
             }
             else if (ks_solver == "lapack")
             {
-                hsolver::DiagoLapack<T> la;
+                hsolver::DiagoLapack<T> la(nlocal, nbands);
                 la.diag(&hmtest, psi, e_solver.data());
             }
     #ifdef __ELPA
             else if (ks_solver == "genelpa")
             {
-                hsolver::DiagoElpa<T> dh;
+                hsolver::DiagoElpa<T> dh(nlocal, nbands);
                 dh.diag(&hmtest, psi, e_solver.data());
             }
     #endif
diff --git a/source/source_lcao/LCAO_set.cpp b/source/source_lcao/LCAO_set.cpp
index c27b4485f6..ad5345585e 100644
--- a/source/source_lcao/LCAO_set.cpp
+++ b/source/source_lcao/LCAO_set.cpp
@@ -234,7 +234,9 @@ void LCAO_domain::init_chg_hr(
                                               ks_solver,
                                               PARAM.globalv.kpar_lcao,
                                               PARAM.globalv.nlocal,
-                                              PARAM.inp.nelec);
+                                              PARAM.inp.nbands,
+                                              PARAM.inp.nelec,
+                                              PARAM.inp.device == "gpu");
     hsolver_lcao_obj.solve(p_hamilt, psi, pelec, dm, chr, nspin, 0);
 }
 
diff --git a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
index 65cbd88608..fc2fdc06b3 100644
--- a/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
+++ b/source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp
@@ -525,7 +525,9 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(
                                                             PARAM.inp.ks_solver,
                                                             PARAM.globalv.kpar_lcao,
                                                             PARAM.globalv.nlocal,
-                                                            PARAM.inp.nelec);
+                                                            PARAM.inp.nbands,
+                                                            PARAM.inp.nelec,
+                                                            PARAM.inp.device == "gpu");
         if (this->nspin_ == 2)
         {
             dynamic_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>*>(this->p_operator)