Try removing template dependency of ESolver (#7014)

* small format changes

* refactor(esolver): extract charge density symmetrization to Symmetry_rho::symmetrize_rho

- Add static method symmetrize_rho() in Symmetry_rho class
- Replace 7 duplicate code blocks with single function call
- Simplify code from 35 lines to 7 lines (80% reduction)
- Improve code readability and maintainability

Modified files:
- source_estate/module_charge/symmetry_rho.h: add static method declaration
- source_estate/module_charge/symmetry_rho.cpp: implement static method
- source_esolver/esolver_ks_lcao.cpp: 2 calls updated
- source_esolver/esolver_ks_pw.cpp: 1 call updated
- source_esolver/esolver_ks_lcao_tddft.cpp: 1 call updated
- source_esolver/esolver_ks_lcaopw.cpp: 1 call updated
- source_esolver/esolver_of.cpp: 1 call updated
- source_esolver/esolver_sdft_pw.cpp: 1 call updated

This refactoring follows the ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

* refactor(esolver): extract DeltaSpin lambda loop to deltaspin_lcao module

- Create new files deltaspin_lcao.h/cpp in module_deltaspin
- Extract DeltaSpin lambda loop logic from ESolver_KS_LCAO
- Simplify code from 18 lines to 1 line in hamilt2rho_single
- Separate LCAO and PW implementations for DeltaSpin

Modified files:
- source_esolver/esolver_ks_lcao.cpp: replace inline code with function call
- source_lcao/module_deltaspin/CMakeLists.txt: add new source file

New files:
- source_lcao/module_deltaspin/deltaspin_lcao.h: function declaration
- source_lcao/module_deltaspin/deltaspin_lcao.cpp: function implementation

This refactoring follows the ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

* refactor(esolver): complete DeltaSpin refactoring in LCAO

- Add init_deltaspin_lcao() function for DeltaSpin initialization
- Add cal_mi_lcao_wrapper() function for magnetic moment calculation
- Refactor all DeltaSpin-related code in esolver_ks_lcao.cpp
- Simplify code from 29 lines to 3 lines (90% reduction)

Modified files:
- source_esolver/esolver_ks_lcao.cpp: replace 3 code blocks with function calls
- source_lcao/module_deltaspin/deltaspin_lcao.h: add 2 new function declarations
- source_lcao/module_deltaspin/deltaspin_lcao.cpp: implement 2 new functions

This completes the DeltaSpin refactoring for LCAO method:
1. init_deltaspin_lcao() - initialize DeltaSpin calculation
2. cal_mi_lcao_wrapper() - calculate magnetic moments
3. run_deltaspin_lambda_loop_lcao() - run lambda loop optimization

All functions follow the ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

* refactor(esolver): extract DFT+U code to dftu_lcao module

- Create new files dftu_lcao.h/cpp in source_lcao directory
- Add init_dftu_lcao() function for DFT+U initialization
- Add finish_dftu_lcao() function for DFT+U finalization
- Simplify code from 32 lines to 2 lines in esolver_ks_lcao.cpp
- Remove conditional checks from ESolver, move them to functions

Modified files:
- source_esolver/esolver_ks_lcao.cpp: replace 2 code blocks with function calls
- source_lcao/CMakeLists.txt: add new source file

New files:
- source_lcao/dftu_lcao.h: function declarations
- source_lcao/dftu_lcao.cpp: function implementations

This refactoring prepares for unifying old and new DFT+U implementations:
- Old DFT+U: source_lcao/module_dftu/
- New DFT+U: source_lcao/module_operator_lcao/op_dftu_lcao.cpp

All functions follow ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

* refactor(esolver): extract diagonalization parameters setup to hsolver module

- Create new files diago_params.h/cpp in source_hsolver directory
- Add setup_diago_params_pw() function for PW diagonalization parameters
- Simplify code from 11 lines to 1 line in esolver_ks_pw.cpp
- Encapsulate diagonalization parameter setup logic

Modified files:
- source_esolver/esolver_ks_pw.cpp: replace inline code with function call
- source_hsolver/CMakeLists.txt: add new source file

New files:
- source_hsolver/diago_params.h: function declaration
- source_hsolver/diago_params.cpp: function implementation

This refactoring follows ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

* fix(deltaspin): add sc_mag_switch check in cal_mi_lcao_wrapper

- Add Input_para parameter to cal_mi_lcao_wrapper function
- Add sc_mag_switch check to avoid calling cal_mi_lcao when DeltaSpin is disabled
- Fix 'atomCounts is not set' error in non-DeltaSpin calculations
- Update function call in esolver_ks_lcao.cpp

This fix resolves the CI/CD failure caused by commit 2a520e3.
The root cause was that cal_mi_lcao_wrapper was called without checking
sc_mag_switch, leading to uninitialized atomCounts error.

Modified files:
- source_esolver/esolver_ks_lcao.cpp: update function call
- source_lcao/module_deltaspin/deltaspin_lcao.h: add parameter
- source_lcao/module_deltaspin/deltaspin_lcao.cpp: add check

This follows the refactoring principle: preserve original condition checks
when extracting code to wrapper functions.

* fix(deltaspin): add #ifdef __LCAO for conditional compilation

- Add #ifdef __LCAO conditional compilation in init_deltaspin_lcao and cal_mi_lcao_wrapper
- Fix parameter order in init_sc call for LCAO and non-LCAO builds
- Fix undefined reference to cal_mi_lcao in non-LCAO build

This fix resolves CI/CD compilation errors in both build_5pt (with __LCAO) and build_1p (without __LCAO) environments.

The The root cause was 1. init_sc has different parameter order in LCAO vs non-LCAO builds
   - LCAO: psi, dm, pelec
   - non-LCAO: psi, pelec
  2. cal_mi_lcao is only defined in LCAO build

Modified files:
- source_hsolver/diago_params.h: add setup_diago_params_sdft declaration
- source_lcao/module_deltaspin/deltaspin_lcao.cpp: add conditional compilation

This follows the refactoring principle: handle conditional compilation properly
when code has different implementations for different build configurations.

* refactor(esolver): extract SDFT diagonalization parameters setup

- Add setup_diago_params_sdft() function for SDFT diagonalization parameters
- Simplify code from 11 lines to 1 line in esolver_sdft_pw.cpp
- Encapsulate diagonalization parameter setup logic for SDFT

Modified files:
- source_esolver/esolver_sdft_pw.cpp: replace inline code with function call
- source_hsolver/diago_params.cpp: add setup_diago_params_sdft implementation

This refactoring follows ESolver cleanup principle:
keep ESolver focused on high-level workflow control.

Note: SDFT has different parameter setup logic compared to PW:
- Different need_subspace condition
- No SCF_ITER setting
- Always set PW_DIAG_NMAX (no nscf check)

* refactor(hamilt): introduce HamiltBase non-template base class

- Create HamiltBase as a non-template base class for Hamilt<T, Device>
- Modify Hamilt<T, Device> to inherit from HamiltBase
- Change ESolver_KS::p_hamilt type from Hamilt<T, Device>* to HamiltBase*
- Add static_cast where needed when passing p_hamilt to functions expecting Hamilt<T, Device>*

This is the first step towards removing template parameters from ESolver.

Modified files:
- source/source_esolver/esolver_ks.h
- source/source_esolver/esolver_ks_lcaopw.cpp
- source/source_esolver/esolver_ks_pw.cpp
- source/source_esolver/esolver_sdft_pw.cpp
- source/source_hamilt/hamilt.h

New files:
- source/source_hamilt/hamilt_base.h

* refactor(esolver): add static_cast for p_hamilt in esolver files

- Add static_cast<hamilt::Hamilt<T>*> when passing p_hamilt to functions
  expecting Hamilt<T, Device>* type
- Split long cast statements into multiple lines for better readability
- Files modified:
  - esolver_ks_pw.cpp: setup_pot, stp.init calls
  - esolver_ks_lcao.cpp: init_chg_hr, hsolver_lcao_obj.solve calls
  - esolver_ks_lcao_tddft.cpp: solve_psi, cal_edm_tddft, matrix calls
  - esolver_gets.cpp: ops access, output_SR call

This follows the HamiltBase refactoring strategy where p_hamilt is stored
as HamiltBase* and cast to Hamilt<T, Device>* when needed.

---------

Co-authored-by: abacus_fixer <mohanchen@pku.eud.cn>

Author: mohanchen

source/source_esolver/esolver_gets.cpp

@@ -108,8 +108,9 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
                                                                                      this->kv,
                                                                                      *(two_center_bundle_.overlap_orb),
                                                                                      orb_.cutoffs());
