Refactor: remove PARAM usage from hsolver.cpp and two redundant reads

Group A (redundant reads -- the value was already available locally):

- hsolver_pw.cpp: `PARAM.globalv.use_uspp` -> `this->use_uspp`. HSolverPW
  already receives use_uspp through its constructor (injected from
  PARAM.globalv.use_uspp at the esolver level), so the direct read was
  duplicating an existing member.
- hsolver_lcao.cpp: `PARAM.inp.nspin` -> `nspin`. HSolverLCAO::solve()
  already takes nspin as an argument; it is now threaded into
  parakSolve() instead of being re-read from PARAM.

Group B (hsolver.cpp, now completely PARAM-free):

set_diagethr_ks() and set_diagethr_sdft() were already pure parameter-based
functions except for two leaked reads of PARAM.inp.scf_thr / PARAM.inp.nelec
in their nscf branches.

- set_diagethr_ks(): added `scf_thr_in`. The nscf branch used
  PARAM.inp.nelec even though the function already had a `nelec_in`
  parameter; both are now `nelec_in`. This is behaviour-preserving because
  the only call site passes PARAM.inp.nelec for that argument.
- set_diagethr_sdft(): added `nelec_in` and `scf_thr_in`.
- Dropped the now-unused parameter.h include and added the <algorithm> /
  <cmath> includes it was relying on transitively.

Call sites in source_esolver/esolver_ks.cpp updated accordingly; these are
the only external callers of the two functions.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: Critsium-xy

source/source_esolver/esolver_ks.cpp

@@ -184,13 +184,14 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter)
     {
         diag_ethr = hsolver::set_diagethr_ks(PARAM.inp.basis_type, PARAM.inp.esolver_type,
           PARAM.inp.calculation, PARAM.inp.init_chg, PARAM.inp.precision, istep, iter,
-          drho, PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nelec);
+          drho, PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nelec, PARAM.inp.scf_thr);
     }
     else if (PARAM.inp.esolver_type == "sdft")
     {
         diag_ethr = hsolver::set_diagethr_sdft(PARAM.inp.basis_type, PARAM.inp.esolver_type,
           PARAM.inp.calculation, PARAM.inp.init_chg, istep, iter, drho,
-          PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nbands, esolver_KS_ne);
+          PARAM.inp.pw_diag_thr, diag_ethr, PARAM.inp.nbands, esolver_KS_ne,
+          PARAM.inp.nelec, PARAM.inp.scf_thr);
     }
 
     // save input charge density (rho)

source/source_hsolver/hsolver.cpp

@@ -1,7 +1,9 @@
 #include "hsolver.h"
 
 #include "source_base/global_function.h"
