diff --git a/source/source_hsolver/diago_cusolvermp.cpp b/source/source_hsolver/diago_cusolvermp.cpp
index b152066c56f..ec6ce1e7917 100644
--- a/source/source_hsolver/diago_cusolvermp.cpp
+++ b/source/source_hsolver/diago_cusolvermp.cpp
@@ -12,11 +12,12 @@ using complex = std::complex<double>;
 namespace hsolver
 {
 template <typename T>
-void DiagoCusolverMP<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in)
+void DiagoCusolverMP<T>::diag(ModuleBase::MatrixBlock<T>& h_mat,
+                              ModuleBase::MatrixBlock<T>& s_mat,
+                              psi::Psi<T>& psi,
+                              Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoCusolverMP", "diag");
-    ModuleBase::MatrixBlock<T> h_mat, s_mat;
-    phm_in->matrix(h_mat, s_mat);
 
     std::vector<Real> eigen(this->nlocal, 0.0);
     std::vector<T> eigenvectors(h_mat.row * h_mat.col);
diff --git a/source/source_hsolver/diago_cusolvermp.h b/source/source_hsolver/diago_cusolvermp.h
index 155a4bf94d3..3794cf1c084 100644
--- a/source/source_hsolver/diago_cusolvermp.h
+++ b/source/source_hsolver/diago_cusolvermp.h
@@ -2,10 +2,11 @@
 #define DIAGO_CUSOLVERMPH
 
 #ifdef __CUSOLVERMP
-#include "source_hamilt/hamilt.h"
 #include "source_base/macros.h"
+#include "source_base/matrix_block.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
 #include "source_hsolver/kernels/cuda/diag_cusolvermp.cuh"
+#include "source_psi/psi.h"
 namespace hsolver
 {
 // DiagoCusolverMP class, for diagonalization using CUSOLVERMP
@@ -22,7 +23,10 @@ class DiagoCusolverMP
     {
     }
     // the diag function for CUSOLVERMP diagonalization
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
 
   private:
     const int nlocal;
diff --git a/source/source_hsolver/diago_elpa.cpp b/source/source_hsolver/diago_elpa.cpp
index 1673bc8091d..eceb85da86c 100644
--- a/source/source_hsolver/diago_elpa.cpp
+++ b/source/source_hsolver/diago_elpa.cpp
@@ -7,9 +7,6 @@
 #include "source_base/tool_title.h"
 #include "source_base/tool_quit.h"
 
-typedef ModuleBase::MatrixBlock<double> matd;
-typedef ModuleBase::MatrixBlock<std::complex<double>> matcd;
-
 namespace hsolver {
 #ifdef __MPI
 template <>
@@ -65,14 +62,12 @@ MPI_Comm DiagoElpa<std::complex<double>>::setmpicomm() {
 #endif
 template <>
 void DiagoElpa<std::complex<double>>::diag(
-    hamilt::Hamilt<std::complex<double>>* phm_in,
+    ModuleBase::MatrixBlock<std::complex<double>>& h_mat,
+    ModuleBase::MatrixBlock<std::complex<double>>& s_mat,
     psi::Psi<std::complex<double>>& psi,
     Real* eigenvalue_in) {
     ModuleBase::TITLE("DiagoElpa", "diag");
 #ifdef __MPI
-    matcd h_mat, s_mat;
-    phm_in->matrix(h_mat, s_mat);
-
     std::vector<double> eigen(this->nlocal, 0.0);
 
     bool isReal = false;
@@ -103,14 +98,12 @@ void DiagoElpa<std::complex<double>>::diag(
 }
 
 template <>
-void DiagoElpa<double>::diag(hamilt::Hamilt<double>* phm_in,
+void DiagoElpa<double>::diag(ModuleBase::MatrixBlock<double>& h_mat,
+                             ModuleBase::MatrixBlock<double>& s_mat,
                              psi::Psi<double>& psi,
                              Real* eigenvalue_in) {
     ModuleBase::TITLE("DiagoElpa", "diag");
 #ifdef __MPI
-    matd h_mat, s_mat;
-    phm_in->matrix(h_mat, s_mat);
-
     std::vector<double> eigen(this->nlocal, 0.0);
 
     bool isReal = true;
diff --git a/source/source_hsolver/diago_elpa.h b/source/source_hsolver/diago_elpa.h
index 09e0dd55e4a..e7e528f380f 100644
--- a/source/source_hsolver/diago_elpa.h
+++ b/source/source_hsolver/diago_elpa.h
@@ -3,8 +3,8 @@
 
 #include "source_base/macros.h"   // GetRealType
 #include "source_base/matrix_block.h"
-#include "source_hamilt/hamilt.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
+#include "source_psi/psi.h"
 
 namespace hsolver
 {
@@ -20,7 +20,10 @@ class DiagoElpa
     /// @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);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
     void diag_pool(ModuleBase::MatrixBlock<T>& h_mat, ModuleBase::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
diff --git a/source/source_hsolver/diago_elpa_native.cpp b/source/source_hsolver/diago_elpa_native.cpp
index a17c192137f..a07f8303202 100644
--- a/source/source_hsolver/diago_elpa_native.cpp
+++ b/source/source_hsolver/diago_elpa_native.cpp
@@ -143,12 +143,13 @@ void DiagoElpaNative<T>::diag_pool(ModuleBase::MatrixBlock<T>& h_mat,
 #endif
 
 template <typename T>
-void DiagoElpaNative<T>::diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in)
+void DiagoElpaNative<T>::diag(ModuleBase::MatrixBlock<T>& h_mat,
+                              ModuleBase::MatrixBlock<T>& s_mat,
+                              psi::Psi<T>& psi,
+                              Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoElpaNative", "diag");
 #ifdef __MPI
-    ModuleBase::MatrixBlock<T> h_mat, s_mat;
-    phm_in->matrix(h_mat, s_mat);
     MPI_Comm COMM_DIAG = setmpicomm(); // set mpi_comm needed
     diag_pool(h_mat, s_mat, psi, eigenvalue_in, COMM_DIAG);
 #else
diff --git a/source/source_hsolver/diago_elpa_native.h b/source/source_hsolver/diago_elpa_native.h
index 5511c8d8956..b938240299d 100644
--- a/source/source_hsolver/diago_elpa_native.h
+++ b/source/source_hsolver/diago_elpa_native.h
@@ -3,8 +3,8 @@
 
 #include "source_base/macros.h"   // GetRealType
 #include "source_base/matrix_block.h"
-#include "source_hamilt/hamilt.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
+#include "source_psi/psi.h"
 
 namespace hsolver
 {
@@ -22,7 +22,10 @@ class DiagoElpaNative
     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);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
     void diag_pool(ModuleBase::MatrixBlock<T>& h_mat, ModuleBase::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
diff --git a/source/source_hsolver/diago_lapack.cpp b/source/source_hsolver/diago_lapack.cpp
index 77d42679f0a..48e9651d067 100644
--- a/source/source_hsolver/diago_lapack.cpp
+++ b/source/source_hsolver/diago_lapack.cpp
@@ -8,9 +8,6 @@
 
 #include <cstring>
 
-typedef ModuleBase::MatrixBlock<double> matd;
-typedef ModuleBase::MatrixBlock<std::complex<double>> matcd;
-
 namespace hsolver
 {
 namespace
@@ -32,13 +29,12 @@ void check_lapack_layout(const ModuleBase::MatrixBlock<T>& h_mat,
 }
 } // namespace
 template <>
-void DiagoLapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>& psi, Real* eigenvalue_in)
+void DiagoLapack<double>::diag(ModuleBase::MatrixBlock<double>& h_mat,
+                               ModuleBase::MatrixBlock<double>& s_mat,
+                               psi::Psi<double>& psi,
+                               Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoLapack", "diag");
-    // Prepare H and S matrix
-    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(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
@@ -51,13 +47,12 @@ void DiagoLapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>&
 }
 
 template <>
-void DiagoLapack<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>>* phm_in,
+void DiagoLapack<std::complex<double>>::diag(ModuleBase::MatrixBlock<std::complex<double>>& h_mat,
+                                             ModuleBase::MatrixBlock<std::complex<double>>& s_mat,
                                              psi::Psi<std::complex<double>>& psi,
                                              Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoLapack", "diag");
-    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(this->nlocal, 0.0);
     check_lapack_layout(h_mat, s_mat, eigen.size());
diff --git a/source/source_hsolver/diago_lapack.h b/source/source_hsolver/diago_lapack.h
index 213d772dbcd..75f98dfc0c0 100644
--- a/source/source_hsolver/diago_lapack.h
+++ b/source/source_hsolver/diago_lapack.h
@@ -10,9 +10,9 @@
 
 #include "source_base/macros.h"   // GetRealType
 #include "source_base/matrix_block.h"
-#include "source_hamilt/hamilt.h"
 #include "source_base/matrix.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
+#include "source_psi/psi.h"
 
 #include <complex>
 #include <utility>
@@ -31,7 +31,10 @@ class DiagoLapack
     /// @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);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
   #ifdef __MPI
     // diagnolization used in parallel-k case
     void diag_pool(ModuleBase::MatrixBlock<T>& h_mat, ModuleBase::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
diff --git a/source/source_hsolver/diago_pexsi.cpp b/source/source_hsolver/diago_pexsi.cpp
index b4fac092bfc..6e310f8e2c9 100644
--- a/source/source_hsolver/diago_pexsi.cpp
+++ b/source/source_hsolver/diago_pexsi.cpp
@@ -9,9 +9,6 @@
 #include "source_basis/module_ao/parallel_orbitals.h"
 #include "module_pexsi/pexsi_solver.h"
 
-typedef ModuleBase::MatrixBlock<double> matd;
-typedef ModuleBase::MatrixBlock<std::complex<double>> matcd;
-
 namespace hsolver
 {
 template <typename T>
@@ -60,11 +57,12 @@ DiagoPexsi<T>::~DiagoPexsi()
 }
 
 template <>
-void DiagoPexsi<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>& psi, double* eigenvalue_in)
+void DiagoPexsi<double>::diag(ModuleBase::MatrixBlock<double>& h_mat,
+                              ModuleBase::MatrixBlock<double>& s_mat,
+                              psi::Psi<double>& psi,
+                              double* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoPEXSI", "diag");
-    matd h_mat, s_mat;
-    phm_in->matrix(h_mat, s_mat);
     int ik = psi.get_current_k();
     this->ps->prepare(this->ParaV->blacs_ctxt,
                       this->ParaV->nb,
@@ -84,7 +82,8 @@ void DiagoPexsi<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>&
 }
 
 template <>
-void DiagoPexsi<std::complex<double>>::diag(hamilt::Hamilt<std::complex<double>>* phm_in,
+void DiagoPexsi<std::complex<double>>::diag(ModuleBase::MatrixBlock<std::complex<double>>& h_mat,
+                                            ModuleBase::MatrixBlock<std::complex<double>>& s_mat,
                                             psi::Psi<std::complex<double>>& psi,
                                             double* eigenvalue_in)
 {
diff --git a/source/source_hsolver/diago_pexsi.h b/source/source_hsolver/diago_pexsi.h
index bbca89f495e..953f2a06571 100644
--- a/source/source_hsolver/diago_pexsi.h
+++ b/source/source_hsolver/diago_pexsi.h
@@ -4,8 +4,9 @@
 #include <vector>
 #include <memory>
 #include "source_base/macros.h"   // GetRealType
-#include "source_hamilt/hamilt.h"
+#include "source_base/matrix_block.h"
 #include "source_basis/module_ao/parallel_orbitals.h"
+#include "source_psi/psi.h"
 #include "module_pexsi/pexsi_solver.h"
 
 namespace hsolver
@@ -24,7 +25,10 @@ class DiagoPexsi
                const int nlocal_in,
                const double nelec_in,
                const int world_nproc_in);
-    void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
     const Parallel_Orbitals* ParaV = nullptr;
     std::vector<T*> DM;
     std::vector<T*> EDM;
diff --git a/source/source_hsolver/diago_scalapack.cpp b/source/source_hsolver/diago_scalapack.cpp
index 32645855af3..644fb83f9a0 100644
--- a/source/source_hsolver/diago_scalapack.cpp
+++ b/source/source_hsolver/diago_scalapack.cpp
@@ -15,9 +15,6 @@
 #include <cassert>
 #include <cstring>
 
-typedef ModuleBase::MatrixBlock<double> matd;
-typedef ModuleBase::MatrixBlock<std::complex<double>> matcd;
-
 namespace hsolver
 {
 namespace
@@ -33,11 +30,12 @@ int blacs_grid_size(const int* const desc)
 } // namespace
 
     template<>
-    void DiagoScalapack<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>& psi, Real* eigenvalue_in)
+    void DiagoScalapack<double>::diag(ModuleBase::MatrixBlock<double>& h_mat,
+    ModuleBase::MatrixBlock<double>& s_mat,
+    psi::Psi<double>& psi,
+    Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoScalapack", "diag");
-    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(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);
@@ -45,11 +43,12 @@ int blacs_grid_size(const int* const desc)
     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)
+    void DiagoScalapack<std::complex<double>>::diag(ModuleBase::MatrixBlock<std::complex<double>>& h_mat,
+    ModuleBase::MatrixBlock<std::complex<double>>& s_mat,
+    psi::Psi<std::complex<double>>& psi,
+    Real* eigenvalue_in)
 {
     ModuleBase::TITLE("DiagoScalapack", "diag");
-    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(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);
diff --git a/source/source_hsolver/diago_scalapack.h b/source/source_hsolver/diago_scalapack.h
index 7cffdf8652f..eb05928dc29 100644
--- a/source/source_hsolver/diago_scalapack.h
+++ b/source/source_hsolver/diago_scalapack.h
@@ -14,7 +14,6 @@
 
 #include "source_base/macros.h"   // GetRealType
 #include "source_base/matrix_block.h"
-#include "source_hamilt/hamilt.h"
 #include "source_psi/psi.h"
 #include "source_base/complexmatrix.h"
 #include "source_base/matrix.h"
@@ -32,7 +31,10 @@ namespace hsolver
     /// @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);
+    void diag(ModuleBase::MatrixBlock<T>& h_mat,
+              ModuleBase::MatrixBlock<T>& s_mat,
+              psi::Psi<T>& psi,
+              Real* eigenvalue_in);
 #ifdef __MPI
     // diagnolization used in parallel-k case
     void diag_pool(ModuleBase::MatrixBlock<T>& h_mat, ModuleBase::MatrixBlock<T>& s_mat, psi::Psi<T>& psi, Real* eigenvalue_in, MPI_Comm& comm);
diff --git a/source/source_hsolver/hsolver_lcao.cpp b/source/source_hsolver/hsolver_lcao.cpp
index 967ed9c153b..79751179d77 100644
--- a/source/source_hsolver/hsolver_lcao.cpp
+++ b/source/source_hsolver/hsolver_lcao.cpp
@@ -120,8 +120,10 @@ void HSolverLCAO<TK>::solve(hamilt::Hamilt<TK>* pHamilt,
             /// update H(k) for each k point
             pHamilt->updateHk(ik);
             psi.fix_k(ik);
+            ModuleBase::MatrixBlock<TK> hk, sk;
+            pHamilt->matrix(hk, sk);
             // solve eigenvector and eigenvalue for H(k)
-            pe.diag(pHamilt, psi, nullptr);
+            pe.diag(hk, sk, psi, nullptr);
         }
         auto _pes = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(pes);
         pes->f_en.eband = pe.totalFreeEnergy;
@@ -140,23 +142,28 @@ void HSolverLCAO<T>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>& psi, do
     ModuleBase::TITLE("HSolverLCAO", "hamiltSolvePsiK");
     ModuleBase::timer::start("HSolverLCAO", "hamiltSolvePsiK");
 
+    // H(k) and S(k) are all the eigensolvers need from the Hamiltonian, so
+    // fetch them once here rather than once inside each solver.
+    ModuleBase::MatrixBlock<T> hk, sk;
+    hm->matrix(hk, sk);
+
     if (this->method == "scalapack_gvx")
     {
 #ifdef __MPI
         DiagoScalapack<T> sa(this->nlocal, this->nbands);
-        sa.diag(hm, psi, eigenvalue);
+        sa.diag(hk, sk, psi, eigenvalue);
 #endif
     }
 #ifdef __ELPA
     else if (this->method == "genelpa")
     {
         DiagoElpa<T> el(this->nlocal, this->nbands);
-        el.diag(hm, psi, eigenvalue);
+        el.diag(hk, sk, psi, eigenvalue);
     }
     else if (this->method == "elpa")
     {
         DiagoElpaNative<T> el(this->nlocal, this->nbands, this->use_gpu);
-        el.diag(hm, psi, eigenvalue);
+        el.diag(hk, sk, psi, eigenvalue);
     }
 #endif
 #ifdef __CUDA
@@ -164,22 +171,20 @@ void HSolverLCAO<T>::hamiltSolvePsiK(hamilt::Hamilt<T>* hm, psi::Psi<T>& psi, do
     {
         // Note: This branch will only be executed in the single-process case
         DiagoCusolver<T> cu(this->nlocal, this->nbands);
-        ModuleBase::MatrixBlock<T> hk, sk;
-        hm->matrix(hk, sk);
         cu.diag(hk, sk, psi, eigenvalue);
     }
 #ifdef __CUSOLVERMP
     else if (this->method == "cusolvermp")
     {
         DiagoCusolverMP<T> cm(this->nlocal, this->nbands);
-        cm.diag(hm, psi, eigenvalue);
+        cm.diag(hk, sk, psi, eigenvalue);
     }
 #endif
 #endif
     else if (this->method == "lapack") // only for single core
     {
         DiagoLapack<T> la(this->nlocal, this->nbands);
-        la.diag(hm, psi, eigenvalue);
+        la.diag(hk, sk, psi, eigenvalue);
     }
     else
     {
@@ -210,6 +215,12 @@ void HSolverLCAO<T>::parakSolve(hamilt::Hamilt<T>* pHamilt,
     int coord_col = k2d.get_p2D_pool()->get_coord_col();
     int ncol_bands_pool
         = numroc_(&(nbands), &(nb2d), &coord_col, &zero, &(k2d.get_p2D_pool()->dim1));
+    /// Parallel_K2D only redistributes H(k)/S(k); updating the Hamiltonian
+    /// for a given k point stays here, where the Hamiltonian is known.
+    auto get_hsk = [pHamilt](int ik, ModuleBase::MatrixBlock<T>& hk, ModuleBase::MatrixBlock<T>& sk) {
+        pHamilt->updateHk(ik);
+        pHamilt->matrix(hk, sk);
+    };
     /// Loop over k points for solve Hamiltonian to charge density
     for (int ik = 0; ik < k2d.get_pKpoints()->get_max_nks_pool(); ++ik)
     {
@@ -235,7 +246,7 @@ void HSolverLCAO<T>::parakSolve(hamilt::Hamilt<T>* pHamilt,
                 ik_kpar[i] = ik + k2d.get_pKpoints()->startk_pool[i];
             }
         }
-        k2d.distribute_hsk(pHamilt, ik_kpar, nrow);
+        k2d.distribute_hsk(get_hsk, ik_kpar, nrow);
         /// global index of k point
         int ik_global = ik + k2d.get_pKpoints()->startk_pool[k2d.get_my_pool()];
         auto psi_pool = psi::Psi<T>(1, ncol_bands_pool, k2d.get_p2D_pool()->nrow, k2d.get_p2D_pool()->nrow, true);
diff --git a/source/source_hsolver/parallel_k2d.cpp b/source/source_hsolver/parallel_k2d.cpp
index b7521fde703..a80205440df 100644
--- a/source/source_hsolver/parallel_k2d.cpp
+++ b/source/source_hsolver/parallel_k2d.cpp
@@ -36,16 +36,15 @@ void Parallel_K2D<TK>::set_para_env(int nks,
 }
 
 template <typename TK>
-void Parallel_K2D<TK>::distribute_hsk(hamilt::Hamilt<TK>* pHamilt,
+void Parallel_K2D<TK>::distribute_hsk(const typename Parallel_K2D<TK>::HskFunc& get_hsk,
                                       const std::vector<int>& ik_kpar,
                                       const int& nw) {
 #ifdef __MPI
     ModuleBase::timer::start("Parallel_K2D", "distribute_hsk");
     for (int ipool = 0; ipool < ik_kpar.size(); ++ipool)
     {
-        pHamilt->updateHk(ik_kpar[ipool]);
         ModuleBase::MatrixBlock<TK> HK_global, SK_global;
-        pHamilt->matrix(HK_global, SK_global);
+        get_hsk(ik_kpar[ipool], HK_global, SK_global);
         if (this->MY_POOL == this->Pkpoints->whichpool[ik_kpar[ipool]]) {
             this->hk_pool.resize(this->P2D_pool->get_local_size(), 0.0);
             this->sk_pool.resize(this->P2D_pool->get_local_size(), 0.0);
diff --git a/source/source_hsolver/parallel_k2d.h b/source/source_hsolver/parallel_k2d.h
index 57e37ef278f..fb4ee104c94 100644
--- a/source/source_hsolver/parallel_k2d.h
+++ b/source/source_hsolver/parallel_k2d.h
@@ -1,12 +1,15 @@
 #ifndef PARALLEL_K2D_H
 #define PARALLEL_K2D_H
 
+#include "source_base/matrix_block.h"
 #include "source_base/parallel_2d.h"
 #include "source_cell/parallel_kpoints.h"
 #ifdef __MPI
 #include "mpi.h"
 #endif
-#include "source_hamilt/hamilt.h"
+
+#include <functional>
+#include <vector>
 
 /***
  * This is a class to realize k-points parallelism in LCAO code.
@@ -32,8 +35,14 @@ class Parallel_K2D {
                     const int& my_rank,
                     const int& nspin);
 
+    /// Supplies H(k) and S(k) for one k point. The caller owns whatever has
+    /// to happen before the blocks are valid (updating the Hamiltonian for
+    /// that k point, for instance); this class only redistributes them.
+    using HskFunc = std::function<
+        void(int ik, ModuleBase::MatrixBlock<TK>& hk, ModuleBase::MatrixBlock<TK>& sk)>;
+
     /// this function distributes the Hk and Sk matrices to hk_pool and sk_pool
-    void distribute_hsk(hamilt::Hamilt<TK>* pHamilt,
+    void distribute_hsk(const HskFunc& get_hsk,
                         const std::vector<int>& ik_kpar,
                         const int& nw);
 
diff --git a/source/source_hsolver/test/diago_lapack_test.cpp b/source/source_hsolver/test/diago_lapack_test.cpp
index e3d358052e1..f5dead8b518 100644
--- a/source/source_hsolver/test/diago_lapack_test.cpp
+++ b/source/source_hsolver/test/diago_lapack_test.cpp
@@ -7,6 +7,7 @@
 #include <cstring>
 #include <vector>
 
+#include "source_base/matrix_block.h"
 #include "source_hsolver/diago_lapack.h"
 
 #define PASSTHRESHOLD 1e-5
@@ -16,8 +17,10 @@
 
 // A hamilt class used for test. It will be removed in the future.
 
+/// Minimal H(k)/S(k) supplier. The LCAO eigensolvers take the matrix blocks
+/// directly, so this test no longer needs a hamilt::Hamilt subclass.
 template <typename T>
-class HamiltTEST : public hamilt::Hamilt<T>
+class HamiltTEST
 {
   public:
     int desc[9];
@@ -25,17 +28,10 @@ class HamiltTEST : public hamilt::Hamilt<T>
     std::vector<T> h_local;
     std::vector<T> s_local;
 
-    void matrix(hamilt::MatrixBlock<T>& hk_in, hamilt::MatrixBlock<T>& sk_in)
-    {
-        hk_in = hamilt::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-        sk_in = hamilt::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-    }
-
-    void constructHamilt(const int iter, const hamilt::MatrixBlock<double> rho)
-    {
-    }
-    void updateHk(const int ik)
+    void matrix(ModuleBase::MatrixBlock<T>& hk_in, ModuleBase::MatrixBlock<T>& sk_in)
     {
+        hk_in = ModuleBase::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
+        sk_in = ModuleBase::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
     }
 };
 
@@ -186,10 +182,12 @@ class DiagoLapackPrepare
         this->pb2d();
         this->print_hs();
 
+        ModuleBase::MatrixBlock<T> h_mat, s_mat;
+        hmtest.matrix(h_mat, s_mat);
         for (int i = 0; i < REPEATRUN; i++)
         {
             hsolver::DiagoLapack<T> dh(nlocal, nbands);
-            dh.diag(&hmtest, psi, e_solver.data());
+            dh.diag(h_mat, s_mat, psi, e_solver.data());
             // dh->diag(&hmtest, psi, e_solver.data());
         }
         // delete dh;
diff --git a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
index 80669e49515..88170826de7 100644
--- a/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_cusolver_test.cpp
@@ -1,3 +1,4 @@
+#include "source_base/matrix_block.h"
 #include "source_hsolver/diago_scalapack.h"
 #include "source_hsolver/test/diago_elpa_utils.h"
 #include "mpi.h"
@@ -31,8 +32,10 @@
  * self-realized functions in source_hsolver/test/diago_elpa_utils.h
  */
 
+/// Minimal H(k)/S(k) supplier. The LCAO eigensolvers take the matrix blocks
+/// directly, so this test no longer needs a hamilt::Hamilt subclass.
 template <typename T>
-class HamiltTEST : public hamilt::Hamilt<T>
+class HamiltTEST
 {
   public:
     int desc[9];
@@ -40,17 +43,10 @@ class HamiltTEST : public hamilt::Hamilt<T>
     std::vector<T> h_local;
     std::vector<T> s_local;
 
-    void matrix(hamilt::MatrixBlock<T>& hk_in, hamilt::MatrixBlock<T>& sk_in)
-    {
-        hk_in = hamilt::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-        sk_in = hamilt::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-    }
-
-    void constructHamilt(const int iter, const hamilt::MatrixBlock<double> rho)
-    {
-    }
-    void updateHk(const int ik)
+    void matrix(ModuleBase::MatrixBlock<T>& hk_in, ModuleBase::MatrixBlock<T>& sk_in)
     {
+        hk_in = ModuleBase::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
+        sk_in = ModuleBase::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
     }
 };
 
@@ -217,17 +213,17 @@ class DiagoPrepare
         {
             hmtest.h_local = this->h_local;
             hmtest.s_local = this->s_local;
+            ModuleBase::MatrixBlock<T> h_mat, s_mat;
+            hmtest.matrix(h_mat, s_mat);
             if (ks_solver == "scalapack_gvx")
             {
                 hsolver::DiagoScalapack<T> dh(nlocal, nbands);
-                dh.diag(&hmtest, psi, e_solver.data());
+                dh.diag(h_mat, s_mat, psi, e_solver.data());
             }
     #ifdef __CUDA
             else if (ks_solver == "cusolver")
                 {
                     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());
                 }
     #endif
diff --git a/source/source_hsolver/test/diago_lcao_test.cpp b/source/source_hsolver/test/diago_lcao_test.cpp
index 60ef9427fd0..f0bbce6574f 100644
--- a/source/source_hsolver/test/diago_lcao_test.cpp
+++ b/source/source_hsolver/test/diago_lcao_test.cpp
@@ -1,3 +1,4 @@
+#include "source_base/matrix_block.h"
 #include "source_hsolver/diago_scalapack.h"
 #include "source_hsolver/diago_lapack.h"
 #include "source_hsolver/test/diago_elpa_utils.h"
@@ -29,8 +30,10 @@
  * self-realized functions in source_hsolver/test/diago_elpa_utils.h
  */
 
+/// Minimal H(k)/S(k) supplier. The LCAO eigensolvers take the matrix blocks
+/// directly, so this test no longer needs a hamilt::Hamilt subclass.
 template <typename T>
-class HamiltTEST : public hamilt::Hamilt<T>
+class HamiltTEST
 {
   public:
     int desc[9];
@@ -38,17 +41,10 @@ class HamiltTEST : public hamilt::Hamilt<T>
     std::vector<T> h_local;
     std::vector<T> s_local;
 
-    void matrix(hamilt::MatrixBlock<T>& hk_in, hamilt::MatrixBlock<T>& sk_in)
-    {
-        hk_in = hamilt::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-        sk_in = hamilt::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-    }
-
-    void constructHamilt(const int iter, const hamilt::MatrixBlock<double> rho)
-    {
-    }
-    void updateHk(const int ik)
+    void matrix(ModuleBase::MatrixBlock<T>& hk_in, ModuleBase::MatrixBlock<T>& sk_in)
     {
+        hk_in = ModuleBase::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
+        sk_in = ModuleBase::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
     }
 };
 
@@ -218,21 +214,23 @@ class DiagoPrepare
         {
             hmtest.h_local = this->h_local;
             hmtest.s_local = this->s_local;
+            ModuleBase::MatrixBlock<T> h_mat, s_mat;
+            hmtest.matrix(h_mat, s_mat);
             if (ks_solver == "scalapack_gvx")
             {
                 hsolver::DiagoScalapack<T> dh(nlocal, nbands);
-                dh.diag(&hmtest, psi, e_solver.data());
+                dh.diag(h_mat, s_mat, psi, e_solver.data());
             }
             else if (ks_solver == "lapack")
             {
                 hsolver::DiagoLapack<T> la(nlocal, nbands);
-                la.diag(&hmtest, psi, e_solver.data());
+                la.diag(h_mat, s_mat, psi, e_solver.data());
             }
     #ifdef __ELPA
             else if (ks_solver == "genelpa")
             {
                 hsolver::DiagoElpa<T> dh(nlocal, nbands);
-                dh.diag(&hmtest, psi, e_solver.data());
+                dh.diag(h_mat, s_mat, psi, e_solver.data());
             }
     #endif
             // dh.diag(&hmtest, psi, e_solver.data());
diff --git a/source/source_hsolver/test/diago_pexsi_test.cpp b/source/source_hsolver/test/diago_pexsi_test.cpp
index 3429e5b3176..14b23d05b0d 100644
--- a/source/source_hsolver/test/diago_pexsi_test.cpp
+++ b/source/source_hsolver/test/diago_pexsi_test.cpp
@@ -1,4 +1,5 @@
 #ifdef __PEXSI
+#include "source_base/matrix_block.h"
 #include "source_hsolver/diago_pexsi.h"
 
 #include "source_base/module_external/scalapack_connector.h"
@@ -24,8 +25,10 @@
 #define PRINT_HS false
 #define REPEATRUN 1
 
+/// Minimal H(k)/S(k) supplier. The LCAO eigensolvers take the matrix blocks
+/// directly, so this test no longer needs a hamilt::Hamilt subclass.
 template <typename T>
-class HamiltTEST : public hamilt::Hamilt<T>
+class HamiltTEST
 {
   public:
     int desc[9];
@@ -33,17 +36,10 @@ class HamiltTEST : public hamilt::Hamilt<T>
     std::vector<T> h_local;
     std::vector<T> s_local;
 
-    void matrix(hamilt::MatrixBlock<T>& hk_in, hamilt::MatrixBlock<T>& sk_in)
-    {
-        hk_in = hamilt::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-        sk_in = hamilt::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
-    }
-
-    void constructHamilt(const int iter, const hamilt::MatrixBlock<double> rho)
-    {
-    }
-    void updateHk(const int ik)
+    void matrix(ModuleBase::MatrixBlock<T>& hk_in, ModuleBase::MatrixBlock<T>& sk_in)
     {
+        hk_in = ModuleBase::MatrixBlock<T>{this->h_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
+        sk_in = ModuleBase::MatrixBlock<T>{this->s_local.data(), (size_t)this->nrow, (size_t)this->ncol, this->desc};
     }
 };
 
@@ -254,7 +250,9 @@ class PexsiPrepare
         {
             hmtest.h_local = this->h_local;
             hmtest.s_local = this->s_local;
-            dh->diag(&hmtest, psi, nullptr);
+            ModuleBase::MatrixBlock<T> h_mat, s_mat;
+            hmtest.matrix(h_mat, s_mat);
+            dh->diag(h_mat, s_mat, psi, nullptr);
 
             // copy the density matrix to dm_local
             dm_local = dh->DM;