diff --git a/source/source_base/complexarray.cpp b/source/source_base/complexarray.cpp
index becb8e84482..17bfa7f9ee4 100644
--- a/source/source_base/complexarray.cpp
+++ b/source/source_base/complexarray.cpp
@@ -85,7 +85,7 @@ void ComplexArray::operator=(const std::complex < double> c){
 	const int size = this->getSize();
 	for (int i = 0; i < size; i++)
 		ptr[i] = c;}
-ComplexArray ComplexArray::operator+(const ComplexArray &cd){
+ComplexArray ComplexArray::operator+(const ComplexArray &cd) const{
 	const int size = this->getSize();
 	assert(size==cd.getSize());
 	ComplexArray cd2(*this);
@@ -98,7 +98,7 @@ void ComplexArray::operator+=(const ComplexArray & cd){
 	for (int i = 0; i < size; i++)
 		ptr[i] += cd.ptr[i];
 }
-ComplexArray ComplexArray::operator-(const ComplexArray &cd){
+ComplexArray ComplexArray::operator-(const ComplexArray &cd) const{
 	const int size = this->getSize();
 	assert(size==cd.getSize());
 	ComplexArray cd2(*this);
@@ -123,7 +123,7 @@ ComplexArray operator*(const double r, const ComplexArray &cd){
 	for (int i = 0; i < size; i++)
 		cd2.ptr[i] *= r;
 	return cd2;}
-ComplexArray ComplexArray::operator*(const double r){
+ComplexArray ComplexArray::operator*(const double r) const{
 	ComplexArray cd2(*this);
 	const int size = this->getSize();
 	for (int i = 0; i < size; i++)
@@ -135,7 +135,7 @@ ComplexArray operator*(const std::complex < double> c, const ComplexArray &cd){
 	for (int i = 0; i < size; i++)
 		cd2.ptr[i] = c * cd.ptr[i];
 	return cd2;}
-ComplexArray ComplexArray::operator*(const std::complex < double> c){
+ComplexArray ComplexArray::operator*(const std::complex < double> c) const{
 	const int size = this->getSize();
 	ComplexArray cd(size);
 	for (int i = 0; i < size; i++)
diff --git a/source/source_base/complexarray.h b/source/source_base/complexarray.h
index 9cddc9a27cb..831f42ee02d 100644
--- a/source/source_base/complexarray.h
+++ b/source/source_base/complexarray.h
@@ -35,17 +35,17 @@ class ComplexArray
 	/// Assignment of scalar:  all entries set to c.
 	void operator=(std::complex <double> c);
 	/// Add two ComplexArray
-	ComplexArray operator+(const ComplexArray &cd);
+	ComplexArray operator+(const ComplexArray &cd) const;
 	/// Accumulate sum of ComplexArray
 	void operator+=(const ComplexArray &cd);
 	/// Subtract two ComplexArray
-	ComplexArray operator-(const ComplexArray &cd);
+	ComplexArray operator-(const ComplexArray &cd) const;
 	/// Accumulate difference of arrays
 	void operator-=(const ComplexArray &cd);
 	/// Scale a ComplexArray by real r
-	ComplexArray operator*(const double r);
+	ComplexArray operator*(const double r) const;
 	/// Scale a ComplexArray by a std::complex number c
-	ComplexArray operator*(const std::complex <double> c);
+	ComplexArray operator*(const std::complex <double> c) const;
 	/// Scale a ComplexArray by real number in place
 	void operator*=(const double r);
 	/// Scale a ComplexArray by std::complex c in place
diff --git a/source/source_base/complexmatrix.cpp b/source/source_base/complexmatrix.cpp
index 7c5f9a64537..217ebaf0dd7 100644
--- a/source/source_base/complexmatrix.cpp
+++ b/source/source_base/complexmatrix.cpp
@@ -438,7 +438,7 @@ std::ostream & ComplexMatrix::print( std::ostream & os, const double threshold_a
 	return os;
 }
 
-bool ComplexMatrix::checkreal(void)
+bool ComplexMatrix::checkreal(void) const
 {
 	const double tiny = 1e-12;
 	for(int i=0;i<this->nr;i++)
diff --git a/source/source_base/complexmatrix.h b/source/source_base/complexmatrix.h
index e6909c9366f..ac09c3ddb0b 100644
--- a/source/source_base/complexmatrix.h
+++ b/source/source_base/complexmatrix.h
@@ -62,7 +62,7 @@ class ComplexMatrix
 	std::ostream & print( std::ostream & os, const double threshold_abs=0.0, const double threshold_imag=0.0 ) const;		// Peize Lin add 2021.09.08
 
 	// check if all the elements are real
-	bool checkreal(void);
+	bool checkreal(void) const;
 
 	using type=std::complex<double>;					// Peiae Lin add 2022.08.08 for template
 };
diff --git a/source/source_base/formatter.h b/source/source_base/formatter.h
index b795979fb6d..5f465e90a61 100644
--- a/source/source_base/formatter.h
+++ b/source/source_base/formatter.h
@@ -66,7 +66,7 @@ class FmtCore
      *  
      * @return std::string 
      */
-    const std::string& fmt() { return fmt_; }
+    const std::string& fmt() const { return fmt_; }
     /**
      * Python-style string functions will be implemented here as toolbox
      */
diff --git a/source/source_base/opt_CG.h b/source/source_base/opt_CG.h
index 5e14975fcb3..ee1218168e3 100644
--- a/source/source_base/opt_CG.h
+++ b/source/source_base/opt_CG.h
@@ -45,11 +45,11 @@ class Opt_CG
                        int& ifPD        // if postive definit
     );
 
-    double get_residual()
+    double get_residual() const
     {
         return sqrt(this->gg_);
     };
-    int get_iter()
+    int get_iter() const
     {
         return this->iter_;
     }
diff --git a/source/source_basis/module_nao/radial_set.cpp b/source/source_basis/module_nao/radial_set.cpp
index 197f39401b8..eade7423dd8 100644
--- a/source/source_basis/module_nao/radial_set.cpp
+++ b/source/source_basis/module_nao/radial_set.cpp
@@ -138,7 +138,7 @@ void RadialSet::indexing()
     }
 }
 
-const NumericalRadial& RadialSet::chi(const int l, const int izeta)
+const NumericalRadial& RadialSet::chi(const int l, const int izeta) const
 {
     int i = index_map_[l * nzeta_max_ + izeta];
 #ifdef __DEBUG
diff --git a/source/source_basis/module_nao/radial_set.h b/source/source_basis/module_nao/radial_set.h
index 2fc53cc2ee2..f03f7f34b82 100644
--- a/source/source_basis/module_nao/radial_set.h
+++ b/source/source_basis/module_nao/radial_set.h
@@ -159,7 +159,7 @@ class RadialSet
     int nzeta_max() const { return nzeta_max_; }
     int nchi() const { return nchi_; }
 
-    const NumericalRadial& chi(const int l, const int izeta);
+    const NumericalRadial& chi(const int l, const int izeta) const;
     const NumericalRadial* cbegin() const { return chi_; }
     const NumericalRadial* cend() const { return chi_ + nchi_; }
     ///@}
diff --git a/source/source_cell/pseudo.cpp b/source/source_cell/pseudo.cpp
index 8b7ba21ca4c..d2025768212 100644
--- a/source/source_cell/pseudo.cpp
+++ b/source/source_cell/pseudo.cpp
@@ -38,7 +38,7 @@ void pseudo::check_betar()
 	}
 }
 
-void pseudo::print_pseudo(std::ofstream& ofs)
+void pseudo::print_pseudo(std::ofstream& ofs) const
 {
 	print_pseudo_vl(ofs);
 	ofs << "\n pseudo : ";
@@ -50,7 +50,7 @@ void pseudo::print_pseudo(std::ofstream& ofs)
 	ofs << "\n ----------------------";
 }
 
-void pseudo::print_pseudo_atom(std::ofstream &ofs)
+void pseudo::print_pseudo_atom(std::ofstream& ofs) const
 {
 	print_pseudo_h(ofs);
 	ofs << "\n pseudo_atom : ";
@@ -66,7 +66,7 @@ void pseudo::print_pseudo_atom(std::ofstream &ofs)
 }
 
 
-void pseudo::print_pseudo_vl(std::ofstream &ofs)
+void pseudo::print_pseudo_vl(std::ofstream& ofs) const
 {
 	ofs << "\n pseudo_vl:";
 	print_pseudo_atom(ofs);
@@ -74,7 +74,7 @@ void pseudo::print_pseudo_vl(std::ofstream &ofs)
 	ofs << "\n ----------------------------------- ";
 }
 
-void pseudo::print_pseudo_h(std::ofstream &ofs)
+void pseudo::print_pseudo_h(std::ofstream& ofs) const
 {
     ofs << "\n pseudo_info :";
     ofs << "\n nv       " << nv;
diff --git a/source/source_cell/pseudo.h b/source/source_cell/pseudo.h
index 17519ddebaf..abd5be8475b 100644
--- a/source/source_cell/pseudo.h
+++ b/source/source_cell/pseudo.h
@@ -82,10 +82,10 @@ class pseudo
      */
     void check_betar();
 
-    void print_pseudo_h(std::ofstream& ofs);
-    void print_pseudo_atom(std::ofstream& ofs);
-    void print_pseudo_vl(std::ofstream& ofs);
-    void print_pseudo(std::ofstream& ofs);
+    void print_pseudo_h(std::ofstream& ofs) const;
+    void print_pseudo_atom(std::ofstream& ofs) const;
+    void print_pseudo_vl(std::ofstream& ofs) const;
+    void print_pseudo(std::ofstream& ofs) const;
 };
 
 #endif // PSEUDO_H
diff --git a/source/source_cell/sep.cpp b/source/source_cell/sep.cpp
index 30840c40bc7..bd55f3c06b7 100644
--- a/source/source_cell/sep.cpp
+++ b/source/source_cell/sep.cpp
@@ -74,7 +74,7 @@ int SepPot::read_sep(std::ifstream& ifs)
     return 0;
 }
 
-void SepPot::print_sep_info(std::ofstream& ofs)
+void SepPot::print_sep_info(std::ofstream& ofs) const
 {
     ofs << "\n sep_vl:";
     ofs << "\n sep_info:";
@@ -84,7 +84,7 @@ void SepPot::print_sep_info(std::ofstream& ofs)
     ofs << "\n strip electron" << strip_elec;
 }
 
-void SepPot::print_sep_vsep(std::ofstream& ofs)
+void SepPot::print_sep_vsep(std::ofstream& ofs) const
 {
     ofs << "\n mesh  " << mesh;
     output::printr1_d(ofs, " r : ", r, mesh);
diff --git a/source/source_cell/sep.h b/source/source_cell/sep.h
index 20c95b528c7..ae6fb2679c7 100644
--- a/source/source_cell/sep.h
+++ b/source/source_cell/sep.h
@@ -29,8 +29,8 @@ class SepPot
     double* rv = nullptr;   /**< sep potential, but rV, unit: Ry */
 
     int read_sep(std::ifstream& is);
-    void print_sep_info(std::ofstream& ofs);
-    void print_sep_vsep(std::ofstream& ofs);
+    void print_sep_info(std::ofstream& ofs) const;
+    void print_sep_vsep(std::ofstream& ofs) const;
 #ifdef __MPI
     void bcast_sep();
 #endif /* ifdef __MPI */
diff --git a/source/source_cell/test/support/mock_unitcell.cpp b/source/source_cell/test/support/mock_unitcell.cpp
index ce8f7460f3a..67fabe5a9fb 100644
--- a/source/source_cell/test/support/mock_unitcell.cpp
+++ b/source/source_cell/test/support/mock_unitcell.cpp
@@ -40,4 +40,4 @@ void UnitCell::setup(const std::string& latname_in,
 
 void cal_nelec(const Atom* atoms, const int& ntype, double& nelec) {}
 
-void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2) {}
+void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2) const {}
diff --git a/source/source_cell/unitcell.cpp b/source/source_cell/unitcell.cpp
index bae4b1b2dd0..a4966530074 100644
--- a/source/source_cell/unitcell.cpp
+++ b/source/source_cell/unitcell.cpp
@@ -469,7 +469,7 @@ void UnitCell::setup(const std::string& latname_in,
 }
 
 
-void UnitCell::compare_atom_labels(const std::string &label1, const std::string &label2)
+void UnitCell::compare_atom_labels(const std::string& label1, const std::string& label2) const
 {
     if (label1!= label2) //'!( "Ag" == "Ag" || "47" == "47" || "Silver" == Silver" )'
     {
diff --git a/source/source_cell/unitcell.h b/source/source_cell/unitcell.h
index 4b0e702a67a..3292a40a065 100644
--- a/source/source_cell/unitcell.h
+++ b/source/source_cell/unitcell.h
@@ -235,7 +235,7 @@ class UnitCell {
 
     /// @brief check consistency between two atom labels from STRU and pseudo or
     /// orb file
-    void compare_atom_labels(const std::string &label1, const std::string &label2);
+    void compare_atom_labels(const std::string& label1, const std::string& label2) const;
     /// @brief get atomCounts, which is a map from element type to atom number
     std::map<int, int> get_atom_Counts() const;
     /// @brief get orbitalCounts, which is a map from element type to orbital
diff --git a/source/source_esolver/esolver_lj.cpp b/source/source_esolver/esolver_lj.cpp
index d6bb6ff21cb..f71698586ea 100644
--- a/source/source_esolver/esolver_lj.cpp
+++ b/source/source_esolver/esolver_lj.cpp
@@ -126,7 +126,7 @@ void ESolver_LJ::runner(UnitCell& ucell, const int istep)
         GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
     }
 
-    double ESolver_LJ::LJ_energy(const double& d, const int& i, const int& j)
+    double ESolver_LJ::LJ_energy(const double& d, const int& i, const int& j) const
     {
         assert(d > 1e-6); // avoid atom overlap
         const double r2 = d * d;
@@ -135,7 +135,7 @@ void ESolver_LJ::runner(UnitCell& ucell, const int istep)
         return lj_c12(i, j) / (r6 * r6) - lj_c6(i, j) / r6;
     }
 
-    ModuleBase::Vector3<double> ESolver_LJ::LJ_force(const ModuleBase::Vector3<double>& dr, const int& i, const int& j)
+    ModuleBase::Vector3<double> ESolver_LJ::LJ_force(const ModuleBase::Vector3<double>& dr, const int& i, const int& j) const
     {
         const double d = dr.norm();
         assert(d > 1e-6); // avoid atom overlap
diff --git a/source/source_esolver/esolver_lj.h b/source/source_esolver/esolver_lj.h
index d8c53616492..ca23a8bed80 100644
--- a/source/source_esolver/esolver_lj.h
+++ b/source/source_esolver/esolver_lj.h
@@ -27,9 +27,9 @@ namespace ModuleESolver
         void after_all_runners(UnitCell& ucell) override;
 
       private:
-        double LJ_energy(const double& d, const int& i, const int& j);
+        double LJ_energy(const double& d, const int& i, const int& j) const;
 
-        ModuleBase::Vector3<double> LJ_force(const ModuleBase::Vector3<double>& dr, const int& i, const int& j);
+        ModuleBase::Vector3<double> LJ_force(const ModuleBase::Vector3<double>& dr, const int& i, const int& j) const;
 
         void LJ_virial(const ModuleBase::Vector3<double>& force, const ModuleBase::Vector3<double>& dtau);
 
diff --git a/source/source_esolver/esolver_of.h b/source/source_esolver/esolver_of.h
index b94019d5fe8..df4b96543c6 100644
--- a/source/source_esolver/esolver_of.h
+++ b/source/source_esolver/esolver_of.h
@@ -99,7 +99,7 @@ class ESolver_OF : public ESolver_FP
     void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
 
     // --------------------- interface to blas --------------------------
-    double inner_product(double* pa, double* pb, int length, double dV = 1)
+    double inner_product(double* pa, double* pb, int length, double dV = 1) const
     {
         double innerproduct = BlasConnector::dot(length, pa, 1, pb, 1);
         innerproduct *= dV;
diff --git a/source/source_hsolver/diago_cusolver.h b/source/source_hsolver/diago_cusolver.h
index b8cd0fccd22..28ebdc852c6 100644
--- a/source/source_hsolver/diago_cusolver.h
+++ b/source/source_hsolver/diago_cusolver.h
@@ -38,7 +38,7 @@ class DiagoCusolver
   private:
 #ifdef __MPI
     // Function to check if ELPA handle needs to be created or reused in MPI settings
-    bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF);
+    bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const;
 #endif
 };
 
diff --git a/source/source_hsolver/diago_elpa.cpp b/source/source_hsolver/diago_elpa.cpp
index 1e9fac1bc68..24f629f51fd 100644
--- a/source/source_hsolver/diago_elpa.cpp
+++ b/source/source_hsolver/diago_elpa.cpp
@@ -211,7 +211,7 @@ void DiagoElpa<double>::diag_pool(hamilt::MatrixBlock<double>& h_mat,
 
 #ifdef __MPI
 template <typename T>
-bool DiagoElpa<T>::ifElpaHandle(const bool& newIteration, const bool& ifNSCF) {
+bool DiagoElpa<T>::ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const {
     int doHandle = false;
     if (newIteration) {
         doHandle = true;
diff --git a/source/source_hsolver/diago_elpa.h b/source/source_hsolver/diago_elpa.h
index e60ec3febd3..c6050af2245 100644
--- a/source/source_hsolver/diago_elpa.h
+++ b/source/source_hsolver/diago_elpa.h
@@ -27,7 +27,7 @@ class DiagoElpa
 
   private:
 #ifdef __MPI
-    bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF);
+    bool ifElpaHandle(const bool& newIteration, const bool& ifNSCF) const;
     static int lastmpinum; // last using mpi;
 #endif
 };
diff --git a/source/source_hsolver/parallel_k2d.h b/source/source_hsolver/parallel_k2d.h
index f0d0d718cb7..aebd683f4bc 100644
--- a/source/source_hsolver/parallel_k2d.h
+++ b/source/source_hsolver/parallel_k2d.h
@@ -44,15 +44,15 @@ class Parallel_K2D {
     /// set the number of k-points
     void set_kpar(int kpar);
     /// get the number of k-points
-    int get_kpar() { return this->kpar_; }
+    int get_kpar() const { return this->kpar_; }
     /// get my pool
-    int get_my_pool() { return this->MY_POOL; }
+    int get_my_pool() const { return this->MY_POOL; }
     /// get pKpoints
-    Parallel_Kpoints* get_pKpoints() { return this->Pkpoints; }
+    Parallel_Kpoints* get_pKpoints() const { return this->Pkpoints; }
     /// get p2D_global
-    Parallel_2D* get_p2D_global() { return this->P2D_global; }
+    Parallel_2D* get_p2D_global() const { return this->P2D_global; }
     /// get p2D_pool
-    Parallel_2D* get_p2D_pool() { return this->P2D_pool; }
+    Parallel_2D* get_p2D_pool() const { return this->P2D_pool; }
 
     /**
      * the local Hk, Sk matrices in POOL_WORLD_K2D
diff --git a/source/source_io/module_output/csr_reader.cpp b/source/source_io/module_output/csr_reader.cpp
index 4a9c9d79791..019d2715196 100644
--- a/source/source_io/module_output/csr_reader.cpp
+++ b/source/source_io/module_output/csr_reader.cpp
@@ -163,7 +163,7 @@ SparseMatrix<T> csrFileReader<T>::getMatrix(int index) const
 
 // function to get matrix using R coordinate
 template <typename T>
-SparseMatrix<T> csrFileReader<T>::getMatrix(int Rx, int Ry, int Rz)
+SparseMatrix<T> csrFileReader<T>::getMatrix(int Rx, int Ry, int Rz) const
 {
     for (int i = 0; i < RCoordinates.size(); i++)
     {
diff --git a/source/source_io/module_output/csr_reader.h b/source/source_io/module_output/csr_reader.h
index e71bcf8823d..6b38db07f3d 100644
--- a/source/source_io/module_output/csr_reader.h
+++ b/source/source_io/module_output/csr_reader.h
@@ -57,7 +57,7 @@ class csrFileReader : public FileReader
     int getNumberOfR() const;
 
     // get sparse matrix of a specific R coordinate
-    SparseMatrix<T> getMatrix(int Rx, int Ry, int Rz);
+    SparseMatrix<T> getMatrix(int Rx, int Ry, int Rz) const;
 
     // get matrix by using index
     SparseMatrix<T> getMatrix(int index) const;
diff --git a/source/source_io/module_unk/berryphase.cpp b/source/source_io/module_unk/berryphase.cpp
index 04c44a29357..a41ef3f6bd0 100644
--- a/source/source_io/module_unk/berryphase.cpp
+++ b/source/source_io/module_unk/berryphase.cpp
@@ -678,7 +678,7 @@ void berryphase::Macroscopic_polarization(const UnitCell& ucell,
 
 std::string berryphase::outFormat(const double polarization,
                                   const double modulus,
-                                  const ModuleBase::Vector3<double> project)
+                                  const ModuleBase::Vector3<double> project) const
 {
     std::stringstream outStr;
     outStr << std::setw(12) << std::fixed << std::setprecision(7) << polarization << "  (mod ";
diff --git a/source/source_io/module_unk/berryphase.h b/source/source_io/module_unk/berryphase.h
index d7871a028c4..1ba7360dc26 100644
--- a/source/source_io/module_unk/berryphase.h
+++ b/source/source_io/module_unk/berryphase.h
@@ -78,7 +78,7 @@ class berryphase
                                   const ModulePW::PW_Basis_K* wfcpw,
                                   const K_Vectors& kv);
 
-    std::string outFormat(const double polarization, const double modulus, const ModuleBase::Vector3<double> project);
+    std::string outFormat(const double polarization, const double modulus, const ModuleBase::Vector3<double> project) const;
 };
 
 #endif
diff --git a/source/source_io/module_wf/get_wf_lcao.cpp b/source/source_io/module_wf/get_wf_lcao.cpp
index 4fd6c9a5f3c..b896b067bd9 100644
--- a/source/source_io/module_wf/get_wf_lcao.cpp
+++ b/source/source_io/module_wf/get_wf_lcao.cpp
@@ -589,14 +589,14 @@ void Get_wf_lcao::prepare_get_wf(std::ofstream& ofs_running)
     ofs_running << std::setprecision(6);
 }
 
-int Get_wf_lcao::globalIndex(int localindex, int nblk, int nprocs, int myproc)
+int Get_wf_lcao::globalIndex(int localindex, int nblk, int nprocs, int myproc) const
 {
     const int iblock = localindex / nblk;
     const int gIndex = (iblock * nprocs + myproc) * nblk + localindex % nblk;
     return gIndex;
 }
 
-int Get_wf_lcao::localIndex(int globalindex, int nblk, int nprocs, int& myproc)
+int Get_wf_lcao::localIndex(int globalindex, int nblk, int nprocs, int& myproc) const
 {
     myproc = int((globalindex % (nblk * nprocs)) / nblk);
     return int(globalindex / (nblk * nprocs)) * nblk + globalindex % nblk;
diff --git a/source/source_io/module_wf/get_wf_lcao.h b/source/source_io/module_wf/get_wf_lcao.h
index 94afb2cb648..6d01e006fb4 100644
--- a/source/source_io/module_wf/get_wf_lcao.h
+++ b/source/source_io/module_wf/get_wf_lcao.h
@@ -56,9 +56,9 @@ class Get_wf_lcao
                     const double* const* const rho,
                     psi::Psi<std::complex<double>>& wfc_g);
 
-    int globalIndex(int localindex, int nblk, int nprocs, int myproc);
+    int globalIndex(int localindex, int nblk, int nprocs, int myproc) const;
 
-    int localIndex(int globalindex, int nblk, int nprocs, int& myproc);
+    int localIndex(int globalindex, int nblk, int nprocs, int& myproc) const;
 
 #ifdef __MPI
     template <typename T>
diff --git a/source/source_lcao/module_deltaspin/spin_constrain.cpp b/source/source_lcao/module_deltaspin/spin_constrain.cpp
index 6b494879755..6b898f34f6e 100644
--- a/source/source_lcao/module_deltaspin/spin_constrain.cpp
+++ b/source/source_lcao/module_deltaspin/spin_constrain.cpp
@@ -34,7 +34,7 @@ double SpinConstrain<TK>::cal_escon()
 }
 
 template <typename TK>
-double SpinConstrain<TK>::get_escon()
+double SpinConstrain<TK>::get_escon() const
 {
     return this->escon_;
 }
@@ -67,7 +67,7 @@ void SpinConstrain<TK>::set_nspin(int nspin_in)
 
 /// get nspin
 template <typename TK>
-int SpinConstrain<TK>::get_nspin()
+int SpinConstrain<TK>::get_nspin() const
 {
     return this->nspin_;
 }
@@ -439,35 +439,35 @@ void SpinConstrain<TK>::set_input_parameters(double sc_thr_in,
 
 /// get sc_thr
 template <typename TK>
-double SpinConstrain<TK>::get_sc_thr()
+double SpinConstrain<TK>::get_sc_thr() const
 {
     return this->sc_thr_;
 }
 
 /// get nsc
 template <typename TK>
-int SpinConstrain<TK>::get_nsc()
+int SpinConstrain<TK>::get_nsc() const
 {
     return this->nsc_;
 }
 
 /// get nsc_min
 template <typename TK>
-int SpinConstrain<TK>::get_nsc_min()
+int SpinConstrain<TK>::get_nsc_min() const
 {
     return this->nsc_min_;
 }
 
 /// get alpha_trial
 template <typename TK>
-double SpinConstrain<TK>::get_alpha_trial()
+double SpinConstrain<TK>::get_alpha_trial() const
 {
     return this->alpha_trial_;
 }
 
 /// get sccut
 template <typename TK>
-double SpinConstrain<TK>::get_sccut()
+double SpinConstrain<TK>::get_sccut() const
 {
     return this->restrict_current_;
 }
@@ -481,7 +481,7 @@ void SpinConstrain<TK>::set_sc_drop_thr(double sc_drop_thr_in)
 
 /// get sc_drop_thr
 template <typename TK>
-double SpinConstrain<TK>::get_sc_drop_thr()
+double SpinConstrain<TK>::get_sc_drop_thr() const
 {
     return this->sc_drop_thr_;
 }
diff --git a/source/source_lcao/module_deltaspin/spin_constrain.h b/source/source_lcao/module_deltaspin/spin_constrain.h
index c900a771f30..224af123fe4 100644
--- a/source/source_lcao/module_deltaspin/spin_constrain.h
+++ b/source/source_lcao/module_deltaspin/spin_constrain.h
@@ -66,7 +66,7 @@ class SpinConstrain
    */
   double cal_escon();
 
-  double get_escon();
+  double get_escon() const;
 
   void run_lambda_loop(int outer_step, 
 		  bool rerun = true);
@@ -181,7 +181,7 @@ class SpinConstrain
     /// set nspin
     void set_nspin(int nspin);
     /// get nspin
-    int get_nspin();
+    int get_nspin() const;
     /// zero atomic magnetic moment
     void zero_Mi();
     /// get decay_grad
@@ -202,17 +202,17 @@ class SpinConstrain
                               double sccut_in,
                               double sc_drop_thr_in);
     /// get sc_thr
-    double get_sc_thr();
+    double get_sc_thr() const;
     /// get nsc
-    int get_nsc();
+    int get_nsc() const;
     /// get nsc_min
-    int get_nsc_min();
+    int get_nsc_min() const;
     /// get alpha_trial
-    double get_alpha_trial();
+    double get_alpha_trial() const;
     /// get sccut
-    double get_sccut();
+    double get_sccut() const;
     /// get sc_drop_thr
-    double get_sc_drop_thr();
+    double get_sc_drop_thr() const;
     /// @brief set orbital parallel info
     void set_ParaV(Parallel_Orbitals* ParaV_in);
     /// @brief set parameters for solver
diff --git a/source/source_md/msst.cpp b/source/source_md/msst.cpp
index cadea1a8f1a..03bd72aa646 100644
--- a/source/source_md/msst.cpp
+++ b/source/source_md/msst.cpp
@@ -236,7 +236,7 @@ void MSST::restart(const std::string& global_readin_dir)
     return;
 }
 
-double MSST::vel_sum()
+double MSST::vel_sum() const
 {
     double vsum = 0;
 
diff --git a/source/source_md/msst.h b/source/source_md/msst.h
index 6b420169556..ca87571c4d6 100644
--- a/source/source_md/msst.h
+++ b/source/source_md/msst.h
@@ -30,7 +30,7 @@ class MSST : public MD_base
      *
      * @return the sum of square of velocities
      */
-    double vel_sum();
+    double vel_sum() const;
 
     /**
      * @brief rescale the lattice and velocities
diff --git a/source/source_pw/module_ofdft/kedf_manager.cpp b/source/source_pw/module_ofdft/kedf_manager.cpp
index 401d843baf8..0916ebc25a1 100644
--- a/source/source_pw/module_ofdft/kedf_manager.cpp
+++ b/source/source_pw/module_ofdft/kedf_manager.cpp
@@ -166,7 +166,7 @@ void KEDF_Manager::get_potential(
  *
  * @return kinetic energy
  */
-double KEDF_Manager::get_energy()
+double KEDF_Manager::get_energy() const
 {
     double kinetic_energy = 0.0;
 
diff --git a/source/source_pw/module_ofdft/kedf_manager.h b/source/source_pw/module_ofdft/kedf_manager.h
index 8ce607cdd2b..af80646d6ad 100644
--- a/source/source_pw/module_ofdft/kedf_manager.h
+++ b/source/source_pw/module_ofdft/kedf_manager.h
@@ -41,7 +41,7 @@ class KEDF_Manager
         ModuleBase::matrix& rpot
     );
 
-    double get_energy();
+    double get_energy() const;
 
     void get_energy_density(
         const double* const* prho,
diff --git a/source/source_pw/module_ofdft/ml_base.cpp b/source/source_pw/module_ofdft/ml_base.cpp
index 4297083ae43..434de3bc55e 100644
--- a/source/source_pw/module_ofdft/ml_base.cpp
+++ b/source/source_pw/module_ofdft/ml_base.cpp
@@ -132,7 +132,7 @@ void ML_Base::NN_forward(const double * const * prho, const ModulePW::PW_Basis *
     }
 }
 
-torch::Tensor ML_Base::get_data(std::string parameter, const int ikernel){
+torch::Tensor ML_Base::get_data(std::string parameter, const int ikernel) const {
 
     if (parameter == "gamma") return torch::tensor(this->gamma, this->device_type);
     if (parameter == "p") return torch::tensor(this->p, this->device_type);
diff --git a/source/source_pw/module_ofdft/ml_base.h b/source/source_pw/module_ofdft/ml_base.h
index c0ec00c11ed..bace3c1745b 100644
--- a/source/source_pw/module_ofdft/ml_base.h
+++ b/source/source_pw/module_ofdft/ml_base.h
@@ -29,7 +29,7 @@ class ML_Base
     void dumpMatrix(std::string filename, const ModuleBase::matrix &data);
 
     int nx_tot = 0; // equal to nx (called by NN)
-    torch::Tensor get_data(std::string parameter, const int ikernel);
+    torch::Tensor get_data(std::string parameter, const int ikernel) const;
 
 protected:
     void updateInput(const double * const * prho, const ModulePW::PW_Basis *pw_rho);
diff --git a/source/source_pw/module_ofdft/ml_tools/data.cpp b/source/source_pw/module_ofdft/ml_tools/data.cpp
index 75ac968f0f8..d5a658415f1 100644
--- a/source/source_pw/module_ofdft/ml_tools/data.cpp
+++ b/source/source_pw/module_ofdft/ml_tools/data.cpp
@@ -29,7 +29,7 @@ void Data::load_data(Input &input, const int ndata, std::string *dir, const torc
     std::cout << "Load data done" << std::endl;
 }
 
-torch::Tensor Data::get_data(std::string parameter, const int ikernel){
+torch::Tensor Data::get_data(std::string parameter, const int ikernel) const {
     if (parameter == "gamma"){
         return this->gamma.reshape({this->nx_tot});
     }
diff --git a/source/source_pw/module_ofdft/ml_tools/data.h b/source/source_pw/module_ofdft/ml_tools/data.h
index 2116709fd32..dbd04545dec 100644
--- a/source/source_pw/module_ofdft/ml_tools/data.h
+++ b/source/source_pw/module_ofdft/ml_tools/data.h
@@ -60,7 +60,7 @@ class Data
     bool* load_tanhq_nl = nullptr;
 
     void load_data(Input &input, const int ndata, std::string *dir, const torch::Device device);
-    torch::Tensor get_data(std::string parameter, const int ikernel);
+    torch::Tensor get_data(std::string parameter, const int ikernel) const;
 
   private:
     void init_label(Input &input);
diff --git a/source/source_pw/module_ofdft/ml_tools/kernel.cpp b/source/source_pw/module_ofdft/ml_tools/kernel.cpp
index c0d23d6a1f0..53038ef58a6 100644
--- a/source/source_pw/module_ofdft/ml_tools/kernel.cpp
+++ b/source/source_pw/module_ofdft/ml_tools/kernel.cpp
@@ -57,7 +57,7 @@ void Kernel::fill_kernel(const int fftdim,
     std::cout << "Fill kernel done" << std::endl;
 }
 
-double Kernel::wt_kernel(double eta, double tf_weight, double vw_weight)
+double Kernel::wt_kernel(double eta, double tf_weight, double vw_weight) const
 {
     if (eta < 0.) 
     {
@@ -104,7 +104,7 @@ double Kernel::wt_kernel(double eta, double tf_weight, double vw_weight)
     }
 }
 
-double Kernel::yukawa_kernel(double eta, double alpha)
+double Kernel::yukawa_kernel(double eta, double alpha) const
 {
     return (eta == 0 && alpha == 0) ? 0. : M_PI / (eta * eta + alpha * alpha / 4.);
 }
diff --git a/source/source_pw/module_ofdft/ml_tools/kernel.h b/source/source_pw/module_ofdft/ml_tools/kernel.h
index 7199f2a6eff..e9e614c5225 100644
--- a/source/source_pw/module_ofdft/ml_tools/kernel.h
+++ b/source/source_pw/module_ofdft/ml_tools/kernel.h
@@ -30,8 +30,8 @@ class Kernel
                      const std::string *cell,
                      const torch::Device device,
                      const std::vector<torch::Tensor> &fft_gg);
-    double wt_kernel(double eta, double tf_weight = 1., double vw_weight = 1.);
-    double yukawa_kernel(double eta, double alpha);
+    double wt_kernel(double eta, double tf_weight = 1., double vw_weight = 1.) const;
+    double yukawa_kernel(double eta, double alpha) const;
     void read_kernel(const int fftdim,
                      const int ndata,
                      const torch::Tensor &rho,
diff --git a/source/source_pw/module_pwdft/onsite_proj_tools.h b/source/source_pw/module_pwdft/onsite_proj_tools.h
index f50efd6701f..586a45fd948 100644
--- a/source/source_pw/module_pwdft/onsite_proj_tools.h
+++ b/source/source_pw/module_pwdft/onsite_proj_tools.h
@@ -75,8 +75,8 @@ class Onsite_Proj_tools
     void cal_stress_dspin(int ik, int npm, FPTYPE* stress, const ModuleBase::Vector3<double>* lambda, const FPTYPE* h_wg);
 
 
-    std::complex<FPTYPE>* get_becp() { return becp; }
-    std::complex<FPTYPE>* get_dbecp() { return dbecp; }
+    std::complex<FPTYPE>* get_becp() const { return becp; }
+    std::complex<FPTYPE>* get_dbecp() const { return dbecp; }
 
   private:
     /**
diff --git a/source/source_pw/module_pwdft/soc.cpp b/source/source_pw/module_pwdft/soc.cpp
index a44a0b6a48f..66df617a2bc 100644
--- a/source/source_pw/module_pwdft/soc.cpp
+++ b/source/source_pw/module_pwdft/soc.cpp
@@ -41,7 +41,7 @@ Soc::~Soc()
     }
 }
 
-double Soc::spinor(const int l, const double j, const int m, const int spin)
+double Soc::spinor(const int l, const double j, const int m, const int spin) const
 {
     if (spin != 0 && spin != 1)
         ModuleBase::WARNING_QUIT("spinor", "spin direction unknown");
@@ -116,7 +116,7 @@ void Soc::rot_ylm(const int lmax)
     return;
 }
 
-int Soc::sph_ind(const int l, const double j, const int m, const int spin)
+int Soc::sph_ind(const int l, const double j, const int m, const int spin) const
 {
     // This function calculates the m index of the spherical harmonic
     // in a spinor with orbital angular momentum l, total angular
diff --git a/source/source_pw/module_pwdft/soc.h b/source/source_pw/module_pwdft/soc.h
index 282b77f604f..25a56a6429d 100644
--- a/source/source_pw/module_pwdft/soc.h
+++ b/source/source_pw/module_pwdft/soc.h
@@ -50,9 +50,9 @@ class Soc
     Soc(){};
     ~Soc();
 
-    double spinor(const int l, const double j, const int m, const int spin);
+    double spinor(const int l, const double j, const int m, const int spin) const;
 
-    int sph_ind(const int l, const double j, const int m, const int spin);
+    int sph_ind(const int l, const double j, const int m, const int spin) const;
 
     void rot_ylm(const int lmax);
     // std::complex<double> **rotylm;
diff --git a/source/source_pw/module_pwdft/structure_factor.h b/source/source_pw/module_pwdft/structure_factor.h
index 6de0a5b0e37..ae10a8d8a9b 100644
--- a/source/source_pw/module_pwdft/structure_factor.h
+++ b/source/source_pw/module_pwdft/structure_factor.h
@@ -60,7 +60,7 @@ class Structure_Factor
                                   int it,
                                   int ia,
                                   const ModulePW::PW_Basis_K* wfc_basis,
-                                  ModuleBase::Vector3<double> q);
+                                  ModuleBase::Vector3<double> q) const;
 
   private:
 
diff --git a/source/source_pw/module_pwdft/structure_factor_k.cpp b/source/source_pw/module_pwdft/structure_factor_k.cpp
index 3ca4980c580..b5d9c282bee 100644
--- a/source/source_pw/module_pwdft/structure_factor_k.cpp
+++ b/source/source_pw/module_pwdft/structure_factor_k.cpp
@@ -147,7 +147,7 @@ std::complex<double>* Structure_Factor::get_skq(int ik,
                                                 const int it,
                                                 const int ia,
                                                 const ModulePW::PW_Basis_K* wfc_basis,
-                                                ModuleBase::Vector3<double> q) // pengfei 2016-11-23
+                                                ModuleBase::Vector3<double> q) const // pengfei 2016-11-23
 {
     const int npw = wfc_basis->npwk[ik];
     std::complex<double> *skq = new std::complex<double>[npw];
diff --git a/source/source_pw/module_stodft/sto_func.cpp b/source/source_pw/module_stodft/sto_func.cpp
index 255e4a77f05..eb30289b3cc 100644
--- a/source/source_pw/module_stodft/sto_func.cpp
+++ b/source/source_pw/module_stodft/sto_func.cpp
@@ -21,7 +21,7 @@ void Sto_Func<REAL>::set_E_range(REAL* Emin_in, REAL* Emax_in)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::root_fd(REAL e)
+REAL Sto_Func<REAL>::root_fd(REAL e) const
 {
     REAL e_mu = (e - mu) / this->tem;
     if (e_mu > 72) {
@@ -32,7 +32,7 @@ REAL Sto_Func<REAL>::root_fd(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nroot_fd(REAL e)
+REAL Sto_Func<REAL>::nroot_fd(REAL e) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -45,7 +45,7 @@ REAL Sto_Func<REAL>::nroot_fd(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::fd(REAL e)
+REAL Sto_Func<REAL>::fd(REAL e) const
 {
     REAL e_mu = (e - mu) / this->tem;
     if (e_mu > 36) {
@@ -56,7 +56,7 @@ REAL Sto_Func<REAL>::fd(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nfd(REAL e)
+REAL Sto_Func<REAL>::nfd(REAL e) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -69,7 +69,7 @@ REAL Sto_Func<REAL>::nfd(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nxfd(REAL rawe)
+REAL Sto_Func<REAL>::nxfd(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -83,7 +83,7 @@ REAL Sto_Func<REAL>::nxfd(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::fdlnfd(REAL e)
+REAL Sto_Func<REAL>::fdlnfd(REAL e) const
 {
     REAL e_mu = (e - mu) / this->tem;
     if (e_mu > 36) {
@@ -102,7 +102,7 @@ REAL Sto_Func<REAL>::fdlnfd(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nfdlnfd(REAL rawe)
+REAL Sto_Func<REAL>::nfdlnfd(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -123,7 +123,7 @@ REAL Sto_Func<REAL>::nfdlnfd(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::n_root_fdlnfd(REAL rawe)
+REAL Sto_Func<REAL>::n_root_fdlnfd(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -144,7 +144,7 @@ REAL Sto_Func<REAL>::n_root_fdlnfd(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nroot_mfd(REAL rawe)
+REAL Sto_Func<REAL>::nroot_mfd(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -157,7 +157,7 @@ REAL Sto_Func<REAL>::nroot_mfd(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::ncos(REAL rawe)
+REAL Sto_Func<REAL>::ncos(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -166,7 +166,7 @@ REAL Sto_Func<REAL>::ncos(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nsin(REAL rawe)
+REAL Sto_Func<REAL>::nsin(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -175,7 +175,7 @@ REAL Sto_Func<REAL>::nsin(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::n_sin(REAL rawe)
+REAL Sto_Func<REAL>::n_sin(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -184,7 +184,7 @@ REAL Sto_Func<REAL>::n_sin(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::gauss(REAL e)
+REAL Sto_Func<REAL>::gauss(REAL e) const
 {
     REAL a = pow((targ_e - e), 2) / 2.0 / pow(sigma, 2);
     if (a > 72) {
@@ -195,7 +195,7 @@ REAL Sto_Func<REAL>::gauss(REAL e)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::ngauss(REAL rawe)
+REAL Sto_Func<REAL>::ngauss(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
@@ -209,7 +209,7 @@ REAL Sto_Func<REAL>::ngauss(REAL rawe)
 }
 
 template <typename REAL>
-REAL Sto_Func<REAL>::nroot_gauss(REAL rawe)
+REAL Sto_Func<REAL>::nroot_gauss(REAL rawe) const
 {
     REAL Ebar = (*Emin + *Emax) / 2;
     REAL DeltaE = (*Emax - *Emin) / 2;
diff --git a/source/source_pw/module_stodft/sto_func.h b/source/source_pw/module_stodft/sto_func.h
index 2eb67dd50ee..6b2aef43bde 100644
--- a/source/source_pw/module_stodft/sto_func.h
+++ b/source/source_pw/module_stodft/sto_func.h
@@ -14,28 +14,28 @@ class Sto_Func
     void set_E_range(REAL* Emin_in, REAL* Emax_in);
 
   public:
-    REAL root_fd(REAL e);
-    REAL fd(REAL e);
-    REAL nroot_fd(REAL e);
-    REAL nfd(REAL e);
-    REAL nxfd(REAL e);
-    REAL fdlnfd(REAL e);
-    REAL nfdlnfd(REAL e);
-    REAL n_root_fdlnfd(REAL e);
-    REAL nroot_mfd(REAL e);
+    REAL root_fd(REAL e) const;
+    REAL fd(REAL e) const;
+    REAL nroot_fd(REAL e) const;
+    REAL nfd(REAL e) const;
+    REAL nxfd(REAL e) const;
+    REAL fdlnfd(REAL e) const;
+    REAL nfdlnfd(REAL e) const;
+    REAL n_root_fdlnfd(REAL e) const;
+    REAL nroot_mfd(REAL e) const;
 
   public:
     REAL t;
-    REAL ncos(REAL e);
-    REAL nsin(REAL e);
-    REAL n_sin(REAL e);
+    REAL ncos(REAL e) const;
+    REAL nsin(REAL e) const;
+    REAL n_sin(REAL e) const;
 
   public:
     REAL sigma;
     REAL targ_e;
-    REAL gauss(REAL e);
-    REAL ngauss(REAL e);
-    REAL nroot_gauss(REAL e);
+    REAL gauss(REAL e) const;
+    REAL ngauss(REAL e) const;
+    REAL nroot_gauss(REAL e) const;
 };
 
 #endif
\ No newline at end of file