fix(deltaspin): reuse pre-initialized phsol in lambda loop for solver convergence

Previously HSolverLCAO was created on stack each inner step in
cal_mw_from_lambda, losing solver state between calls. Now phsol is
initialized once in ESolver_KS_LCAO::before_scf and passed through
init_sc -> set_solver_parameters -> SpinConstrain::phsol, matching
zdy-tmp behavior.

Changes:
- ESolver_KS_LCAO::phsol member: persistent solver pointer
- before_scf: create phsol once, pass to init_sc
- hamilt2rho_single: reuse phsol instead of stack solver
- init_sc/set_solver_parameters: accept phsol_in parameter
- cal_mw_from_lambda: use this->phsol when available

Author: dyzheng

source/source_esolver/esolver_ks_lcao.cpp

@@ -153,13 +153,19 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     // since it depends on ionic positions
     this->deepks.build_overlap(ucell, orb_, pv, gd, *(two_center_bundle_.overlap_orb_alpha), PARAM.inp);
 
-    // 10) prepare sc calculation
+    // 10) initialize HSolver once (reuse across SCF iterations and lambda loop)
+    if (this->phsol == nullptr)
+    {
+        this->phsol = new hsolver::HSolverLCAO<TK>(&(this->pv), PARAM.inp.ks_solver);
+    }
+
+    // 11) prepare sc calculation
     if (PARAM.inp.sc_mag_switch)
     {
         spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
         sc.init_sc(PARAM.inp.sc_thr, PARAM.inp.nsc, PARAM.inp.nsc_min, PARAM.inp.alpha_trial,
                    PARAM.inp.sccut, PARAM.inp.sc_drop_thr, ucell, &(this->pv),
-                   PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec);
+                   PARAM.inp.nspin, this->kv, this->p_hamilt, this->psi, this->dmat.dm, this->pelec, nullptr, this->phsol);
         // Set lambda update strategy
         if (PARAM.inp.sc_lambda_strategy == "linear_response")
         {
@@ -460,8 +466,9 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2rho_single(UnitCell& ucell, int istep, int
     // 3) run Hsolver
     if (!skip_solve)
     {
-        hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-        hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), this->psi[0], this->pelec, *this->dmat.dm, 
+        hsolver::HSolverLCAO<TK>* hsolver_lcao_obj
+            = static_cast<hsolver::HSolverLCAO<TK>*>(this->phsol);
+        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.h

@@ -96,6 +96,9 @@ class ESolver_KS_LCAO : public ESolver_KS
     friend class LR::ESolver_LR<double, double>;
     friend class LR::ESolver_LR<std::complex<double>, double>;
 
+    //! Persistent solver pointer (for SpinConstrain reuse)
+    void* phsol = nullptr;
+
     // Temporarily store the stress to unify the interface with PW,
     // because it's hard to seperate force and stress calculation in LCAO.
     ModuleBase::matrix scs;

source/source_lcao/module_deltaspin/CMakeLists.txt

@@ -12,6 +12,7 @@ list(APPEND objects
     template_helpers.cpp
     lambda_update_strategies.cpp
     lambda_strategy_integration.cpp
+    lambda_solvers.cpp
     deltaspin_lcao.cpp
 )
 

source/source_lcao/module_deltaspin/cal_mw_from_lambda.cpp

@@ -122,7 +122,6 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int
     {
         psi::Psi<std::complex<double>>* psi_t = static_cast<psi::Psi<std::complex<double>>*>(this->psi);
         hamilt::Hamilt<std::complex<double>>* hamilt_t = static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt);
-        hsolver::HSolverLCAO<std::complex<double>> hsolver_t(this->ParaV, PARAM.inp.ks_solver);
         if (PARAM.inp.nspin == 2)
         {
             dynamic_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, double>>*>(this->p_operator)->update_lambda();
@@ -131,7 +130,19 @@ void spinconstrain::SpinConstrain<std::complex<double>>::cal_mw_from_lambda(int
         {
             dynamic_cast<hamilt::DeltaSpin<hamilt::OperatorLCAO<std::complex<double>, std::complex<double>>>*>(this->p_operator)->update_lambda();
         }
-        hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, PARAM.inp.nspin, true);
+        // Use the pre-initialized solver from ESolver (same as zdy-tmp)
+        if (this->phsol != nullptr)
+        {
+            hsolver::HSolverLCAO<std::complex<double>>* hsolver_t
+                = static_cast<hsolver::HSolverLCAO<std::complex<double>>*>(this->phsol);
+            hsolver_t->solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, PARAM.inp.nspin, true);
+        }
+        else
+        {
+            // Fallback: create solver on stack (less efficient but functional)
+            hsolver::HSolverLCAO<std::complex<double>> hsolver_t(this->ParaV, PARAM.inp.ks_solver);
+            hsolver_t.solve(hamilt_t, psi_t[0], this->pelec, *this->dm_, *this->pelec->charge, PARAM.inp.nspin, true);
+        }
         elecstate::calculate_weights(this->pelec->ekb, this->pelec->wg, this->pelec->klist,
                                      this->pelec->eferm, this->pelec->f_en, this->pelec->nelec_spin, this->pelec->skip_weights);
         elecstate::calEBand(this->pelec->ekb, this->pelec->wg, this->pelec->f_en);

source/source_lcao/module_deltaspin/init_sc.cpp

@@ -18,7 +18,8 @@ void spinconstrain::SpinConstrain<TK>::init_sc(double sc_thr_in,
 		elecstate::DensityMatrix<TK, double>* dm_in, // mohan add 2025-11-03
 #endif
 		elecstate::ElecState* pelec_in,
-		ModulePW::PW_Basis_K* pw_wfc_in)
+		ModulePW::PW_Basis_K* pw_wfc_in,
+		void* phsol_in)
 {
     this->set_input_parameters(sc_thr_in, nsc_in, nsc_min_in, alpha_trial_in, sccut_in, sc_drop_thr_in);
     this->set_atomCounts(ucell.get_atom_Counts());
@@ -33,7 +34,7 @@ void spinconstrain::SpinConstrain<TK>::init_sc(double sc_thr_in,
     this->pw_wfc_ = pw_wfc_in;
     this->set_decay_grad();
     if(ParaV_in != nullptr) this->set_ParaV(ParaV_in);
-    this->set_solver_parameters(kv_in, p_hamilt_in, psi_in, pelec_in);
+    this->set_solver_parameters(kv_in, p_hamilt_in, psi_in, pelec_in, phsol_in);
 #ifdef __LCAO
     this->dm_ = dm_in; // mohan add 2025-11-03
 #endif