refactor the call of initialization of charge density (#6744)

* remove some codes in esolver_ks_lcao and add LCAO_set.cpp

* add LCAO_set.h and LCAO_set.cpp

* finish updating esolver

* update LCAO_set, fix ref of pointer of psi

* update esolver

* delete useless head files

* update esolver

* change effective to eff for short

* add dm2rho in init_dm

* delete useless calculate_weights

* fix a potential bug, the bcast of descf

* move init_rho to KS:before_all_runners

* fix the input parametes of init_rho

* small update of esolver_of.cpp

* fix tests

* update eff in ML_KEDF

* fix a potential bug in allocate() function in charge.cpp

* move set_rhopw from elecstate to before charge.allocate

* fix bug about two allocations of charge class

* fix a bug in OFDFT

* move the common parts in esolver_ks and esolver_of to esolver_fp

* small updates of fp and ks esolvers

* update double xc

* fix bug about allocation of charge density

Author: mohanchen

source/Makefile.Objects

@@ -239,7 +239,6 @@ OBJS_ELECSTAT=elecstate.o\
     potential_types.o\
     pot_sep.o\
     pot_local.o\
-    pot_local_paw.o\
     H_Hartree_pw.o\
     H_TDDFT_pw.o\
     pot_xc.o\

source/source_esolver/esolver_double_xc.cpp

@@ -51,7 +51,6 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
         this->pelec_base = new elecstate::ElecStateLCAO<TK>(&(this->chr_base), // use which parameter?
                                                        &(this->kv),
                                                        this->kv.get_nks(),
-                                                       this->pw_rho,
                                                        this->pw_big);
     }    
 
@@ -87,7 +86,10 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
     this->dmat_base.allocate_dm(&this->kv, &this->pv, PARAM.inp.nspin);
 
     // 10) inititlize the charge density
-    this->chr_base.allocate(PARAM.inp.nspin);
+	this->chr_base.set_rhopw(this->pw_rhod); // mohan add 20251130
+	this->chr_base.allocate(PARAM.inp.nspin);
+	this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
+	this->chr_base.check_rho();
 
     // 11) initialize the potential
     if (this->pelec_base->pot == nullptr)

source/source_esolver/esolver_fp.cpp

@@ -41,24 +41,56 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
     ModuleBase::TITLE("ESolver_FP", "before_all_runners");
 
-    //! read pseudopotentials
+    //! 1) read pseudopotentials
     elecstate::read_pseudo(GlobalV::ofs_running, ucell);
 
-    // setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials
+    //! 2) setup pw_rho, pw_rhod, pw_big, sf, and read_pseudopotentials
     pw::setup_pwrho(ucell, PARAM.globalv.double_grid, this->pw_rho_flag, 
-      this->pw_rho, this->pw_rhod, this->pw_big, 
-      this->classname, inp);
+      this->pw_rho, this->pw_rhod, this->pw_big, this->classname, inp);
 
-    // setup structure factors
+    //! 3) setup structure factors
     this->sf.set(this->pw_rhod, inp.nbspline);
 
-    // write geometry file
+    //! 4)  write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
       ucell, "# Generated by ABACUS ModuleIO::CifParser", "data_?");
 
-    // init charge extrapolation
+    //! 5) init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);
 
+    //! 6) symmetry analysis should be performed every time the cell is changed
+    if (ModuleSymmetry::Symmetry::symm_flag == 1)
+    {
+        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
+        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
+    }
+
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
+
+    //! 7) setup k points in the Brillouin zone according to symmetry.
+    this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
+
+    //! 8) print information
+    ModuleIO::print_parameters(ucell, this->kv, inp);
+
+    //! 9) parallel of FFT grid
+    this->Pgrid.init(this->pw_rhod->nx, this->pw_rhod->ny, this->pw_rhod->nz,
+            this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
+
+    //! 10) calculate the structure factor
+    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
+
+    //! 11) setup the xc functional
+    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
+
+    //! 11) initialize the charge density, we need to first set xc_type,
+    // then we can call chr.allocate()
+	this->chr.set_rhopw(this->pw_rhod); // mohan add 20251130
+	this->chr.allocate(inp.nspin); // mohan move this from setup_estate_pw, 20251128
+
+
     return;
 }
 

