source_estate: pass Charge its nspin, nelec and verbosity explicitly, and add the two seams charge_test needed elsewhere

charge_test needed the macro for three separate reasons, in three different
modules. All three are fixed here; the module-crossing ones are called out below
because they touch modules that currently have no access hacks at all.

1) source_estate -- Charge reads three INPUT values (19 test writes, 3 keys)

  charge.cpp read nspin, nelec and test_charge from the global singleton in five
  methods. The obvious move -- give renormalize_rho() an nelec parameter --
  cascades badly: of its six production call sites, rho_tau_lcao.cpp has no PARAM
  reference at all, so the value would have to be threaded through
  LCAO_domain::dm2rho (six more sites) and then ElecStateLCAO::dm2rho into
  hsolver_lcao.cpp -- about fifteen sites across three modules, and it would put a
  global read into a file that is clean today.

  Charge already keeps nspin as state set by allocate(), so nelec is kept the same
  way. allocate() takes test_charge and nelec_in; renormalize_rho() and
  check_rho() read this->nelec, so their signatures are untouched and nothing
  cascades. save_rho_before_sum_band() takes nspin, init_final_scf() takes nspin
  and test_charge.

  allocate() has nine production call sites, not the two a grep on guessed object
  names suggested; the others reach it through chg1_, chr_base, chg_gs,
  ptemp_rho_, chg_01 and chg_drho[d]. Seven of the nine now pass an injected inp_
  or the source Charge's own nspin/nelec, so only veff_dh.cpp adds PARAM
  references, and it already read PARAM.inp.nspin on the same line.

  Two const getters, get_allocate_rho() and get_allocate_rho_final_scf(), cover
  the 14 read-only accesses to the two allocation flags.

2) source_hamilt/module_xc -- XC_Functional::func_type and ked_flag (6 writes)

  Both are private statics with public getters and no setters. Added
  set_func_type() and set_ked_flag(), matching the get_func_type() /
  get_ked_flag() pair and the existing set_hybrid_alpha() / set_hse_omega() style
  in the same header. The test's out-of-line definitions of the two statics need
  no seam -- a definition does not require access.

3) source_basis/module_pw -- PW_Basis::distribute_r() and distribute_g()

  No new API. The test hand-rolled initgrids / distribute_r / initparameters /
  distribute_g; pw_basis.h documents the public sequence as initgrids /
  initparameters / setuptransform(), and setuptransform() performs both
  distribute calls itself. The test now follows the documented order.

Three EXPECT_EQ(PARAM.input.nelec, 8) assertions guarded the 8.0 that the
renormalize/check assertions below them expect; they now assert the same thing
about the value the test configures.

The five test_charge/nspin/nelec reads left in charge.cpp are inside atomic_rho(),
which charge_test does not call, so injecting them would cost budget without
removing a macro.

No production logic changed. No expected value changed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: Critsium-xy

source/source_esolver/esolver_dfpt_pw.cpp

@@ -45,7 +45,7 @@ class XC_First_Order_FDM : public ModuleDFPT::XC_First_Order
         fdm_ = new elecstate::PotXC_FDM(rho_basis, chg0, ucell);
         chg1_ = new Charge();
         chg1_->set_rhopw(rho_basis);
-        chg1_->allocate(chg0->nspin, false);
+        chg1_->allocate(chg0->nspin, false, /*test_charge=*/0, chg0->nelec);
         veff_1_.create(chg0->nspin, chg0->nrxx);
     }
 

source/source_esolver/esolver_double_xc.cpp

@@ -92,7 +92,7 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(BaseCell& basecell, const Inpu
     // 10) inititlize the charge density
     this->chr_base.set_rhopw(this->pw_rhod);           // mohan add 20251130
     const bool kin_den = this->chr_base.kin_density(); // mohan add 20251202
-    this->chr_base.allocate(this->inp_->nspin, kin_den);
+    this->chr_base.allocate(this->inp_->nspin, kin_den, this->inp_->test_charge, this->inp_->nelec);
     this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
     this->chr_base.check_rho();
 

source/source_esolver/esolver_fp.cpp

@@ -148,7 +148,7 @@ void ESolver_FP::before_all_runners(BaseCell& basecell, const Input_para& inp)
     // then we can call chr.allocate()
 	this->chr.set_rhopw(this->pw_rhod); // mohan add 20251130
     const bool kin_den = this->chr.kin_density(); // mohan add 20251202
-	this->chr.allocate(inp.nspin, kin_den); // mohan move this from setup_estate_pw, 20251128
+	this->chr.allocate(inp.nspin, kin_den, inp.test_charge, inp.nelec); // mohan move this from setup_estate_pw, 20251128
 
 
     return;

source/source_esolver/esolver_ks.cpp

@@ -210,7 +210,7 @@ void ESolver_KS::iter_init(UnitCell& ucell, const int istep, const int iter)
     }
 
     // save input charge density (rho)
-    this->chr.save_rho_before_sum_band();
+    this->chr.save_rho_before_sum_band(this->inp_->nspin);
 }
 
 void ESolver_KS::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver)

