Update some names of deepks tests and refactor some ofdft codes (#6754)

* 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

* delete istep in init_pot and init_scf

* add an input parameter in chr.allocate

* update directory names of deepks

* update some small places

* add of_print_info in ofdft directory and delete the function in esolver_of

* move ofdft print_info to module_ofdft

* fix bugs in tests

* add header

* small update of timer in ofdft

* fix a bug related to timer

* update

* update

Author: mohanchen

source/Makefile.Objects

@@ -738,6 +738,7 @@ OBJS_SRCPW=H_Ewald_pw.o\
     stress_func_onsite.o\
     stress_pw.o\
     of_stress_pw.o\
+    of_print_info.o\
     symmetry_rho.o\
     symmetry_rhog.o\
     setup_psi_pw.o\

source/source_esolver/esolver_double_xc.cpp

@@ -87,7 +87,8 @@ void ESolver_DoubleXC<TK, TR>::before_all_runners(UnitCell& ucell, const Input_p
 
     // 10) inititlize the charge density
 	this->chr_base.set_rhopw(this->pw_rhod); // mohan add 20251130
-	this->chr_base.allocate(PARAM.inp.nspin);
+    const bool kin_den = this->chr_base.kin_density(); // mohan add 20251202
+	this->chr_base.allocate(PARAM.inp.nspin, kin_den);
 	this->chr_base.init_rho(ucell, this->Pgrid, this->sf.strucFac, ucell.symm, &this->kv);
 	this->chr_base.check_rho();
 
@@ -156,7 +157,7 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
 	}
 
     XC_Functional::set_xc_type(PARAM.inp.deepks_out_base);
-    this->pelec_base->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    this->pelec_base->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 
     // DMR should be same size with Hamiltonian(R)

source/source_esolver/esolver_fp.cpp

@@ -88,7 +88,8 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
     //! 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
+    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
 
 
     return;

source/source_esolver/esolver_ks.cpp

@@ -62,22 +62,18 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     // cell_factor
     this->ppcell.cell_factor = inp.cell_factor;
 
-    //! 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
+    //! 3) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
     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();
 
-    //! 5) setup plane wave for electronic wave functions
+    //! 4) setup plane wave for electronic wave functions
     pw::setup_pwwfc(inp, ucell, *this->pw_rho, this->kv, this->pw_wfc);
 
-    //! 6) read in charge density, mohan add 2025-11-28
+    //! 5) 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

source/source_esolver/esolver_ks_lcao.cpp

@@ -173,7 +173,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
         LCAO_domain::dm2rho(this->dmat.dm->get_DMR_vector(), PARAM.inp.nspin, this->pelec->charge, true);
     }
     // 12.2) init_scf, should be before_scf? mohan add 2025-03-10
-    this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
     // 13) initalize DM(R), which has the same size with Hamiltonian(R)
     auto* hamilt_lcao = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);

source/source_esolver/esolver_ks_lcao_tddft.cpp

@@ -133,7 +133,7 @@ void ESolver_KS_LCAO_TDDFT<TR, Device>::runner(UnitCell& ucell, const int istep)
                                         GlobalV::ofs_running,
                                         GlobalV::ofs_warning);
             // need to test if correct when estep>0