source/source_esolver/esolver_gets.cpp

@@ -53,7 +53,6 @@ void ESolver_GetS::before_all_runners(UnitCell& ucell, const Input_para& inp)
         this->pelec = new elecstate::ElecStateLCAO<std::complex<double>>(&(this->chr), // use which parameter?
                                                                          &(this->kv),
                                                                          this->kv.get_nks(),
-                                                                         this->pw_rho,
                                                                          this->pw_big);
     }
 

source/source_esolver/esolver_ks.cpp

@@ -62,46 +62,26 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     // cell_factor
     this->ppcell.cell_factor = inp.cell_factor;
 
-    //! 3) setup charge mixing
+    //! 3) setup Exc for the first element '0' (all elements have same exc) 
+//    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+//    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
+
+    //! 4) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
-
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
-
-    //! 4) setup Exc for the first element '0' (all elements have same exc) 
-    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
-    GlobalV::ofs_running<<XC_Functional::output_info()<<std::endl;
-    
-    //! 5) setup the charge mixing parameters
     p_chgmix->set_mixing(inp.mixing_mode, inp.mixing_beta, inp.mixing_ndim,
       inp.mixing_gg0, inp.mixing_tau, inp.mixing_beta_mag, inp.mixing_gg0_mag,
       inp.mixing_gg0_min, inp.mixing_angle, inp.mixing_dmr, ucell.omega, ucell.tpiba);
-
     p_chgmix->init_mixing();
 
-    //! 6) symmetry analysis should be performed every time the cell is changed
-    if (ModuleSymmetry::Symmetry::symm_flag == 1)
-    {
-        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
-        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
-    }
-
-    //! 7) setup k points in the Brillouin zone according to symmetry.
-    this->kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
-
-    //! 8) print information
-    ModuleIO::print_parameters(ucell, this->kv, inp);
-
-    //! 9) setup plane wave for electronic wave functions
+    //! 5) setup plane wave for electronic wave functions
     pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);
 
-    //! 10) parallel of FFT grid 
-	Pgrid.init(this->pw_rhod->nx, this->pw_rhod->ny, this->pw_rhod->nz,
-			this->pw_rhod->nplane, this->pw_rhod->nrxx, pw_big->nbz, pw_big->bz);
-
-    //! 11) calculate the structure factor
-    this->sf.setup(&ucell, Pgrid, this->pw_rhod);
+    //! 6) read in charge density, mohan add 2025-11-28
+    //! Inititlize the charge density.
+    this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv, this->pw_wfc);
+    this->chr.check_rho(); // check the rho
+  
 }
 
 template <typename T, typename Device>
@@ -115,20 +95,14 @@ void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const i
     this->hamilt2rho_single(ucell, istep, iter, diag_ethr);
 
     // 2) for MPI: STOGROUP? need to rewrite
-    //<Temporary> It may be changed when more clever parallel algorithm is
-    // put forward.
-    // When parallel algorithm for bands are adopted. Density will only be
-    // treated in the first group.
-    //(Different ranks should have abtained the same, but small differences
-    // always exist in practice.)
+    //<Temporary> It may be changed when more clever parallel algorithm is put forward.
+    // When parallel algorithm for bands are adopted. Density will only be treated in the first group.
+    //(Different ranks should have abtained the same, but small differences always exist in practice.)
     // Maybe in the future, density and wavefunctions should use different
     // parallel algorithms, in which they do not occupy all processors, for
     // example wavefunctions uses 20 processors while density uses 10.
     if (PARAM.globalv.ks_run)
     {
-        // double drho = this->estate.caldr2();
-        // EState should be used after it is constructed.
-
         drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
         hsolver_error = 0.0;
         if (iter == 1 && PARAM.inp.calculation != "nscf")
@@ -140,23 +114,16 @@ void ESolver_KS<T, Device>::hamilt2rho(UnitCell& ucell, const int istep, const i
             // so a more precise HSolver should be executed.
             if (hsolver_error > drho)
             {
-                diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running,
-                                                     PARAM.inp.basis_type,
-                                                     PARAM.inp.esolver_type,
-                                                     PARAM.inp.precision,
-                                                     hsolver_error,
-                                                     drho,
-                                                     diag_ethr,
-                                                     PARAM.inp.nelec);
+                diag_ethr = hsolver::reset_diag_ethr(GlobalV::ofs_running, PARAM.inp.basis_type,
+                            PARAM.inp.esolver_type, PARAM.inp.precision, hsolver_error,
+                            drho, diag_ethr, PARAM.inp.nelec);
 
                 this->hamilt2rho_single(ucell, istep, iter, diag_ethr);
 
                 drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
 
                 hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type,
-                                                          PARAM.inp.esolver_type,
-                                                          diag_ethr,
-                                                          PARAM.inp.nelec);
+                                PARAM.inp.esolver_type, diag_ethr, PARAM.inp.nelec);
             }
         }
     }