-            dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>*>(this->p_hamilt->ops)
-                ->contributeHR();
+            auto* hamilt_ptr = static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt);
+            auto* ops_ptr = dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>*>(hamilt_ptr->ops);
+            ops_ptr->contributeHR();
         }
         else
         {
@@ -119,13 +120,16 @@ void ESolver_GetS::runner(UnitCell& ucell, const int istep)
                                                                                   this->kv,
                                                                                   *(two_center_bundle_.overlap_orb),
                                                                                   orb_.cutoffs());
-            dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(this->p_hamilt->ops)->contributeHR();
+            auto* hamilt_ptr = static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt);
+            auto* ops_ptr = dynamic_cast<hamilt::OperatorLCAO<std::complex<double>, double>*>(hamilt_ptr->ops);
+            ops_ptr->contributeHR();
         }
     }
 
     const std::string fn = PARAM.globalv.global_out_dir + "sr_nao.csr";
 
-    ModuleIO::output_SR(pv, gd, this->p_hamilt, fn);
+    auto* hamilt_ptr = static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt);
+    ModuleIO::output_SR(pv, gd, hamilt_ptr, fn);
 
     if (PARAM.inp.out_mat_r)
     {

source/source_esolver/esolver_ks.h

@@ -7,6 +7,7 @@
 #include "source_estate/module_charge/charge_mixing.h" // use charge mixing
 #include "source_psi/psi.h" // use electronic wave functions
 #include "source_hamilt/hamilt.h" // use Hamiltonian
+#include "source_hamilt/hamilt_base.h" // use Hamiltonian base class
 #include "source_lcao/module_dftu/dftu.h" // mohan add 20251107
 
 namespace ModuleESolver
@@ -47,8 +48,8 @@ class ESolver_KS : public ESolver_FP
     //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
     virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override;
 
-    //! Hamiltonian
-    hamilt::Hamilt<T, Device>* p_hamilt = nullptr;
+    //! Hamiltonian (base class pointer, actual type determined at runtime)
+    hamilt::HamiltBase* p_hamilt = nullptr;
 
     //! PW for wave functions, only used in KSDFT, not in OFDFT
     ModulePW::PW_Basis_K* pw_wfc = nullptr;

source/source_esolver/esolver_ks_lcao.cpp

@@ -179,7 +179,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         {
             //! 13.1.2) init charge density from Hamiltonian matrix file
             LCAO_domain::init_chg_hr<TK, TR>(PARAM.globalv.global_readin_dir, PARAM.inp.nspin,
-                this->p_hamilt, ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm,
+                static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), ucell, &(this->pv), this->psi[0], this->pelec, *this->dmat.dm,
                 this->chr, PARAM.inp.ks_solver);
         }
     }