-            this->pelec->init_scf(totstep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+            this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
             if (totstep <= PARAM.inp.td_tend + 1)
             {

source/source_esolver/esolver_of.cpp

@@ -13,6 +13,8 @@
 #include "source_estate/cal_ux.h"
 #include "source_pw/module_pwdft/forces.h"
 #include "source_pw/module_ofdft/of_stress_pw.h"
+// mohan add
+#include "source_pw/module_ofdft/of_print_info.h"
 
 namespace ModuleESolver
 {
@@ -92,7 +94,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     this->init_elecstate(ucell);
 
     // calculate the total local pseudopotential in real space
-    this->pelec->init_scf(0, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs
+    this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs
 
     // liuyu move here 2023-10-09
     // D in uspp need vloc, thus behind init_scf()
@@ -206,20 +208,19 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
 
         // Refresh the arrays
         delete this->psi_;
-        this->psi_ = new psi::Psi<double>(1, 
-                                          PARAM.inp.nspin, 
-                                          this->pw_rho->nrxx,
-                                          this->pw_rho->nrxx,
-                                          true);
+        this->psi_ = new psi::Psi<double>(1, PARAM.inp.nspin, 
+                                          this->pw_rho->nrxx, this->pw_rho->nrxx, true);
+
         for (int is = 0; is < PARAM.inp.nspin; ++is)
         {
             this->pphi_[is] = this->psi_->get_pointer(is);
         }
 
         delete this->ptemp_rho_;
         this->ptemp_rho_ = new Charge();
-        this->ptemp_rho_->set_rhopw(this->pw_rho);
-        this->ptemp_rho_->allocate(PARAM.inp.nspin);
+		this->ptemp_rho_->set_rhopw(this->pw_rho);
+		const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
+		this->ptemp_rho_->allocate(PARAM.inp.nspin, kin_den);
 
         for (int is = 0; is < PARAM.inp.nspin; ++is)
         {
@@ -234,7 +235,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
         }
     }
 
-    this->pelec->init_scf(istep, ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
+    this->pelec->init_scf(ucell, Pgrid, sf.strucFac, locpp.numeric, ucell.symm);
 
     Symmetry_rho srho;
     for (int is = 0; is < PARAM.inp.nspin; is++)
@@ -431,7 +432,8 @@ bool ESolver_OF::check_exit(bool& conv_esolver)
     conv_esolver = (this->of_conv_ == "energy" && energyConv) || (this->of_conv_ == "potential" && potConv)
                          || (this->of_conv_ == "both" && potConv && energyConv);
 
-    this->print_info(conv_esolver);
+    OFDFT::print_info(this->iter_, this->iter_time, this->energy_current_, this->energy_last_, 
+                      this->normdLdphi_, this->pelec, this->kedf_manager_, conv_esolver);
 
     if (conv_esolver || this->iter_ >= this->max_iter_)
     {
@@ -567,11 +569,8 @@ void ESolver_OF::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress)
 {
     ModuleBase::matrix kinetic_stress_;
     kinetic_stress_.create(3, 3);
-    this->kedf_manager_->get_stress(ucell.omega,
-                                    this->chr.rho,
-                                    this->pphi_,
-                                    this->pw_rho,
-                                    kinetic_stress_); // kinetic stress
+    this->kedf_manager_->get_stress(ucell.omega, this->chr.rho,
+                         this->pphi_, this->pw_rho, kinetic_stress_); // kinetic stress
 
     OF_Stress_PW ss(this->pelec, this->pw_rho);
     ss.cal_stress(stress, kinetic_stress_, ucell, &ucell.symm, this->locpp, &sf, &kv);

source/source_esolver/esolver_of.h

@@ -98,9 +98,6 @@ class ESolver_OF : public ESolver_FP
     void check_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
     void test_direction(double* dEdtheta, double** ptemp_phi, UnitCell& ucell);
 
-    // --------------------- output the necessary information -----------
-    void print_info(const bool conv_esolver);
-
     // --------------------- interface to blas --------------------------
     double inner_product(double* pa, double* pb, int length, double dV = 1)
     {

source/source_esolver/esolver_of_tool.cpp

@@ -86,7 +86,8 @@ void ESolver_OF::allocate_array()
     delete this->ptemp_rho_;
     this->ptemp_rho_ = new Charge();
     this->ptemp_rho_->set_rhopw(this->pw_rho);
-    this->ptemp_rho_->allocate(PARAM.inp.nspin);
+    const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
+    this->ptemp_rho_->allocate(PARAM.inp.nspin, kin_den);
 
     this->theta_ = new double[PARAM.inp.nspin];
     this->pdLdphi_ = new double*[PARAM.inp.nspin];
@@ -378,140 +379,13 @@ void ESolver_OF::test_direction(double* dEdtheta, double** ptemp_phi, UnitCell&
             Parallel_Reduce::reduce_all(pseudopot_energy);
             temp_energy += kinetic_energy + pseudopot_energy;
             GlobalV::ofs_warning << i << "    " << dEdtheta[0] << "    " << temp_energy << std::endl;
-            if (this->theta_[0] == 0) {
-                std::cout << "dEdtheta    " << dEdtheta[0] << std::endl;
-}
-        }
+			if (this->theta_[0] == 0) 
+			{
+				std::cout << "dEdtheta    " << dEdtheta[0] << std::endl;
+			}
+		}
         exit(0);
     }
 }
 
-/**
- * @brief Print nessecary information to the screen,
- * and write the components of the total energy into running_log.
- */
-void ESolver_OF::print_info(const bool conv_esolver)
-{
-    if (this->iter_ == 0)
-    {
-        std::cout << " ============================= Running OFDFT "
-                     "=============================="
-                  << std::endl;
-        std::cout << " ITER       ETOT/eV           EDIFF/eV        EFERMI/eV    POTNORM   TIME/s"
-                  << std::endl;
-    }
-
-    std::map<std::string, std::string> prefix_map = {
-        {"cg1", "CG"},
-        {"cg2", "CG"},
-        {"tn", "TN"}
-    };
-    std::string iteration = prefix_map[PARAM.inp.of_method] + std::to_string(this->iter_);
-#ifdef __MPI
-    double duration = (double)(MPI_Wtime() - this->iter_time);
-#else
-    double duration
-        = (std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now() - this->iter_time)).count()
-          / static_cast<double>(1e6);
-#endif
-    std::cout << " " << std::setw(8) << iteration
-              << std::setw(18) << std::scientific << std::setprecision(8) << this->energy_current_ * ModuleBase::Ry_to_eV
-              << std::setw(18) << (this->energy_current_ - this->energy_last_) * ModuleBase::Ry_to_eV
-              << std::setw(13) << std::setprecision(4) << this->pelec->eferm.get_efval(0) * ModuleBase::Ry_to_eV
-              << std::setw(13) << std::setprecision(4) << this->normdLdphi_
-              << std::setw(6) << std::fixed << std::setprecision(2) << duration << std::endl;
-
-    GlobalV::ofs_running << std::setprecision(12);
-    GlobalV::ofs_running << std::setiosflags(std::ios::right);
-
-    GlobalV::ofs_running << "\nIter" << this->iter_ << ": the norm of potential is " << this->normdLdphi_ << std::endl;
-
-    std::vector<std::string> titles;
-    std::vector<double> energies_Ry;
-    std::vector<double> energies_eV;
-	if ((PARAM.inp.out_band[0] > 0 && 
-				((this->iter_ + 1) % PARAM.inp.out_band[0] == 0 || 
-				 conv_esolver || 
-				 this->iter_ == PARAM.inp.scf_nmax)) || 
-			PARAM.inp.init_chg == "file")
-    {
-        titles.push_back("E_Total");
-        energies_Ry.push_back(this->pelec->f_en.etot);
-        titles.push_back("E_Kinetic");
-        energies_Ry.push_back(this->pelec->f_en.ekinetic);
-        titles.push_back("E_Hartree");
-        energies_Ry.push_back(this->pelec->f_en.hartree_energy);
-        titles.push_back("E_xc");
-        energies_Ry.push_back(this->pelec->f_en.etxc - this->pelec->f_en.etxcc);
-        titles.push_back("E_LocalPP");
-        energies_Ry.push_back(this->pelec->f_en.e_local_pp);
-        titles.push_back("E_Ewald");
-        energies_Ry.push_back(this->pelec->f_en.ewald_energy);
-
-        this->kedf_manager_->record_energy(titles, energies_Ry);
-        
-        std::string vdw_method = PARAM.inp.vdw_method;
-        if (vdw_method == "d2") // Peize Lin add 2014-04, update 2021-03-09
-        {
-            titles.push_back("E_vdwD2");
-            energies_Ry.push_back(this->pelec->f_en.evdw);
-        }
-        else if (vdw_method == "d3_0" || vdw_method == "d3_bj") // jiyy add 2019-05, update 2021-05-02
-        {
-            titles.push_back("E_vdwD3");
-            energies_Ry.push_back(this->pelec->f_en.evdw);
-        }
-        if (PARAM.inp.imp_sol)
-        {
-            titles.push_back("E_sol_el");
-            energies_Ry.push_back(this->pelec->f_en.esol_el);
-            titles.push_back("E_sol_cav");
-            energies_Ry.push_back(this->pelec->f_en.esol_cav);
-        }
-        if (PARAM.inp.efield_flag)
-        {
-            titles.push_back("E_efield");
-            energies_Ry.push_back(elecstate::Efield::etotefield);
-        }
-        if (PARAM.inp.gate_flag)
-        {
-            titles.push_back("E_gatefield");
-            energies_Ry.push_back(elecstate::Gatefield::etotgatefield);
-        }
-    }
-    else
-    {
-        titles.push_back("E_Total");
-        energies_Ry.push_back(this->pelec->f_en.etot);
-    }
-
-    if (PARAM.globalv.two_fermi)
-    {
-        titles.push_back("E_Fermi_up");
-        energies_Ry.push_back(this->pelec->eferm.get_efval(0));
-        titles.push_back("E_Fermi_dw");
-        energies_Ry.push_back(this->pelec->eferm.get_efval(1));
-    }
-    else
-    {
-        titles.push_back("E_Fermi");
-        energies_Ry.push_back(this->pelec->eferm.get_efval(0));
-    }
-    energies_eV.resize(energies_Ry.size());
-    std::transform(energies_Ry.begin(), energies_Ry.end(), energies_eV.begin(), [](double energy) {
-        return energy * ModuleBase::Ry_to_eV;
-    });
-    FmtTable table(/*titles=*/{"Energy", "Rydberg", "eV"}, 
-                   /*nrows=*/titles.size(), 
-                   /*formats=*/{"%20s", "%20.12f", "%20.12f"}, 0);
-    table << titles << energies_Ry << energies_eV;
-    GlobalV::ofs_running << table.str() << std::endl;
-
-    // reset the iter_time for the next iteration
-#ifdef __MPI
-    this->iter_time = MPI_Wtime();
-#else
-    this->iter_time = std::chrono::system_clock::now();
-#endif
-}
 } // namespace ModuleESolver

source/source_esolver/lcao_others.cpp

@@ -219,7 +219,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     // pelec should be initialized before these calculations
-    this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
+    this->pelec->init_scf(ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
 
     // self consistent calculations for electronic ground state
     if (cal_type == "get_pchg")