Refactor: inject solver config into HSolver/PEXSI instead of reading PARAM

Group C of the source_hsolver PARAM removal: the parameters that were read
from the PARAM global mid-algorithm are now injected through constructors,
following the existing explicit-scalar style of the HSolverPW / HSolverLCAO
constructors.

HSolverPW gains nbands, diago_smooth_ethr, pw_diag_ndim, diag_subspace and
nb2d. Note that PARAM.inp.nbands (the global band count used by bpcg's
init_iter) is deliberately kept distinct from psi.get_nbands(), which is the
local count under band parallelism.

HSolverPW_SDFT gains ks_run, all_ks_run and bndpar, and forwards the five new
base-class values.

HSolverLIP gains use_uspp.

HSolverLCAO gains kpar_lcao, plus nlocal and nelec for the pexsi branch.

PEXSI chain: the two reads in module_pexsi/pexsi_solver.cpp could not simply
be lifted to its caller, because DiagoPexsi did not have nelec either, so the
whole chain is threaded in one go:

- PEXSI_Solver::prepare() takes nlocal and nelec and stores them, alongside
  the nb / nrow / ncol it already received.
- DiagoPexsi's constructor takes nspin, nlocal and nelec. The nspin == 4 ->
  single-density-matrix collapse was computed identically in the constructor
  and the destructor; it is now stored once as nspin_dm, which also removes
  the destructor's dependency on PARAM still holding the same nspin.
- HSolverLCAO passes its solve()-argument nspin plus the injected nlocal and
  nelec down to DiagoPexsi.

Two incidental fixes in code that had to be touched anyway:

- DiagoPexsi::diag() declared a std::vector<double> eigen(nlocal) that was
  never used; removed. This was the only nlocal read in that function.
- The constructor read pexsi_mu through this->ps before ps was assigned. It
  is a static member so this happened to work, but it is UB on a null
  unique_ptr; it now reads pexsi::PEXSI_Solver::pexsi_mu directly.

Also fixes a pre-existing bug in test_hsolver_pw.cpp: the HSolverPW fixture
passed an extra `false` after method_in, so every argument from nspin_in
onwards was shifted by one (nspin_in received use_uspp, diag_thr_in received
PW_DIAG_NMAX, need_subspace_in received PW_DIAG_THR, and so on). Adding the
new parameters made the call fail to compile, which surfaced it. The affected
tests only construct the object or exercise the early npw_total < nbands
guard, so no assertion depended on the shifted values.

PARAM occurrences in source_hsolver production code: 83 -> 60. Remaining are
the six dense LCAO diagonalizers sharing nlocal/nbands (group D) and
DiagoIterAssist's basis_type/calculation switches (group E).

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

Author: Critsium-xy

source/source_esolver/esolver_ks_lcao.cpp

@@ -434,7 +434,11 @@ 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::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv),
+                                                  PARAM.inp.ks_solver,
+                                                  PARAM.globalv.kpar_lcao,
+                                                  PARAM.globalv.nlocal,
+                                                  PARAM.inp.nelec);
         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