source/source_esolver/esolver_lr_lcao_tddft.cpp

@@ -823,7 +823,7 @@ void ModuleESolver::ESolver_LR<T, TR>::read_ks_chg(Charge& chg_gs)
 {
     chg_gs.set_rhopw(this->pw_rho);
     const bool kin_den = chg_gs.kin_density(); // mohan add 20251202
-    chg_gs.allocate(this->nspin, kin_den);
+    chg_gs.allocate(this->nspin, kin_den, this->inp_->test_charge, this->inp_->nelec);
     GlobalV::ofs_running << " try to read charge from file : ";
     for (int is = 0; is < this->nspin; ++is)
     {

source/source_esolver/esolver_of.cpp

@@ -220,7 +220,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
         this->ptemp_rho_ = new Charge();
 		this->ptemp_rho_->set_rhopw(this->pw_rho);
 		const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
-		this->ptemp_rho_->allocate(this->inp_->nspin, kin_den);
+		this->ptemp_rho_->allocate(this->inp_->nspin, kin_den, this->inp_->test_charge, this->inp_->nelec);
 
         for (int is = 0; is < this->inp_->nspin; ++is)
         {

source/source_esolver/esolver_of_tool.cpp

@@ -90,7 +90,7 @@ void ESolver_OF::allocate_array()
     this->ptemp_rho_ = new Charge();
     this->ptemp_rho_->set_rhopw(this->pw_rho);
     const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
-    this->ptemp_rho_->allocate(this->inp_->nspin, kin_den);
+    this->ptemp_rho_->allocate(this->inp_->nspin, kin_den, this->inp_->test_charge, this->inp_->nelec);
 
     this->theta_ = new double[this->inp_->nspin];
     this->pdLdphi_ = new double*[this->inp_->nspin];

source/source_estate/module_charge/charge.cpp

@@ -86,7 +86,7 @@ void Charge::destroy()
     }
 }
 
-void Charge::allocate(const int& nspin_in, const bool kin_den)
+void Charge::allocate(const int& nspin_in, const bool kin_den, const int test_charge, const double& nelec_in)
 {
     ModuleBase::TITLE("Charge", "allocate");
 
@@ -110,8 +110,9 @@ void Charge::allocate(const int& nspin_in, const bool kin_den)
 
     //  mohan add 2021-02-20
     this->nspin = nspin_in;
+    this->nelec = nelec_in;
 
-    if (PARAM.inp.test_charge > 1)
+    if (test_charge > 1)
     {
         std::cout << "\n spin_number = " << nspin << " real_point_number = " << nrxx << std::endl;
     }
@@ -217,7 +218,7 @@ void Charge::renormalize_rho()
     const double sr = this->sum_rho();
     GlobalV::ofs_warning << std::setprecision(15);
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_warning, "charge before normalized", sr);
-    const double normalize_factor = PARAM.inp.nelec / sr;
+    const double normalize_factor = this->nelec / sr;
 
     for (int is = 0; is < nspin; is++)
     {
@@ -651,9 +652,9 @@ void Charge::atomic_rho(const int spin_number_need,
     return;
 }
 
-void Charge::save_rho_before_sum_band()
+void Charge::save_rho_before_sum_band(const int nspin)
 {
-    for (int is = 0; is < PARAM.inp.nspin; is++)
+    for (int is = 0; is < nspin; is++)
     {
         ModuleBase::GlobalFunc::DCOPY(rho[is], rho_save[is], this->rhopw->nrxx);
         if (XC_Functional::get_ked_flag())
@@ -687,7 +688,7 @@ void Charge::check_rho()
     {
         double ne = 0.0;
         ne = this->cal_rho2ne(rho[0]);
-        if (std::abs(ne - PARAM.inp.nelec) > 1.0e-6)
+        if (std::abs(ne - this->nelec) > 1.0e-6)
         {
             ModuleBase::WARNING("Charge", "Charge is not equal to the number of electrons!");
         }
@@ -709,31 +710,31 @@ void Charge::check_rho()
             ModuleBase::WARNING_QUIT("Charge", "Number of spin-up electrons set in starting magnetization exceeds all available.");
         }
         // for total charge
-        if (std::abs(ne_up + ne_dn - PARAM.inp.nelec) > 1.0e-6)
+        if (std::abs(ne_up + ne_dn - this->nelec) > 1.0e-6)
         {
             ModuleBase::WARNING("Charge", "Charge is not equal to the number of electrons!");
         }
     }
 }
 
 // LiuXh add 20180619
-void Charge::init_final_scf()
+void Charge::init_final_scf(const int nspin, const int test_charge)
 {
     ModuleBase::TITLE("Charge", "init_after_scf");
 
     assert(allocate_rho_final_scf == false);
-    if (PARAM.inp.test_charge > 1)
+    if (test_charge > 1)
     {
-        std::cout << "\n spin_number = " << PARAM.inp.nspin << " real_point_number = " << this->rhopw->nrxx << std::endl;
+        std::cout << "\n spin_number = " << nspin << " real_point_number = " << this->rhopw->nrxx << std::endl;
     }
 
     // allocate memory
-    rho = new double*[PARAM.inp.nspin];
-    rhog = new std::complex<double>*[PARAM.inp.nspin];
-    rho_save = new double*[PARAM.inp.nspin];
-    rhog_save = new std::complex<double>*[PARAM.inp.nspin];
+    rho = new double*[nspin];
+    rhog = new std::complex<double>*[nspin];
+    rho_save = new double*[nspin];
+    rhog_save = new std::complex<double>*[nspin];
 
-    for (int is = 0; is < PARAM.inp.nspin; is++)
+    for (int is = 0; is < nspin; is++)
     {
         rho[is] = new double[this->rhopw->nrxx];
         rhog[is] = new std::complex<double>[this->rhopw->npw];
@@ -745,10 +746,10 @@ void Charge::init_final_scf()
         ModuleBase::GlobalFunc::ZEROS(rhog_save[is], this->rhopw->npw);
     }
 
-    ModuleBase::Memory::record("Chg::rho", sizeof(double) * PARAM.inp.nspin * this->rhopw->nrxx);
-    ModuleBase::Memory::record("Chg::rho_save", sizeof(double) * PARAM.inp.nspin * this->rhopw->nrxx);
-    ModuleBase::Memory::record("Chg::rhog", sizeof(double) * PARAM.inp.nspin * this->rhopw->npw);
-    ModuleBase::Memory::record("Chg::rhog_save", sizeof(double) * PARAM.inp.nspin * this->rhopw->npw);
+    ModuleBase::Memory::record("Chg::rho", sizeof(double) * nspin * this->rhopw->nrxx);
+    ModuleBase::Memory::record("Chg::rho_save", sizeof(double) * nspin * this->rhopw->nrxx);
+    ModuleBase::Memory::record("Chg::rhog", sizeof(double) * nspin * this->rhopw->npw);
+    ModuleBase::Memory::record("Chg::rhog_save", sizeof(double) * nspin * this->rhopw->npw);
 
     this->rho_core = new double[this->rhopw->nrxx]; // core charge in real space
     ModuleBase::GlobalFunc::ZEROS(rho_core, this->rhopw->nrxx);

source/source_estate/module_charge/charge.h

@@ -86,7 +86,11 @@ class Charge
     // mohan add 2025-12-02
     bool kin_density() const;
 
-    void allocate(const int &nspin_in, const bool kin_den);
+    /// @param test_charge INPUT verbosity for charge diagnostics
+    /// @param nelec target number of electrons; kept as state, like nspin,
+    ///        because renormalize_rho() and check_rho() need it long after
+    ///        allocate() has returned
+    void allocate(const int &nspin_in, const bool kin_den, const int test_charge, const double& nelec_in);
 
     void atomic_rho(const int spin_number_need,
                     const double& omega,
@@ -102,7 +106,8 @@ class Charge
 
     double sum_rho() const;
 
-    void save_rho_before_sum_band();
+    /// @param nspin number of spin channels to copy
+    void save_rho_before_sum_band(const int nspin);
 
 	// for non-linear core correction
     void non_linear_core_correction
@@ -121,7 +126,25 @@ class Charge
 
     void check_rho(); // to check whether the charge density is normal
 
-    void init_final_scf(); //LiuXh add 20180619
+    /// @param nspin number of spin channels to allocate
+    /// @param test_charge INPUT verbosity for charge diagnostics
+    void init_final_scf(const int nspin, const int test_charge); //LiuXh add 20180619
+
+    //====================================================================
+    // Test seam: the two allocation flags are read-only state that
+    // charge_test checks after each allocate()/init_final_scf() call.
+    //====================================================================
+
+    /// @brief whether allocate() has claimed the rho arrays
+    bool get_allocate_rho() const
+    {
+        return allocate_rho;
+    }
+    /// @brief whether init_final_scf() has claimed the rho arrays
+    bool get_allocate_rho_final_scf() const
+    {
+        return allocate_rho_final_scf;
+    }
 
 	public:
     /**
@@ -152,6 +175,7 @@ class Charge
     int nxyz = 0; // total number of r vectors
     int ngmc=0; // number of g vectors in this processor
     int nspin=0; // number of spins
+    double nelec=0.0; // target number of electrons, set by allocate()
     ModulePW::PW_Basis* rhopw = nullptr;// When double_grid is used, rhopw = rhodpw (dense grid)
     bool cal_elf = false; // whether to calculate electron localization function (ELF)
 

source/source_estate/module_pot/pot_xc_fdm.cpp

@@ -49,7 +49,7 @@ void PotXC_FDM::cal_v_eff(
 
 	Charge chg_01;
 	chg_01.set_rhopw(chg_1->rhopw);
-	chg_01.allocate(chg_1->nspin, chg_01.kin_density());
+	chg_01.allocate(chg_1->nspin, chg_01.kin_density(), /*test_charge=*/0, chg_1->nelec);
 
 	for(int ir=0; ir<chg_01.nrxx; ++ir)
 	{