@@ -382,7 +382,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
     if (!skip_solve)
     {
         hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, *this->dmat.dm, 
+        hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, 
           this->chr, PARAM.inp.nspin, skip_charge);
     }
 

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -235,7 +235,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
                                                       PARAM.inp.nbands,
                                                       PARAM.globalv.nlocal,
                                                       this->kv.get_nks(),
-                                                      this->p_hamilt,
+                                                      static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                                       this->pv,
                                                       this->psi,
                                                       this->psi_laststep,
@@ -255,7 +255,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
                                                   PARAM.inp.nbands,
                                                   PARAM.globalv.nlocal,
                                                   this->kv.get_nks(),
-                                                  this->p_hamilt,
+                                                  static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                                   this->pv,
                                                   this->psi,
                                                   this->psi_laststep,
@@ -277,7 +277,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
         {
             bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
             hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver);
-            hsolver_lcao_obj.solve(this->p_hamilt,
+            hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                    this->psi[0],
                                    this->pelec,
                                    *this->dmat.dm,
@@ -342,11 +342,11 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::iter_finish(UnitCell& ucell,
     {
         if (use_tensor && use_lapack)
         {
-            elecstate::cal_edm_tddft_tensor_lapack<Device>(this->pv, this->dmat, this->kv, this->p_hamilt);
+            elecstate::cal_edm_tddft_tensor_lapack<Device>(this->pv, this->dmat, this->kv, static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt));
         }
         else
         {
-            elecstate::cal_edm_tddft(this->pv, this->dmat, this->kv, this->p_hamilt);
+            elecstate::cal_edm_tddft(this->pv, this->dmat, this->kv, static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt));
         }
     }
 }
@@ -416,7 +416,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::store_h_s_psi(UnitCell& ucell,
             this->p_hamilt->updateHk(ik);
             hamilt::MatrixBlock<std::complex<double>> h_mat;
             hamilt::MatrixBlock<std::complex<double>> s_mat;
-            this->p_hamilt->matrix(h_mat, s_mat);
+            static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt)->matrix(h_mat, s_mat);
 
             // Store H and S matrices to Hk_laststep and Sk_laststep
             if (use_tensor && use_lapack)

source/source_esolver/esolver_ks_lcaopw.cpp

@@ -146,7 +146,7 @@ namespace ModuleESolver
         bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
         hsolver::HSolverLIP<T> hsolver_lip_obj(this->pw_wfc);
-        hsolver_lip_obj.solve(this->p_hamilt, this->stp.psi_t[0], this->pelec, 
+        hsolver_lip_obj.solve(static_cast<hamilt::Hamilt<T>*>(this->p_hamilt), this->stp.psi_t[0], this->pelec, 
           *this->psi_local, skip_charge,ucell.tpiba,ucell.nat);
 
         // add exx

source/source_esolver/esolver_ks_pw.cpp

@@ -128,10 +128,10 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     // init DFT+U is done in "before_all_runners" in LCAO basis. This should be refactored, mohan note 2025-11-06
     pw::setup_pot(istep, ucell, this->kv, this->sf, this->pelec, this->Pgrid,
               this->chr, this->locpp, this->ppcell, this->dftu, this->vsep_cell,
-              this->stp.psi_t, this->p_hamilt, this->pw_wfc, this->pw_rhod, PARAM.inp);
+              this->stp.psi_t, static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt), this->pw_wfc, this->pw_rhod, PARAM.inp);
 
     // setup psi (electronic wave functions)