-#include "source_io/module_parameter/parameter.h"
+
+#include <algorithm>
+#include <cmath>
 
 namespace hsolver
 {
@@ -16,7 +18,8 @@ double set_diagethr_ks(const std::string basis_type,
                        const double drho,
                        const double pw_diag_thr_init,
                        const double diag_ethr_in,
-                       const double nelec_in)
+                       const double nelec_in,
+                       const double scf_thr_in)
 {
     double res_diag_ethr = diag_ethr_in;
 
@@ -27,7 +30,7 @@ double set_diagethr_ks(const std::string basis_type,
         {
             if (res_diag_ethr - 1e-2 > -1e-5)
             {
-                res_diag_ethr = std::max(1e-13, 0.1 * std::min(1e-2, PARAM.inp.scf_thr / PARAM.inp.nelec));
+                res_diag_ethr = std::max(1e-13, 0.1 * std::min(1e-2, scf_thr_in / nelec_in));
             }
         }
         else if (iter == 1)
@@ -97,15 +100,17 @@ double set_diagethr_sdft(const std::string basis_type,
                          const double pw_diag_thr_init,
                          const double diag_ethr_in,
                          const int nband_in,
-                         const double stoiter_ks_ne_in)
+                         const double stoiter_ks_ne_in,
+                         const double nelec_in,
+                         const double scf_thr_in)
 {
     double res_diag_ethr = diag_ethr_in;
 
     if (basis_type == "pw" && esolver_type == "sdft")
     {
         if (calculation_in == "nscf")
         {
-            res_diag_ethr = std::max(std::min(1e-5, 0.1 * PARAM.inp.scf_thr / std::max(1.0, PARAM.inp.nelec)), 1e-12);
+            res_diag_ethr = std::max(std::min(1e-5, 0.1 * scf_thr_in / std::max(1.0, nelec_in)), 1e-12);
         }
         else if (iter == 1)
         {
@@ -124,7 +129,7 @@ double set_diagethr_sdft(const std::string basis_type,
         }
         else
         {
-            if (nband_in > 0 && stoiter_ks_ne_in > 1e-6) //PARAM.inp.nbands > 0 && this->stoiter.KS_ne > 1e-6
+            if (nband_in > 0 && stoiter_ks_ne_in > 1e-6)
             {
                 res_diag_ethr = std::min(res_diag_ethr, 0.1 * drho / std::max(1.0, stoiter_ks_ne_in));
             }

source/source_hsolver/hsolver.h

@@ -17,7 +17,8 @@ double set_diagethr_ks(const std::string basis_type,
                        const double drho,
                        const double pw_diag_thr_init,
                        const double diag_ethr_in,
-                       const double nelec_in);
+                       const double nelec_in,
+                       const double scf_thr_in);
 
 double set_diagethr_sdft(const std::string basis_type,
                          const std::string esolver_type,
@@ -29,7 +30,9 @@ double set_diagethr_sdft(const std::string basis_type,
                          const double pw_diag_thr_init,
                          const double diag_ethr_in,
                          const int nband_in,
-                         const double stoiter_ks_ne_in);
+                         const double stoiter_ks_ne_in,
+                         const double nelec_in,
+                         const double scf_thr_in);
 
 
 // reset diagethr according to drho and hsolver_error

source/source_hsolver/hsolver_lcao.cpp

@@ -62,7 +62,7 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
         if (PARAM.globalv.kpar_lcao > 1
             && (this->method == "genelpa" || this->method == "elpa" || this->method == "scalapack_gvx" || this->method == "lapack"))
         {
-            this->parakSolve(pHamilt, psi, pes, PARAM.globalv.kpar_lcao);
+            this->parakSolve(pHamilt, psi, pes, PARAM.globalv.kpar_lcao, nspin);
         } else
     #endif
         if (PARAM.globalv.kpar_lcao == 1)
@@ -192,7 +192,8 @@ template <typename T, typename Device>
 void HSolverLCAO<T, Device>::parakSolve(hamilt::Hamilt<T>* pHamilt,
                                         psi::Psi<T>& psi,
                                         elecstate::ElecState* pes,
-                                        int kpar)
+                                        const int kpar,
+                                        const int nspin)
 {
 #ifdef __MPI
     ModuleBase::timer::start("HSolverLCAO", "parakSolve");
@@ -202,7 +203,7 @@ void HSolverLCAO<T, Device>::parakSolve(hamilt::Hamilt<T>* pHamilt,
     int nks = psi.get_nk();
     int nrow = this->ParaV->get_global_row_size();
     int nb2d = this->ParaV->get_block_size();
-    k2d.set_para_env(psi.get_nk(), nrow, nb2d, GlobalV::NPROC, GlobalV::MY_RANK, PARAM.inp.nspin);
+    k2d.set_para_env(psi.get_nk(), nrow, nb2d, GlobalV::NPROC, GlobalV::MY_RANK, nspin);
     /// set psi_pool
     const int zero = 0;
     int coord_col = k2d.get_p2D_pool()->get_coord_col();

source/source_hsolver/hsolver_lcao.h

@@ -28,7 +28,11 @@ class HSolverLCAO
   private:
     void hamiltSolvePsiK(hamilt::Hamilt<TK>* hm, psi::Psi<TK>& psi, double* eigenvalue); // for kpar_lcao == 1
 
-    void parakSolve(hamilt::Hamilt<TK>* pHamilt, psi::Psi<TK>& psi, elecstate::ElecState* pes, int kpar); // for kpar_lcao > 1
+    void parakSolve(hamilt::Hamilt<TK>* pHamilt,
+                    psi::Psi<TK>& psi,
+                    elecstate::ElecState* pes,
+                    const int kpar,
+                    const int nspin); // for kpar_lcao > 1
 
     // The solving algorithm using cusolver is different from others, so a separate function is needed
     void parakSolve_cusolver(hamilt::Hamilt<TK>* pHamilt,

source/source_hsolver/hsolver_pw.cpp

@@ -212,7 +212,7 @@ void HSolverPW<T, Device>::solve(hamilt::Hamilt<T, Device>* pHamilt,
     elecstate::calEBand(_pes_pw->ekb,_pes_pw->wg,_pes_pw->f_en);
     if (skip_charge)
     {
-        if (PARAM.globalv.use_uspp)
+        if (this->use_uspp)
         {
             reinterpret_cast<elecstate::ElecStatePW<T, Device>*>(pes)->cal_becsum(psi);
         }