diff --git a/python/pyabacus/src/py_diago_dav_subspace.hpp b/python/pyabacus/src/py_diago_dav_subspace.hpp
index bd8bbb3e41..9990acc07b 100644
--- a/python/pyabacus/src/py_diago_dav_subspace.hpp
+++ b/python/pyabacus/src/py_diago_dav_subspace.hpp
@@ -113,23 +113,21 @@ class PyDiagoDavSubspace
         auto hpsi_func = [mm_op] (
             std::complex<double> *psi_in,
             std::complex<double> *hpsi_out, 
-            const int nband_in,
-            const int nbasis_in, 
-            const int band_index1,
-            const int band_index2
+            const int ld_psi,
+            const int nvec
         ) {
             // Note: numpy's py::array_t is row-major, but
             //       our raw pointer-array is column-major
-            py::array_t<std::complex<double>, py::array::f_style> psi({nbasis_in, band_index2 - band_index1 + 1});
+            py::array_t<std::complex<double>, py::array::f_style> psi({ld_psi, nvec});
             py::buffer_info psi_buf = psi.request();
             std::complex<double>* psi_ptr = static_cast<std::complex<double>*>(psi_buf.ptr);
-            std::copy(psi_in + band_index1 * nbasis_in, psi_in + (band_index2 + 1) * nbasis_in, psi_ptr);
+            std::copy(psi_in, psi_in + nvec * ld_psi, psi_ptr);
 
             py::array_t<std::complex<double>, py::array::f_style> hpsi = mm_op(psi);
 
             py::buffer_info hpsi_buf = hpsi.request();
             std::complex<double>* hpsi_ptr = static_cast<std::complex<double>*>(hpsi_buf.ptr);
-            std::copy(hpsi_ptr, hpsi_ptr + (band_index2 - band_index1 + 1) * nbasis_in, hpsi_out);
+            std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
         };
 
         obj = std::make_unique<hsolver::Diago_DavSubspace<std::complex<double>, base_device::DEVICE_CPU>>(
diff --git a/python/pyabacus/src/py_diago_david.hpp b/python/pyabacus/src/py_diago_david.hpp
index 2db5284f4b..cb57451d91 100644
--- a/python/pyabacus/src/py_diago_david.hpp
+++ b/python/pyabacus/src/py_diago_david.hpp
@@ -111,23 +111,21 @@ class PyDiagoDavid
         auto hpsi_func = [mm_op] (
             std::complex<double> *psi_in,
             std::complex<double> *hpsi_out, 
-            const int nband_in, 
-            const int nbasis_in, 
-            const int band_index1, 
-            const int band_index2
+            const int ld_psi, 
+            const int nvec
         ) {
             // Note: numpy's py::array_t is row-major, but
             //       our raw pointer-array is column-major
-            py::array_t<std::complex<double>, py::array::f_style> psi({nbasis_in, band_index2 - band_index1 + 1});
+            py::array_t<std::complex<double>, py::array::f_style> psi({ld_psi, nvec});
             py::buffer_info psi_buf = psi.request();
             std::complex<double>* psi_ptr = static_cast<std::complex<double>*>(psi_buf.ptr);
-            std::copy(psi_in + band_index1 * nbasis_in, psi_in + (band_index2 + 1) * nbasis_in, psi_ptr);
+            std::copy(psi_in, psi_in + nvec * ld_psi, psi_ptr);
 
             py::array_t<std::complex<double>, py::array::f_style> hpsi = mm_op(psi);
 
             py::buffer_info hpsi_buf = hpsi.request();
             std::complex<double>* hpsi_ptr = static_cast<std::complex<double>*>(hpsi_buf.ptr);
-            std::copy(hpsi_ptr, hpsi_ptr + (band_index2 - band_index1 + 1) * nbasis_in, hpsi_out);
+            std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
         };
 
         auto spsi_func = [this] (
diff --git a/source/module_hsolver/diago_dav_subspace.cpp b/source/module_hsolver/diago_dav_subspace.cpp
index bcb9780271..6256d6597c 100644
--- a/source/module_hsolver/diago_dav_subspace.cpp
+++ b/source/module_hsolver/diago_dav_subspace.cpp
@@ -124,7 +124,8 @@ int Diago_DavSubspace<T, Device>::diag_once(const HPsiFunc& hpsi_func,
 
     // compute h*psi_in_iter
     // NOTE: bands after the first n_band should yield zero
-    hpsi_func(this->psi_in_iter, this->hphi, this->nbase_x, this->dim, 0, this->nbase_x - 1);
+    // hphi[:, 0:nbase_x] = H * psi_in_iter[:, 0:nbase_x]
+    hpsi_func(this->psi_in_iter, this->hphi, this->dim, this->nbase_x);
 
     // at this stage, notconv = n_band and nbase = 0
     // note that nbase of cal_elem is an inout parameter: nbase := nbase + notconv
@@ -421,7 +422,8 @@ void Diago_DavSubspace<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
     }
 
     // update hpsi[:, nbase:nbase+notconv]
-    hpsi_func(psi_iter, &hphi[nbase * this->dim], this->nbase_x, this->dim, nbase, nbase + notconv - 1);
+    // hpsi[:, nbase:nbase+notconv] = H * psi_iter[:, nbase:nbase+notconv]
+    hpsi_func(psi_iter + nbase * dim, hphi + nbase * this->dim, this->dim, notconv);
 
     ModuleBase::timer::tick("Diago_DavSubspace", "cal_grad");
     return;
@@ -886,7 +888,8 @@ void Diago_DavSubspace<T, Device>::diagH_subspace(T* psi_pointer, // [in] & [out
 
     {
         // do hPsi for all bands
-        hpsi_func(psi_pointer, hphi, n_band, dmax, 0, nstart - 1);
+        // hphi[:, 0:nstart] = H * psi_pointer[:, 0:nstart]
+        hpsi_func(psi_pointer, hphi, dmax, nstart);
 
         gemm_op<T, Device>()(ctx,
                              'C',
diff --git a/source/module_hsolver/diago_dav_subspace.h b/source/module_hsolver/diago_dav_subspace.h
index 78553f0174..8e26a08b45 100644
--- a/source/module_hsolver/diago_dav_subspace.h
+++ b/source/module_hsolver/diago_dav_subspace.h
@@ -31,7 +31,8 @@ class Diago_DavSubspace : public DiagH<T, Device>
 
     virtual ~Diago_DavSubspace() override;
 
-    using HPsiFunc = std::function<void(T*, T*, const int, const int, const int, const int)>;
+    // See diago_david.h for information on the HPsiFunc function type
+    using HPsiFunc = std::function<void(T*, T*, const int, const int)>;
 
     int diag(const HPsiFunc& hpsi_func,
              T* psi_in,
diff --git a/source/module_hsolver/diago_david.cpp b/source/module_hsolver/diago_david.cpp
index 8d0a2a8235..0590a810cd 100644
--- a/source/module_hsolver/diago_david.cpp
+++ b/source/module_hsolver/diago_david.cpp
@@ -152,7 +152,7 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
                                      const SPsiFunc& spsi_func,
                                      const int dim,
                                      const int nband,
-                                     const int ldPsi,
+                                     const int ld_psi,
                                      T *psi_in,
                                      Real* eigenvalue_in,
                                      const Real david_diag_thr,
@@ -191,20 +191,20 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
         if(this->use_paw)
         {
 #ifdef USE_PAW
-            GlobalC::paw_cell.paw_nl_psi(1, reinterpret_cast<const std::complex<double>*> (psi_in + m*ldPsi),
+            GlobalC::paw_cell.paw_nl_psi(1, reinterpret_cast<const std::complex<double>*> (psi_in + m*ld_psi),
                 reinterpret_cast<std::complex<double>*>(&this->spsi[m * dim]));
 #endif
         }
         else
         {
-            // phm_in->sPsi(psi_in + m*ldPsi, &this->spsi[m * dim], dim, dim, 1);
-            spsi_func(psi_in + m*ldPsi,&this->spsi[m*dim],dim,dim,1);
+            // phm_in->sPsi(psi_in + m*ld_psi, &this->spsi[m * dim], dim, dim, 1);
+            spsi_func(psi_in + m*ld_psi,&this->spsi[m*dim],dim,dim,1);
         }
     }
     // begin SchmidtOrth
     for (int m = 0; m < nband; m++)
     {
-        syncmem_complex_op()(this->ctx, this->ctx, basis + dim*m, psi_in + m*ldPsi, dim);
+        syncmem_complex_op()(this->ctx, this->ctx, basis + dim*m, psi_in + m*ld_psi, dim);
 
         this->SchmidtOrth(dim,
                          nband,
@@ -230,7 +230,9 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
     // end of SchmidtOrth and calculate H|psi>
     // hpsi_info dav_hpsi_in(&basis, psi::Range(true, 0, 0, nband - 1), this->hpsi);
     // phm_in->ops->hPsi(dav_hpsi_in);
-    hpsi_func(basis, hpsi, nbase_x, dim, 0, nband - 1);
+    // hpsi[:, 0:nband] = H basis[:, 0:nband]
+    // slice index in this piece of code is in C manner. i.e. 0:id stands for [0,id)
+    hpsi_func(basis, hpsi, dim, nband);
 
     this->cal_elem(dim, nbase, nbase_x, this->notconv, this->hpsi, this->spsi, this->hcc, this->scc);
 
@@ -287,7 +289,7 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
 
             // update eigenvectors of Hamiltonian
 
-            setmem_complex_op()(this->ctx, psi_in, 0, nband * ldPsi);
+            setmem_complex_op()(this->ctx, psi_in, 0, nband * ld_psi);
             //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
             gemm_op<T, Device>()(this->ctx,
                                       'N',
@@ -302,7 +304,7 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
                                       nbase_x,
                                       this->zero,
                                       psi_in,       // C dim * nband
-                                      ldPsi
+                                      ld_psi
             );
 
             if (!this->notconv || (dav_iter == david_maxiter))
@@ -322,7 +324,7 @@ int DiagoDavid<T, Device>::diag_once(const HPsiFunc& hpsi_func,
                               nbase_x,
                               eigenvalue_in,
                               psi_in,
-                              ldPsi,
+                              ld_psi,
                               this->hpsi,
                               this->spsi,
                               this->hcc,
@@ -601,7 +603,8 @@ void DiagoDavid<T, Device>::cal_grad(const HPsiFunc& hpsi_func,
     //                       psi::Range(true, 0, nbase, nbase + notconv - 1),
     //                       &hpsi[nbase * dim]); // &hp(nbase, 0)
     // phm_in->ops->hPsi(dav_hpsi_in);
-    hpsi_func(basis, &hpsi[nbase * dim], nbase_x, dim, nbase, nbase + notconv - 1);
+    // hpsi[:, nbase:nbase+notcnv] = H basis[:, nbase:nbase+notcnv]
+    hpsi_func(basis + nbase * dim, hpsi + nbase * dim, dim, notconv);
 
     delmem_complex_op()(this->ctx, lagrange);
     delmem_complex_op()(this->ctx, vc_ev_vector);
@@ -785,7 +788,7 @@ void DiagoDavid<T, Device>::diag_zhegvx(const int& nbase,
  * @param nbase_x The maximum dimension of the reduced basis set.
  * @param eigenvalue_in Pointer to the array of eigenvalues.
  * @param psi_in Pointer to the array of wavefunctions.
- * @param ldPsi The leading dimension of the wavefunction array.
+ * @param ld_psi The leading dimension of the wavefunction array.
  * @param hpsi Pointer to the output array for the updated basis set.
  * @param spsi Pointer to the output array for the updated basis set (nband-th column).
  * @param hcc Pointer to the output array for the updated reduced Hamiltonian.
@@ -800,7 +803,7 @@ void DiagoDavid<T, Device>::refresh(const int& dim,
                                          const int nbase_x, // maximum dimension of the reduced basis set
                                          const Real* eigenvalue_in,
                                          const T *psi_in,
-                                         const int ldPsi,
+                                         const int ld_psi,
                                          T* hpsi,
                                          T* spsi,
                                          T* hcc,
@@ -866,7 +869,7 @@ void DiagoDavid<T, Device>::refresh(const int& dim,
 
     for (int m = 0; m < nband; m++)
     {
-        syncmem_complex_op()(this->ctx, this->ctx, basis + dim*m,psi_in + m*ldPsi, dim);
+        syncmem_complex_op()(this->ctx, this->ctx, basis + dim*m,psi_in + m*ld_psi, dim);
         /*for (int ig = 0; ig < npw; ig++)
             basis(m, ig) = psi(m, ig);*/
     }
@@ -1149,15 +1152,13 @@ void DiagoDavid<T, Device>::planSchmidtOrth(const int nband, std::vector<int>& p
 /**
  * @brief Performs iterative diagonalization using the David algorithm.
  * 
- * @warning Please see docs of `HPsiFunc` for more information.
- * @warning Please adhere strictly to the requirements of the function pointer
- * @warning for the hpsi mat-vec interface; it may seem counterintuitive.
+ * @warning Please see docs of `HPsiFunc` for more information about the hpsi mat-vec interface.
  * 
  * @tparam T The type of the elements in the matrix.
  * @tparam Device The device type (CPU or GPU).
  * @param hpsi_func The function object that computes the matrix-blockvector product H * psi.
  * @param spsi_func The function object that computes the matrix-blockvector product overlap S * psi.
- * @param ldPsi The leading dimension of the psi_in array.
+ * @param ld_psi The leading dimension of the psi_in array.
  * @param psi_in The input wavefunction.
  * @param eigenvalue_in The array to store the eigenvalues.
  * @param david_diag_thr The convergence threshold for the diagonalization.
@@ -1172,7 +1173,7 @@ void DiagoDavid<T, Device>::planSchmidtOrth(const int nband, std::vector<int>& p
 template <typename T, typename Device>
 int DiagoDavid<T, Device>::diag(const HPsiFunc& hpsi_func,
                                 const SPsiFunc& spsi_func,
-                                const int ldPsi,
+                                const int ld_psi,
                                 T *psi_in,
                                 Real* eigenvalue_in,
                                 const Real david_diag_thr,
@@ -1187,7 +1188,7 @@ int DiagoDavid<T, Device>::diag(const HPsiFunc& hpsi_func,
     int sum_dav_iter = 0;
     do
     {
-        sum_dav_iter += this->diag_once(hpsi_func, spsi_func, dim, nband, ldPsi, psi_in, eigenvalue_in, david_diag_thr, david_maxiter);
+        sum_dav_iter += this->diag_once(hpsi_func, spsi_func, dim, nband, ld_psi, psi_in, eigenvalue_in, david_diag_thr, david_maxiter);
         ++ntry;
     } while (!check_block_conv(ntry, this->notconv, ntry_max, notconv_max));
 
diff --git a/source/module_hsolver/diago_david.h b/source/module_hsolver/diago_david.h
index a032583779..6fed3fffeb 100644
--- a/source/module_hsolver/diago_david.h
+++ b/source/module_hsolver/diago_david.h
@@ -38,50 +38,48 @@ class DiagoDavid : public DiagH<T, Device>
      * this function computes the product of the Hamiltonian matrix H and a blockvector X.
      * 
      * Called as follows:
-     * hpsi(X, HX, nvec, dim, id_start, id_end)
-     * Result is stored in HX.
-     * HX = H * X[id_start:id_end]
+     * hpsi(X, HX, ld, nvec) where X and HX are (ld, nvec)-shaped blockvectors.
+     * Result HX = H * X is stored in HX.
      *
      * @param[out] X      Head address of input blockvector of type `T*`.
-     * @param[in]  HX     Where to write output blockvector of type `T*`.
-     * @param[in]  nvec   Number of eigebpairs, i.e. number of vectors in a block.
-     * @param[in]  dim    Dimension of matrix.
-     * @param[in]  id_start Start index of blockvector.
-     * @param[in]  id_end   End index of blockvector.
+     * @param[in]  HX     Head address of output blockvector of type `T*`.
+     * @param[in]  ld     Leading dimension of blockvector.
+     * @param[in]  nvec   Number of vectors in a block.
      * 
-     * @warning HX is the exact address to store output H*X[id_start:id_end];
-     * @warning while X is the head address of input blockvector, \b without offset.
-     * @warning Calling function should pass X and HX[offset] as arguments,
-     * @warning where offset is usually id_start * leading dimension.
+     * @warning X and HX are the exact address to read input X and store output H*X,
+     * @warning both of size ld * nvec.
      */
-    using HPsiFunc = std::function<void(T*, T*, const int, const int, const int, const int)>;
+    using HPsiFunc = std::function<void(T*, T*, const int, const int)>;
 
     /**
      * @brief A function type representing the SX function.
+     * 
+     * nrow is leading dimension of spsi, npw is leading dimension of psi, nbands is number of vecs
      *
      * This function type is used to define a matrix-blockvector operator S.
      * For generalized eigenvalue problem HX = λSX,
      * this function computes the product of the overlap matrix S and a blockvector X.
      *
-     * @param[in]   X     Pointer to the input array.
-     * @param[out] SX     Pointer to the output array.
-     * @param[in] nrow    Dimension of SX: nbands * nrow.
-     * @param[in] npw     Number of plane waves.
-     * @param[in] nbands  Number of bands.
+     * @param[in]   X       Pointer to the input blockvector.
+     * @param[out] SX       Pointer to the output blockvector.
+     * @param[in] ld_spsi   Leading dimension of spsi. Dimension of SX: nbands * nrow.
+     * @param[in] ld_psi    Leading dimension of psi. Number of plane waves.
+     * @param[in] nbands    Number of vectors.
      * 
-     * @note called as spsi(in, out, dim, dim, 1)
+     * @note called like spsi(in, out, dim, dim, 1)
      */
     using SPsiFunc = std::function<void(T*, T*, const int, const int, const int)>;
 
-    int diag(const HPsiFunc& hpsi_func,           // function void hpsi(T*, T*, const int, const int, const int, const int) 
-             const SPsiFunc& spsi_func,           // function void spsi(T*, T*, const int, const int, const int) 
-                      const int ldPsi,            // Leading dimension of the psi input
-                      T *psi_in,                  // Pointer to eigenvectors
-                      Real* eigenvalue_in,        // Pointer to store the resulting eigenvalues
-                      const Real david_diag_thr,  // Convergence threshold for the Davidson iteration
-                      const int david_maxiter,    // Maximum allowed iterations for the Davidson method
-                      const int ntry_max = 5,     // Maximum number of diagonalization attempts (default is 5)
-                      const int notconv_max = 0); // Maximum number of allowed non-converged eigenvectors
+    int diag(
+      const HPsiFunc& hpsi_func,  // function void hpsi(T*, T*, const int, const int) 
+      const SPsiFunc& spsi_func,  // function void spsi(T*, T*, const int, const int, const int) 
+      const int ld_psi,           // Leading dimension of the psi input
+      T *psi_in,                  // Pointer to eigenvectors
+      Real* eigenvalue_in,        // Pointer to store the resulting eigenvalues
+      const Real david_diag_thr,  // Convergence threshold for the Davidson iteration
+      const int david_maxiter,    // Maximum allowed iterations for the Davidson method
+      const int ntry_max = 5,     // Maximum number of diagonalization attempts (5 by default)
+      const int notconv_max = 0); // Maximum number of allowed non-converged eigenvectors
 
   private:
     bool use_paw = false;
@@ -130,7 +128,7 @@ class DiagoDavid : public DiagH<T, Device>
                   const SPsiFunc& spsi_func,
                   const int dim,
                   const int nband,
-                  const int ldPsi,
+                  const int ld_psi,
                   T *psi_in,
                   Real* eigenvalue_in,
                   const Real david_diag_thr,
@@ -163,7 +161,7 @@ class DiagoDavid : public DiagH<T, Device>
                  const int nbase_x,
                  const Real* eigenvalue,
                  const T *psi_in,
-                 const int ldPsi,
+                 const int ld_psi,
                  T* hpsi,
                  T* spsi,
                  T* hcc,
@@ -171,12 +169,12 @@ class DiagoDavid : public DiagH<T, Device>
                  T* vcc);
 
     void SchmidtOrth(const int& dim,
-                    const int nband,
-                    const int m,
-                    const T* spsi,
-                    T* lagrange_m,
-                    const int mm_size,
-                    const int mv_size);
+                     const int nband,
+                     const int m,
+                     const T* spsi,
+                     T* lagrange_m,
+                     const int mm_size,
+                     const int mv_size);
 
     void planSchmidtOrth(const int nband, std::vector<int>& pre_matrix_mm_m, std::vector<int>& pre_matrix_mv_m);
 
diff --git a/source/module_hsolver/hsolver_pw.cpp b/source/module_hsolver/hsolver_pw.cpp
index 7b8e5fc0a6..678ec0adde 100644
--- a/source/module_hsolver/hsolver_pw.cpp
+++ b/source/module_hsolver/hsolver_pw.cpp
@@ -434,18 +434,17 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
     else if (this->method == "dav_subspace")
     {
         auto ngk_pointer = psi.get_ngk_pointer();
-        auto hpsi_func = [hm, ngk_pointer](T* psi_in,
-                                           T* hpsi_out,
-                                           const int nband_in,
-                                           const int nbasis_in,
-                                           const int band_index1,
-                                           const int band_index2) {
+        // hpsi_func (X, HX, ld, nvec) -> HX = H(X), X and HX blockvectors of size ld x nvec
+        auto hpsi_func = [hm, ngk_pointer](T *psi_in,
+                                           T *hpsi_out,
+                                           const int ld_psi,
+                                           const int nvec) {
             ModuleBase::timer::tick("DavSubspace", "hpsi_func");
 
             // Convert "pointer data stucture" to a psi::Psi object
-            auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nband_in, nbasis_in, ngk_pointer);
+            auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nvec, ld_psi, ngk_pointer);
 
-            psi::Range bands_range(true, 0, band_index1, band_index2);
+            psi::Range bands_range(true, 0, 0, nvec-1);
 
             using hpsi_info = typename hamilt::Operator<T, Device>::hpsi_info;
             hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
@@ -486,25 +485,23 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
         // dimensions of matrix to be solved
         const int dim = psi.get_current_nbas(); /// dimension of matrix
         const int nband = psi.get_nbands();     /// number of eigenpairs sought
-        const int ldPsi = psi.get_nbasis();     /// leading dimension of psi
+        const int ld_psi = psi.get_nbasis();     /// leading dimension of psi
 
         // Davidson matrix-blockvector functions
 
         auto ngk_pointer = psi.get_ngk_pointer();
         /// wrap hpsi into lambda function, Matrix \times blockvector
-        /// hpsi(X, HX, nband, dim, band_index1, band_index2)
-        auto hpsi_func = [hm, ngk_pointer](T* psi_in,
-                                           T* hpsi_out,
-                                           const int nband_in,
-                                           const int nbasis_in,
-                                           const int band_index1,
-                                           const int band_index2) {
+        // hpsi_func (X, HX, ld, nvec) -> HX = H(X), X and HX blockvectors of size ld x nvec
+        auto hpsi_func = [hm, ngk_pointer](T *psi_in,
+                                           T *hpsi_out,
+                                           const int ld_psi,
+                                           const int nvec) {
             ModuleBase::timer::tick("David", "hpsi_func");
 
-            // Convert "pointer data stucture" to a psi::Psi object
-            auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nband_in, nbasis_in, ngk_pointer);
+            // Convert pointer of psi_in to a psi::Psi object
+            auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nvec, ld_psi, ngk_pointer);
 
-            psi::Range bands_range(true, 0, band_index1, band_index2);
+            psi::Range bands_range(true, 0, 0, nvec-1);
 
             using hpsi_info = typename hamilt::Operator<T, Device>::hpsi_info;
             hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
@@ -515,14 +512,16 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
 
         /// wrap spsi into lambda function, Matrix \times blockvector
         /// spsi(X, SX, nrow, npw, nbands)
+        /// nrow is leading dimension of spsi, npw is leading dimension of psi, nbands is number of vecs
         auto spsi_func = [hm](const T* psi_in, T* spsi_out,
-                               const int nrow,  // dimension of spsi: nbands * nrow
-                               const int npw,   // number of plane waves
-                               const int nbands // number of bands
+                               const int ld_spsi,  // Leading dimension of spsi. Dimension of SX: nbands * nrow.
+                               const int ld_psi,   // Leading dimension of psi. Number of plane waves.
+                               const int nvec      // Number of vectors(bands)
                             ){
             ModuleBase::timer::tick("David", "spsi_func");
             // sPsi determines S=I or not by GlobalV::use_uspp inside
-            hm->sPsi(psi_in, spsi_out, nrow, npw, nbands);
+            // sPsi(psi, spsi, nrow, npw, nbands)
+            hm->sPsi(psi_in, spsi_out, ld_spsi, ld_psi, nvec);
             ModuleBase::timer::tick("David", "spsi_func");
         };
 
@@ -535,7 +534,7 @@ void HSolverPW<T, Device>::hamiltSolvePsiK(hamilt::Hamilt<T, Device>* hm,
         // do diag and add davidson iteration counts up to avg_iter
         DiagoIterAssist<T, Device>::avg_iter += static_cast<double>(david.diag(hpsi_func,
                                                                                spsi_func,
-                                                                               ldPsi,
+                                                                               ld_psi,
                                                                                psi.get_pointer(),
                                                                                eigenvalue,
                                                                                david_diag_thr,
diff --git a/source/module_hsolver/test/diago_david_float_test.cpp b/source/module_hsolver/test/diago_david_float_test.cpp
index 20d50d4a60..69b9b898c3 100644
--- a/source/module_hsolver/test/diago_david_float_test.cpp
+++ b/source/module_hsolver/test/diago_david_float_test.cpp
@@ -89,7 +89,7 @@ class DiagoDavPrepare
 
 		const int dim = phi.get_current_nbas() ;
 		const int nband = phi.get_nbands();
-		const int ldPsi =phi.get_nbasis();
+		const int ld_psi =phi.get_nbasis();
 		hsolver::DiagoDavid<std::complex<float>> dav(precondition, nband, dim, order, false, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<float>>::PW_DIAG_NMAX = maxiter;
@@ -108,11 +108,10 @@ class DiagoDavPrepare
 
 		
 		auto hpsi_func = [phm](std::complex<float>* psi_in,std::complex<float>* hpsi_out,
-					const int nband_in, const int nbasis_in,
-                    const int band_index1, const int band_index2)
+					const int ld_psi, const int nvec)
                     {
-                        auto psi_iter_wrapper = psi::Psi<std::complex<float>>(psi_in, 1, nband_in, nbasis_in, nullptr);
-                        psi::Range bands_range(1, 0, band_index1, band_index2);
+                        auto psi_iter_wrapper = psi::Psi<std::complex<float>>(psi_in, 1, nvec, ld_psi, nullptr);
+                        psi::Range bands_range(1, 0, 0, nvec-1);
                         using hpsi_info = typename hamilt::Operator<std::complex<float>>::hpsi_info;
                         hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
                         phm->ops->hPsi(info);
@@ -120,7 +119,7 @@ class DiagoDavPrepare
 		auto spsi_func = [phm](const std::complex<float>* psi_in, std::complex<float>* spsi_out,const int nrow, const int npw,  const int nbands){
 			phm->sPsi(psi_in, spsi_out, nrow, npw, nbands);
 		};
-		dav.diag(hpsi_func,spsi_func, ldPsi, phi.get_pointer(), en, eps, maxiter);
+		dav.diag(hpsi_func,spsi_func, ld_psi, phi.get_pointer(), en, eps, maxiter);
 
 #ifdef __MPI		
 		end = MPI_Wtime();
diff --git a/source/module_hsolver/test/diago_david_real_test.cpp b/source/module_hsolver/test/diago_david_real_test.cpp
index a5be15ba70..d8585210d2 100644
--- a/source/module_hsolver/test/diago_david_real_test.cpp
+++ b/source/module_hsolver/test/diago_david_real_test.cpp
@@ -48,8 +48,10 @@ void lapackEigen(int& npw, std::vector<double>& hm, double* e, bool outtime = fa
     double* work2 = new double[lwork];
     dsyev_(&tmp_c1, &tmp_c2, &npw, tmp.data(), &npw, e, work2, &lwork, &info);
     end = clock();
-    if (info) std::cout << "ERROR: Lapack solver, info=" << info << std::endl;
-    if (outtime) std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl;
+    if (info) { std::cout << "ERROR: Lapack solver, info=" << info << std::endl;
+}
+    if (outtime) { std::cout << "Lapack Run time: " << (double)(end - start) / CLOCKS_PER_SEC << " S" << std::endl;
+}
     delete[] work2;
 }
 class DiagoDavPrepare
@@ -73,7 +75,8 @@ class DiagoDavPrepare
         //calculate eigenvalues by LAPACK;
         double* e_lapack = new double[npw];
         double* ev;
-        if (mypnum == 0) lapackEigen(npw, DIAGOTEST::hmatrix_d, e_lapack, DETAILINFO);
+        if (mypnum == 0) { lapackEigen(npw, DIAGOTEST::hmatrix_d, e_lapack, DETAILINFO);
+}
 
         //do Diago_David::diag()
         double* en = new double[npw];
@@ -88,7 +91,7 @@ class DiagoDavPrepare
 
 		const int dim = phi.get_current_nbas();
         const int nband = phi.get_nbands();
-        const int ldPsi = phi.get_nbasis();
+        const int ld_psi = phi.get_nbasis();
         hsolver::DiagoDavid<double> dav(precondition, nband, dim, order, false, comm_info);
 
         hsolver::DiagoIterAssist<double>::PW_DIAG_NMAX = maxiter;
@@ -107,11 +110,10 @@ class DiagoDavPrepare
 
         
         auto hpsi_func = [phm](double* psi_in,double* hpsi_out,
-					const int nband_in, const int nbasis_in,
-                    const int band_index1, const int band_index2)
+					const int ld_psi, const int nvec)
                     {
-                        auto psi_iter_wrapper = psi::Psi<double>(psi_in, 1, nband_in, nbasis_in, nullptr);
-                        psi::Range bands_range(1, 0, band_index1, band_index2);
+                        auto psi_iter_wrapper = psi::Psi<double>(psi_in, 1, nvec, ld_psi, nullptr);
+                        psi::Range bands_range(true, 0, 0, nvec-1);
                         using hpsi_info = typename hamilt::Operator<double>::hpsi_info;
                         hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
                         phm->ops->hPsi(info);
@@ -119,7 +121,7 @@ class DiagoDavPrepare
         auto spsi_func = [phm](const double* psi_in, double* spsi_out,const int nrow, const int npw,  const int nbands){
 			phm->sPsi(psi_in, spsi_out, nrow, npw, nbands);
 		};
-        dav.diag(hpsi_func,spsi_func, ldPsi, phi.get_pointer(), en, eps, maxiter);
+        dav.diag(hpsi_func,spsi_func, ld_psi, phi.get_pointer(), en, eps, maxiter);
 
 #ifdef __MPI		
         end = MPI_Wtime();
@@ -131,7 +133,8 @@ class DiagoDavPrepare
 
         if (mypnum == 0)
         {
-            if (DETAILINFO) std::cout << "diag Run time: " << use_time << std::endl;
+            if (DETAILINFO) { std::cout << "diag Run time: " << use_time << std::endl;
+}
             for (int i = 0;i < nband;i++)
             {
                 EXPECT_NEAR(en[i], e_lapack[i], CONVTHRESHOLD);
@@ -148,8 +151,9 @@ class DiagoDavTest : public ::testing::TestWithParam<DiagoDavPrepare> {};
 TEST_P(DiagoDavTest, RandomHamilt)
 {
     DiagoDavPrepare ddp = GetParam();
-    if (DETAILINFO && ddp.mypnum == 0) std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity="
+    if (DETAILINFO && ddp.mypnum == 0) { std::cout << "npw=" << ddp.npw << ", nband=" << ddp.nband << ", sparsity="
         << ddp.sparsity << ", eps=" << ddp.eps << std::endl;
+}
 
     HPsi<double> hpsi(ddp.nband, ddp.npw, ddp.sparsity);
     DIAGOTEST::hmatrix_d = hpsi.hamilt();
@@ -236,7 +240,8 @@ int main(int argc, char** argv)
 
     testing::InitGoogleTest(&argc, argv);
     ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
-    if (myrank != 0) delete listeners.Release(listeners.default_result_printer());
+    if (myrank != 0) { delete listeners.Release(listeners.default_result_printer());
+}
 
     int result = RUN_ALL_TESTS();
     if (myrank == 0 && result != 0)
diff --git a/source/module_hsolver/test/diago_david_test.cpp b/source/module_hsolver/test/diago_david_test.cpp
index 79b099c6f4..03fa7947a5 100644
--- a/source/module_hsolver/test/diago_david_test.cpp
+++ b/source/module_hsolver/test/diago_david_test.cpp
@@ -91,7 +91,7 @@ class DiagoDavPrepare
 
 		const int dim = phi.get_current_nbas();
 		const int nband = phi.get_nbands();
-		const int ldPsi = phi.get_nbasis();
+		const int ld_psi = phi.get_nbasis();
 		hsolver::DiagoDavid<std::complex<double>> dav(precondition, nband, dim, order, false, comm_info);
 
 		hsolver::DiagoIterAssist<std::complex<double>>::PW_DIAG_NMAX = maxiter;
@@ -110,11 +110,10 @@ class DiagoDavPrepare
 
 		
 		auto hpsi_func = [phm](std::complex<double>* psi_in,std::complex<double>* hpsi_out,
-					const int nband_in, const int nbasis_in,
-                    const int band_index1, const int band_index2)
+					const int ld_psi, const int nvec)
                     {
-                        auto psi_iter_wrapper = psi::Psi<std::complex<double>>(psi_in, 1, nband_in, nbasis_in, nullptr);
-                        psi::Range bands_range(1, 0, band_index1, band_index2);
+                        auto psi_iter_wrapper = psi::Psi<std::complex<double>>(psi_in, 1, nvec, ld_psi, nullptr);
+                        psi::Range bands_range(1, 0, 0, nvec-1);
                         using hpsi_info = typename hamilt::Operator<std::complex<double>>::hpsi_info;
                         hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
                         phm->ops->hPsi(info);
@@ -122,7 +121,7 @@ class DiagoDavPrepare
 		auto spsi_func = [phm](const std::complex<double>* psi_in, std::complex<double>* spsi_out,const int nrow, const int npw,  const int nbands){
 			phm->sPsi(psi_in, spsi_out, nrow, npw, nbands);
 		};
-		dav.diag(hpsi_func,spsi_func, ldPsi, phi.get_pointer(), en, eps, maxiter);
+		dav.diag(hpsi_func,spsi_func, ld_psi, phi.get_pointer(), en, eps, maxiter);
 
 #ifdef __MPI		
 		end = MPI_Wtime();
diff --git a/source/module_hsolver/test/hsolver_pw_sup.h b/source/module_hsolver/test/hsolver_pw_sup.h
index bca46f585b..fa5cf58610 100644
--- a/source/module_hsolver/test/hsolver_pw_sup.h
+++ b/source/module_hsolver/test/hsolver_pw_sup.h
@@ -153,9 +153,9 @@ DiagoDavid<T, Device>::~DiagoDavid() {
 }
 
 template <typename T, typename Device>
-int DiagoDavid<T, Device>::diag(const std::function<void(T*, T*, const int, const int, const int, const int)>& hpsi_func,
+int DiagoDavid<T, Device>::diag(const std::function<void(T*, T*, const int, const int)>& hpsi_func,
                                 const std::function<void(T*, T*, const int, const int, const int)>& spsi_func,
-                                const int ldPsi,
+                                const int ld_psi,
                                 T *psi_in,
                                 Real* eigenvalue_in,
                                 const Real david_diag_thr,
diff --git a/source/module_lr/hsolver_lrtd.cpp b/source/module_lr/hsolver_lrtd.cpp
index 0cc95fb8f6..0e954a2cd9 100644
--- a/source/module_lr/hsolver_lrtd.cpp
+++ b/source/module_lr/hsolver_lrtd.cpp
@@ -14,7 +14,8 @@ namespace LR
     inline void print_eigs(const std::vector<T>& eigs, const std::string& label = "", const double factor = 1.0)
     {
         std::cout << label << std::endl;
-        for (auto& e : eigs)std::cout << e * factor << " ";
+        for (auto& e : eigs) {std::cout << e * factor << " ";
+}
         std::cout << std::endl;
     }
     template<typename T, typename Device>
@@ -84,13 +85,11 @@ namespace LR
                 auto hpsi_func = [pHamilt](
                     T* psi_in,
                     T* hpsi_out,
-                    const int nband_in,
-                    const int nbasis_in,
-                    const int band_index1,
-                    const int band_index2)
+                    const int ld_psi,
+                    const int nvec)
                     {
-                        auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nband_in, nbasis_in, nullptr);
-                        psi::Range bands_range(true, 0, band_index1, band_index2);
+                        auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nvec, ld_psi, nullptr);
+                        psi::Range bands_range(true, 0, 0, nvec-1);
                         using hpsi_info = typename hamilt::Operator<T, Device>::hpsi_info;
                         hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
                         pHamilt->ops->hPsi(info);
@@ -118,16 +117,14 @@ namespace LR
                     false, //always do the subspace diag (check the implementation)
                     comm_info);
 
-                std::function<void(T*, T*, const int, const int, const int, const int)> hpsi_func = [pHamilt](
+                auto hpsi_func = [pHamilt](
                     T* psi_in,
                     T* hpsi_out,
-                    const int nband_in,
-                    const int nbasis_in,
-                    const int band_index1,
-                    const int band_index2)
+                    const int ld_psi,
+                    const int nvec)
                     {
-                        auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nband_in, nbasis_in, nullptr);
-                        psi::Range bands_range(true, 0, band_index1, band_index2);
+                        auto psi_iter_wrapper = psi::Psi<T, Device>(psi_in, 1, nvec, ld_psi, nullptr);
+                        psi::Range bands_range(true, 0, 0, nvec-1);
                         using hpsi_info = typename hamilt::Operator<T, Device>::hpsi_info;
                         hpsi_info info(&psi_iter_wrapper, bands_range, hpsi_out);
                         pHamilt->ops->hPsi(info);