@@ -329,7 +329,11 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::hamilt2rho_single(UnitCell& ucell,
         if (this->psi != nullptr)
         {
             bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
-            hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv, PARAM.inp.ks_solver);
+            hsolver::HSolverLCAO<std::complex<double>> hsolver_lcao_obj(&this->pv,
+                                                                        PARAM.inp.ks_solver,
+                                                                        PARAM.globalv.kpar_lcao,
+                                                                        PARAM.globalv.nlocal,
+                                                                        PARAM.inp.nelec);
             hsolver_lcao_obj.solve(static_cast<hamilt::Hamilt<std::complex<double>>*>(this->p_hamilt),
                                    this->psi[0],
                                    this->pelec,

source/source_esolver/esolver_ks_lcaopw.cpp

@@ -133,7 +133,7 @@ namespace ModuleESolver
         hsolver::DiagoIterAssist<T>::PW_DIAG_NMAX = PARAM.inp.pw_diag_nmax;
         bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
-        hsolver::HSolverLIP<T> hsolver_lip_obj(this->pw_wfc);
+        hsolver::HSolverLIP<T> hsolver_lip_obj(this->pw_wfc, PARAM.globalv.use_uspp);
         hsolver_lip_obj.solve(static_cast<hamilt::Hamilt<T>*>(this->p_hamilt), *this->stp.template get_psi_t<T, base_device::DEVICE_CPU>(), this->pelec, 
           *this->psi_local, skip_charge,ucell.tpiba,ucell.nat);
 

source/source_esolver/esolver_ks_pw.cpp

@@ -222,6 +222,11 @@ void ESolver_KS_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, const int iste
                                                      hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX,
                                                      hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR,
                                                      hsolver::DiagoIterAssist<T, Device>::need_subspace,
+                                                     PARAM.inp.nbands,
+                                                     PARAM.inp.diago_smooth_ethr,
+                                                     PARAM.inp.pw_diag_ndim,
+                                                     PARAM.inp.diag_subspace,
+                                                     PARAM.inp.nb2d,
                                                      PARAM.inp.use_k_continuity);
 
         hsolver_pw_obj.solve(static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt), *this->stp.template get_psi_t<T, Device>(), this->pelec, this->pelec->ekb.c,

source/source_esolver/esolver_sdft_pw.cpp

@@ -163,7 +163,15 @@ void ESolver_SDFT_PW<T, Device>::hamilt2rho_single(UnitCell& ucell, int istep, i
                                                            hsolver::DiagoIterAssist<T, Device>::SCF_ITER,
                                                            hsolver::DiagoIterAssist<T, Device>::PW_DIAG_NMAX,
                                                            hsolver::DiagoIterAssist<T, Device>::PW_DIAG_THR,
-                                                           hsolver::DiagoIterAssist<T, Device>::need_subspace);
+                                                           hsolver::DiagoIterAssist<T, Device>::need_subspace,
+                                                           PARAM.inp.nbands,
+                                                           PARAM.inp.diago_smooth_ethr,
+                                                           PARAM.inp.pw_diag_ndim,
+                                                           PARAM.inp.diag_subspace,
+                                                           PARAM.inp.nb2d,
+                                                           PARAM.globalv.ks_run,
+                                                           PARAM.globalv.all_ks_run,
+                                                           PARAM.inp.bndpar);
 
     hsolver_pw_sdft_obj.solve(ucell,
                               static_cast<hamilt::Hamilt<T, Device>*>(this->p_hamilt),

source/source_hsolver/diago_pexsi.cpp

@@ -1,6 +1,5 @@
 #include <mpi.h>
 #include <complex>
-#include "source_io/module_parameter/parameter.h"
 #include <memory>
 #ifdef __PEXSI
 #include "diago_pexsi.h"
@@ -19,25 +18,27 @@ template <typename T>
 std::vector<double> DiagoPexsi<T>::mu_buffer;
 
 template <typename T>
-DiagoPexsi<T>::DiagoPexsi(const Parallel_Orbitals* ParaV_in)
+DiagoPexsi<T>::DiagoPexsi(const Parallel_Orbitals* ParaV_in,
+                          const int nspin_in,
+                          const int nlocal_in,
+                          const double nelec_in)
 {
-    int nspin = PARAM.inp.nspin;
-    if (PARAM.inp.nspin == 4)
-    {
-        nspin = 1;
-    }
-    mu_buffer.resize(nspin);
-    for (int i = 0; i < nspin; i++)
+    this->nspin_dm = (nspin_in == 4) ? 1 : nspin_in;
+    this->nlocal = nlocal_in;
+    this->nelec = nelec_in;
+
+    mu_buffer.resize(this->nspin_dm);
+    for (int i = 0; i < this->nspin_dm; i++)
     {
-        mu_buffer[i] = this->ps->pexsi_mu;
+        mu_buffer[i] = pexsi::PEXSI_Solver::pexsi_mu;
     }
 
     this->ParaV = ParaV_in;
     this->ps = std::make_unique<pexsi::PEXSI_Solver>();
 
-    this->DM.resize(nspin);
-    this->EDM.resize(nspin);
-    for (int i = 0; i < nspin; i++)
+    this->DM.resize(this->nspin_dm);
+    this->EDM.resize(this->nspin_dm);
+    for (int i = 0; i < this->nspin_dm; i++)
     {
         this->DM[i] = new T[ParaV->nrow * ParaV->ncol];
         this->EDM[i] = new T[ParaV->nrow * ParaV->ncol];
@@ -48,12 +49,7 @@ DiagoPexsi<T>::DiagoPexsi(const Parallel_Orbitals* ParaV_in)
 template <typename T>
 DiagoPexsi<T>::~DiagoPexsi()
 {
-    int nspin = PARAM.inp.nspin;
-    if (PARAM.inp.nspin == 4)
-    {
-        nspin = 1;
-    }
-    for (int i = 0; i < nspin; i++)
+    for (int i = 0; i < this->nspin_dm; i++)
     {
         delete[] this->DM[i];
         delete[] this->EDM[i];
@@ -67,12 +63,13 @@ void DiagoPexsi<double>::diag(hamilt::Hamilt<double>* phm_in, psi::Psi<double>&
     ModuleBase::TITLE("DiagoPEXSI", "diag");
     matd h_mat, s_mat;
     phm_in->matrix(h_mat, s_mat);
-    std::vector<double> eigen(PARAM.globalv.nlocal, 0.0);
     int ik = psi.get_current_k();
     this->ps->prepare(this->ParaV->blacs_ctxt,
                       this->ParaV->nb,
                       this->ParaV->nrow,
                       this->ParaV->ncol,
+                      this->nlocal,
+                      this->nelec,
                       h_mat.p,
                       s_mat.p,
                       DM[ik],

source/source_hsolver/diago_pexsi.h

@@ -19,7 +19,7 @@ class DiagoPexsi
     static std::vector<double> mu_buffer;
 
   public:
-    DiagoPexsi(const Parallel_Orbitals* ParaV_in);
+    DiagoPexsi(const Parallel_Orbitals* ParaV_in, const int nspin_in, const int nlocal_in, const double nelec_in);
     void diag(hamilt::Hamilt<T>* phm_in, psi::Psi<T>& psi, Real* eigenvalue_in);
     const Parallel_Orbitals* ParaV = nullptr;
     std::vector<T*> DM;
@@ -29,6 +29,14 @@ class DiagoPexsi
     double totalFreeEnergy;
     std::unique_ptr<pexsi::PEXSI_Solver> ps;
     ~DiagoPexsi();
+
+  private:
+    /// number of density matrices to keep: nspin, except that nspin == 4 is
+    /// treated as a single (spinor) density matrix
+    int nspin_dm = 1;
+    /// global dimension of the NAO Hamiltonian
+    int nlocal = 0;
+    double nelec = 0.0;
 };
 } // namespace hsolver
 

source/source_hsolver/hsolver_lcao.cpp

@@ -33,7 +33,6 @@
 #include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_estate/module_dm/density_matrix.h"
 #include "source_hsolver/parallel_k2d.h"
-#include "source_io/module_parameter/parameter.h"
 
 namespace hsolver
 {
@@ -59,13 +58,13 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
             this->parakSolve_cusolver(pHamilt, psi, pes);
         }else 
     #endif
-        if (PARAM.globalv.kpar_lcao > 1
+        if (this->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, nspin);
+            this->parakSolve(pHamilt, psi, pes, this->kpar_lcao, nspin);
         } else
     #endif
-        if (PARAM.globalv.kpar_lcao == 1)
+        if (this->kpar_lcao == 1)
         {
             /// Loop over k points for solve Hamiltonian to eigenpairs(eigenvalues and eigenvectors).
             for (int ik = 0; ik < psi.get_nk(); ++ik)
@@ -113,7 +112,7 @@ void HSolverLCAO<TK, Device>::solve(hamilt::Hamilt<TK>* pHamilt,
     else if (this->method == "pexsi")
     {
 #ifdef __PEXSI // other purification methods should follow this routine
-        DiagoPexsi<TK> pe(ParaV);
+        DiagoPexsi<TK> pe(ParaV, nspin, this->nlocal, this->nelec);
         for (int ik = 0; ik < psi.get_nk(); ++ik)
         {
             /// update H(k) for each k point

source/source_hsolver/hsolver_lcao.h

@@ -15,7 +15,12 @@ template <typename TK, typename Device = base_device::DEVICE_CPU>
 class HSolverLCAO
 {
   public:
-    HSolverLCAO(const Parallel_Orbitals* ParaV_in, std::string method_in) : ParaV(ParaV_in), method(method_in) {};
+    HSolverLCAO(const Parallel_Orbitals* ParaV_in,
+                const std::string method_in,
+                const int kpar_lcao_in,
+                const int nlocal_in,
+                const double nelec_in)
+        : ParaV(ParaV_in), method(method_in), kpar_lcao(kpar_lcao_in), nlocal(nlocal_in), nelec(nelec_in) {};
 
     void solve(hamilt::Hamilt<TK>* pHamilt,
                psi::Psi<TK>& psi,
@@ -40,8 +45,12 @@ class HSolverLCAO
                              elecstate::ElecState* pes);
 
     const Parallel_Orbitals* ParaV = nullptr;
-    
+
     const std::string method;
+
+    const int kpar_lcao; // number of pools for LCAO diagonalization
+    const int nlocal;    // global dimension of the NAO Hamiltonian, only used by the pexsi branch
+    const double nelec;  // total number of electrons, only used by the pexsi branch
 };
 
 } // namespace hsolver

source/source_hsolver/hsolver_lcaopw.cpp

@@ -7,7 +7,6 @@
 #include "source_estate/elecstate_pw.h"
 #include "source_pw/module_pwdft/hamilt_pw.h"
 #include "source_hsolver/diago_iter_assist.h"
-#include "source_io/module_parameter/parameter.h"
 #include "source_estate/elecstate_tools.h"
 #include "source_hamilt/module_xc/exx_info.h"
 
@@ -106,7 +105,7 @@ void HSolverLIP<T>::solve(hamilt::Hamilt<T>* pHamilt, // ESolver_KS_PW::p_hamilt
     elecstate::calEBand(pes->ekb,pes->wg,pes->f_en);
     if (skip_charge)
     {
-        if (PARAM.globalv.use_uspp)
+        if (this->use_uspp)
         {
             reinterpret_cast<elecstate::ElecStatePW<T>*>(pes)->cal_becsum(psi);
         }