Feature: Add LCAO subspace solver for SCF acceleration in MD/relax

- Add HSolverLCAOSubspace class in hsolver module for subspace diagonalization
- Integrate subspace solver into ESolver_KS_LCAO for MD/relax acceleration
- Add input parameters: lcao_subspace_persistent, lcao_subspace_clear_thr
- First SCF step uses full diagonalization, subsequent steps use subspace solver
- Subspace cache updated after each SCF iteration using current wavefunctions
- Add cal_PI_sub method to dspin_lcao for projector overlap computation
- Add LambdaSolver abstraction with BFGS, ChiGuided, Subspace, FDCG solvers

Author: dyzheng

source/source_esolver/esolver_ks_lcao.cpp

@@ -212,6 +212,40 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         rdmft_solver.update_ion(ucell, *(this->pw_rho), this->locpp.vloc, this->sf.strucFac);
     }
 
+    // 18) Manage subspace solver cache for MD/relax acceleration
+    if (PARAM.inp.lcao_subspace_persistent && istep > 0)
+    {
+        // Check atomic displacement
+        double max_disp = 0.0;
+        for (int iat = 0; iat < ucell.nat; iat++)
+        {
+            const auto& tau = ucell.get_tau(iat);
+            double dx = tau.x - this->last_atom_positions_[iat].x;
+            double dy = tau.y - this->last_atom_positions_[iat].y;
+            double dz = tau.z - this->last_atom_positions_[iat].z;
+            double disp = std::sqrt(dx*dx + dy*dy + dz*dz) * ucell.lat0;
+            if (disp > max_disp) max_disp = disp;
+        }
+
+        // Clear cache if displacement exceeds threshold or cache is invalid
+        if (this->subspace_solver_ && 
+            (PARAM.inp.lcao_subspace_clear_thr > 0.0 && 
+             max_disp > PARAM.inp.lcao_subspace_clear_thr))
+        {
+            this->subspace_solver_->clear_subspace();
+            GlobalV::ofs_running << " >> Subspace cache cleared: max atomic displacement = " 
+                                 << max_disp << " Bohr (threshold: " 
+                                 << PARAM.inp.lcao_subspace_clear_thr << " Bohr)" << std::endl;
+        }
+    }
+
+    // Save current atom positions for next step comparison
+    this->last_atom_positions_.resize(ucell.nat);
+    for (int iat = 0; iat < ucell.nat; iat++)
+    {
+        this->last_atom_positions_[iat] = ucell.get_tau(iat);
+    }
+
     ModuleBase::timer::end("ESolver_KS_LCAO", "before_scf");
     return;
 }
@@ -404,9 +438,31 @@ 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, 
-          this->chr, PARAM.inp.nspin, skip_charge);
+        // Try subspace solver first (for SCF acceleration after first step)
+        bool subspace_used = false;
+        if constexpr (std::is_same_v<TK, std::complex<double>>)
+        {
+            if (this->subspace_solver_ && this->subspace_solver_->has_subspace())
+            {
+                this->subspace_solver_->solve(
+                    static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt), 
+                    this->psi[0], 
+                    this->pelec, 
+                    *this->dmat.dm,
+                    this->chr, 
+                    PARAM.inp.nspin, 
+                    skip_charge);
+                subspace_used = true;
+            }
+        }
+
+        // Fall back to standard HSolverLCAO if subspace not available
+        if (!subspace_used)
+        {
+            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, 
+              this->chr, PARAM.inp.nspin, skip_charge);
+        }
     }
 
     // 4) EXX
@@ -484,6 +540,32 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
       this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 
       hamilt_lcao, this->orb_, this->deepks, 
       this->exx_nao, iter, istep, conv_esolver, this->scf_ene_thr);
+
+    // Update subspace cache after each SCF iteration (for next iteration acceleration)
+    // First SCF step (istep=0, iter=1): full diagonalization already done by HSolverLCAO
+    //   -> update_subspace_cache() uses the converged wavefunctions to build initial cache
+    // Subsequent SCF steps: subspace solver used, then cache updated with new wavefunctions
+    if constexpr (std::is_same_v<TK, std::complex<double>>)
+    {
+        if (PARAM.inp.lcao_subspace_persistent)
+        {
+            if (!this->subspace_solver_)
+            {
+                this->subspace_solver_.reset(
+                    new hsolver::HSolverLCAOSubspace(&(this->pv), PARAM.inp.ks_solver));
+                this->subspace_solver_->set_persistent(true);
+            }
+
+            // Update subspace cache using current wavefunctions
+            // lambda_ref is zero because we don't use DeltaSpin perturbation here
+            std::vector<ModuleBase::Vector3<double>> lambda_ref(ucell.nat, {0.0, 0.0, 0.0});
+            this->subspace_solver_->update_subspace_cache(
+                static_cast<hamilt::Hamilt<TK>*>(this->p_hamilt),
+                this->psi[0],
+                this->pelec,
+                lambda_ref);
+        }
+    }
 }
 
 template <typename TK, typename TR>

source/source_esolver/esolver_ks_lcao.h

@@ -11,6 +11,7 @@
 #include "source_lcao/setup_exx.h" // for exx, mohan add 20251008
 #include "source_lcao/module_rdmft/rdmft.h" // rdmft
 #include "source_lcao/setup_dm.h" // mohan add 2025-10-30
+#include "source_hsolver/hsolver_lcao_subspace.h" // subspace solver for SCF acceleration
 
 #include <memory>
 
@@ -92,6 +93,10 @@ class ESolver_KS_LCAO : public ESolver_KS
     //! For RDMFT calculations, added by jghan, 2024-03-16 
     rdmft::RDMFT<TK, TR> rdmft_solver;
 
+    //! Subspace solver for accelerating SCF in MD/relax ionic steps
+    std::unique_ptr<hsolver::HSolverLCAOSubspace> subspace_solver_;
+    std::vector<ModuleBase::Vector3<double>> last_atom_positions_;
+
     //! For linear-response TDDFT
     friend class LR::ESolver_LR<double, double>;
     friend class LR::ESolver_LR<std::complex<double>, double>;

source/source_hsolver/CMakeLists.txt

@@ -20,13 +20,15 @@ if(ENABLE_LCAO)
   if (ENABLE_MPI)
     list(APPEND objects
         hsolver_lcao.cpp
+        hsolver_lcao_subspace.cpp
         diago_scalapack.cpp
         parallel_k2d.cpp
         diago_lapack.cpp
     )
   else ()
     list(APPEND objects
         hsolver_lcao.cpp
+        hsolver_lcao_subspace.cpp
         diago_lapack.cpp
     )
   endif ()