@@ -168,22 +135,17 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     ModuleBase::TITLE("ESolver_KS", "runner");
     ModuleBase::timer::tick(this->classname, "runner");
 
-    //----------------------------------------------------------------
     // 1) before_scf (electronic iteration loops)
-    //----------------------------------------------------------------
     this->before_scf(ucell, istep);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
-    //----------------------------------------------------------------
     // 2) SCF iterations
-    //----------------------------------------------------------------
     bool conv_esolver = false;
     this->niter = this->maxniter;
     this->diag_ethr = PARAM.inp.pw_diag_thr;
     this->scf_nmax_flag = false; // mohan add 2025-09-21
     for (int iter = 1; iter <= this->maxniter; ++iter)
 	{
-        // mohan add 2025-09-21
 		if(iter == this->maxniter)
 		{
 			this->scf_nmax_flag=true;

source/source_esolver/esolver_ks_lcao.cpp

@@ -21,13 +21,6 @@
 #include "source_lcao/rho_tau_lcao.h" // mohan add 20251024
 #include "source_lcao/LCAO_set.h" // mohan add 20251111
 
-
-// tmp
-#include "source_psi/setup_psi.h" // use Setup_Psi
-#include "source_io/read_wfc_nao.h" // use read_wfc_nao
-#include "source_estate/elecstate_tools.h" // use fixed_weights
-
-
 namespace ModuleESolver
 {
 
@@ -61,7 +54,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     {
         // TK stands for double and std::complex<double>?
         this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), &(this->kv),
-          this->kv.get_nks(), this->pw_rho, this->pw_big);
+          this->kv.get_nks(), this->pw_big);
     }
 
     // 3) read LCAO orbitals/projectors and construct the interpolation tables.
@@ -458,8 +451,13 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
             }
             this->dftu.cal_energy_correction(ucell, istep);
         }
-        this->dftu.output(ucell);
-    }
+		this->dftu.output(ucell);
+		// use the converged occupation matrix for next MD/Relax SCF calculation
+		if (conv_esolver)
+		{
+			this->dftu.initialed_locale = true;
+		}
+	}
 
     // 2) for deepks, calculate delta_e, output labels during electronic steps
     this->deepks.delta_e(ucell, this->kv, this->orb_, this->pv, this->gd, dm_vec, this->pelec->f_en, PARAM.inp);
@@ -487,12 +485,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
         }
     }
 