-    this->stp.init(this->p_hamilt);
+    this->stp.init(static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt));
 
     //! Setup EXX helper for Hamiltonian and psi
     exx_helper.before_scf(this->p_hamilt, this->stp.psi_t, PARAM.inp);
@@ -188,7 +188,7 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
                                                      hsolver::DiagoIterAssist<T, Device>::need_subspace,
                                                      PARAM.inp.use_k_continuity);
 
-        hsolver_pw_obj.solve(this->p_hamilt, this->stp.psi_t[0], this->pelec, this->pelec->ekb.c,
+        hsolver_pw_obj.solve(static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt), this->stp.psi_t[0], this->pelec, this->pelec->ekb.c,
           GlobalV::RANK_IN_POOL, GlobalV::NPROC_IN_POOL, skip_charge, ucell.tpiba, ucell.nat);
     }
 

source/source_esolver/esolver_sdft_pw.cpp

@@ -167,7 +167,7 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
                                                            hsolver::DiagoIterAssist<T, Device>::need_subspace);
 
     hsolver_pw_sdft_obj.solve(ucell,
-                              this->p_hamilt,
+                              static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt),
                               this->stp.psi_t[0],
                               this->stp.psi_cpu[0],
                               this->pelec,
@@ -291,7 +291,7 @@ void ESolver_SDFT_PW<T, Device>::after_all_runners(UnitCell& ucell)
                                               this->pw_wfc,
                                               this->stp.psi_t,
                                               &this->ppcell,
-                                              this->p_hamilt,
+                                              static_cast<hamilt::Hamilt<std::complex<double>, Device>*>(this->p_hamilt),
                                               this->stoche,
                                               &stowf);
         sto_elecond.decide_nche(PARAM.inp.cond_dt, 1e-8, this->nche_sto, PARAM.inp.emin_sto, PARAM.inp.emax_sto);

source/source_hamilt/hamilt.h

@@ -7,21 +7,28 @@
 #include "matrixblock.h"
 #include "source_psi/psi.h"
 #include "operator.h"
+#include "hamilt_base.h"
 
 namespace hamilt
 {
 
 template <typename T, typename Device = base_device::DEVICE_CPU>
-class Hamilt
+class Hamilt : public HamiltBase
 {
   public:
     virtual ~Hamilt(){};
 
     /// for target K point, update consequence of hPsi() and matrix()
-    virtual void updateHk(const int ik){return;}
+    void updateHk(const int ik) override { return; }
 
     /// refresh status of Hamiltonian, for example, refresh H(R) and S(R) in LCAO case
-    virtual void refresh(bool yes = true){return;}
+    void refresh(bool yes = true) override { return; }
+
+    /// get the class name
+    std::string get_classname() const override { return classname; }
+
+    /// get the operator chain
+    void* get_ops() override { return static_cast<void*>(ops); }
 
     /// core function: for solving eigenvalues of Hamiltonian with iterative method
 	virtual void hPsi(

source/source_hamilt/hamilt_base.h

@@ -0,0 +1,52 @@
+#ifndef HAMILT_BASE_H
+#define HAMILT_BASE_H
+
+#include <string>
+
+namespace hamilt
+{
+
+/**
+ * @brief Base class for Hamiltonian
+ *
+ * This is a non-template base class for Hamilt<T, Device>.
+ * It provides a common interface for all Hamiltonian types,
+ * allowing ESolver to manage Hamiltonian without template parameters.
+ */
+class HamiltBase
+{
+  public:
+    virtual ~HamiltBase() {}
+
+    /**
+     * @brief Update Hamiltonian for a specific k-point
+     *
+     * @param ik k-point index
+     */
+    virtual void updateHk(const int ik) { return; }
+
+    /**
+     * @brief Refresh the status of Hamiltonian
+     *
+     * @param yes whether to refresh
+     */
+    virtual void refresh(bool yes = true) { return; }
+
+    /**
+     * @brief Get the class name
+     *
+     * @return class name
+     */
+    virtual std::string get_classname() const { return "none"; }
+
+    /**
+     * @brief Get the operator chain (as void* to avoid template)
+     *
+     * @return pointer to operator chain
+     */
+    virtual void* get_ops() { return nullptr; }
+};
+
+} // namespace hamilt
+
+#endif // HAMILT_BASE_H