-    // use the converged occupation matrix for next MD/Relax SCF calculation
-    if (PARAM.inp.dft_plus_u && conv_esolver)
-    {
-        this->dftu.initialed_locale = true;
-    }
-
     // control the output related to the finished iteration
     ModuleIO::ctrl_iter_lcao<TK, TR>(ucell, PARAM.inp, this->kv, this->pelec, *this->dmat.dm,
       this->pv, this->gd, this->psi, this->chr, this->p_chgmix, 

source/source_esolver/esolver_ks_pw.cpp

@@ -80,7 +80,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
     //! Call before_all_runners() of ESolver_KS
     ESolver_KS<T, Device>::before_all_runners(ucell, inp);
 
-    //! setup and allocation for pelec, charge density, potentials, etc. 
+    //! setup and allocation for pelec, potentials, etc. 
     elecstate::setup_estate_pw<T, Device>(ucell, this->kv, this->sf, this->pelec, this->chr,
       this->locpp, this->ppcell, this->vsep_cell, this->pw_wfc, this->pw_rho,
       this->pw_rhod, this->pw_big, this->solvent, inp);
@@ -286,7 +286,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     // pp projectors, liuyu 2023-10-24
     if (PARAM.globalv.use_uspp)
     {
-        ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
+        ModuleBase::matrix veff = this->pelec->pot->get_eff_v();
         this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
     }
 

source/source_esolver/esolver_of.cpp

@@ -73,40 +73,15 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
+//  XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
     int func_type = XC_Functional::get_func_type();
     if (func_type > 2)
     {
         ModuleBase::WARNING_QUIT("esolver_of", "meta-GGA and Hybrid functionals are not supported by OFDFT.");
     }
 
-    // symmetry analysis should be performed every time the cell is changed
-    if (ModuleSymmetry::Symmetry::symm_flag == 1)
-    {
-        ucell.symm.analy_sys(ucell.lat, ucell.st, ucell.atoms, GlobalV::ofs_running);
-        ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SYMMETRY");
-    }
-
-    // Setup the k points according to symmetry.
-    kv.set(ucell,ucell.symm, inp.kpoint_file, inp.nspin, ucell.G, ucell.latvec, GlobalV::ofs_running);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT K-POINTS");
-
-    // print information
-    // mohan add 2021-01-30
-    ModuleIO::print_parameters(ucell, kv, inp);
-
-    // initialize the real-space uniform grid for FFT and parallel
-    // distribution of plane waves
-    Pgrid.init(pw_rho->nx,
-                        pw_rho->ny,
-                        pw_rho->nz,
-                        pw_rho->nplane,
-                        pw_rho->nrxx,
-                        pw_big->nbz,
-                        pw_big->bz); // mohan add 2010-07-22, update 2011-05-04
-    // Calculate Structure factor
-    sf.setup(&ucell, Pgrid, pw_rho);
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT BASIS");
+    this->chr.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
+    this->chr.check_rho(); // check the rho
 
     // initialize local pseudopotential
     this->locpp.init_vloc(ucell,pw_rho);
@@ -122,7 +97,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     // liuyu move here 2023-10-09
     // D in uspp need vloc, thus behind init_scf()
     // calculate the effective coefficient matrix for non-local pseudopotential projectors
-    ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
+    ModuleBase::matrix veff = this->pelec->pot->get_eff_v();
 
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT POTENTIAL");
 
@@ -214,6 +189,8 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
     //! 1) call before_scf() of ESolver_FP
     ESolver_FP::before_scf(ucell, istep);
 
+
+
     if (ucell.cell_parameter_updated)
     {
         this->dV_ = ucell.omega / this->pw_rho->nxyz;
@@ -317,10 +294,10 @@ void ESolver_OF::update_potential(UnitCell& ucell)
     this->kedf_manager_->get_potential(this->chr.rho,
                                        this->pphi_,
                                        this->pw_rho,
-                                       this->pelec->pot->get_effective_v()); // KEDF potential
+                                       this->pelec->pot->get_eff_v()); // KEDF potential
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
-        const double* vr_eff = this->pelec->pot->get_effective_v(is);
+        const double* vr_eff = this->pelec->pot->get_eff_v(is);
         for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
         {
             this->pdEdphi_[is][ir] = vr_eff[ir];
@@ -516,9 +493,9 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell, const bool conv_eso
     this->kedf_manager_->get_potential(this->chr.rho,
                                        this->pphi_,
                                        this->pw_rho,
-                                       this->pelec->pot->get_effective_v()); // KEDF potential
+                                       this->pelec->pot->get_eff_v()); // KEDF potential
 
-        const double* vr_eff = this->pelec->pot->get_effective_v(0);
+        const double* vr_eff = this->pelec->pot->get_eff_v(0);
         for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
         {
             this->pdEdphi_[0][ir] = vr_eff[ir];