Update ESolver, delete pelec->charge, instead using chr defined in ESolver_fp directly (deepmodeling#5963)

* reorganize esolver_ks_lcao.h

* delete some useless files in esolver_ks.h

* change pelec->charge to this->chr

* will delete pelec->charge in near future, we just directly use chr defined in esolver_fp

* fix problem in ML_KEDF

* update this->chr

Author: mohanchen

source/module_esolver/esolver_fp.cpp

@@ -23,15 +23,16 @@ namespace ModuleESolver
 
 ESolver_FP::ESolver_FP()
 {
-    // pw_rho = new ModuleBase::PW_Basis();
-    // LCAO basis doesn't support GPU acceleration on FFT currently
     std::string fft_device = PARAM.inp.device;
+
+    // LCAO basis doesn't support GPU acceleration on FFT currently
     if(PARAM.inp.basis_type == "lcao")
     {
         fft_device = "cpu";
     }
+
     pw_rho = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
-    if ( PARAM.globalv.double_grid)
+    if (PARAM.globalv.double_grid)
     {
         pw_rhod = new ModulePW::PW_Basis_Big(fft_device, PARAM.inp.precision);
     }
@@ -44,6 +45,7 @@ ESolver_FP::ESolver_FP()
     pw_big = static_cast<ModulePW::PW_Basis_Big*>(pw_rhod);
     pw_big->setbxyz(PARAM.inp.bx, PARAM.inp.by, PARAM.inp.bz);
     sf.set(pw_rhod, PARAM.inp.nbspline);
+
 }
 
 ESolver_FP::~ESolver_FP()
@@ -140,7 +142,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
     // 2) write fermi energy
     ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);
 
-    // 3) update delta rho for charge extrapolation
+    // 3) update delta_rho for charge extrapolation
     CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
 
     if (istep % PARAM.inp.out_interval == 0)
@@ -153,13 +155,13 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 double* data = nullptr;
                 if (PARAM.inp.dm_to_rho)
                 {
-                    data = this->pelec->charge->rho[is];
-                    this->pw_rhod->real2recip(this->pelec->charge->rho[is], this->pelec->charge->rhog[is]);
+                    data = this->chr.rho[is];
+                    this->pw_rhod->real2recip(this->chr.rho[is], this->chr.rhog[is]);
                 }
                 else
                 {
-                    data = this->pelec->charge->rho_save[is];
-                    this->pw_rhod->real2recip(this->pelec->charge->rho_save[is], this->pelec->charge->rhog_save[is]);
+                    data = this->chr.rho_save[is];
+                    this->pw_rhod->real2recip(this->chr.rho_save[is], this->chr.rhog_save[is]);
                 }
                 std::string fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_CHG.cube";
                 ModuleIO::write_vdata_palgrid(Pgrid,
@@ -176,7 +178,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 {
                     fn =PARAM.globalv.global_out_dir + "/SPIN" + std::to_string(is + 1) + "_TAU.cube";
                     ModuleIO::write_vdata_palgrid(Pgrid,
-                                                  this->pelec->charge->kin_r_save[is],
+                                                  this->chr.kin_r_save[is],
                                                   is,
                                                   PARAM.inp.nspin,
                                                   istep,
@@ -217,7 +219,7 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 fn,
                 istep,
                 this->pw_rhod,
-                this->pelec->charge,
+                &this->chr,
                 &(ucell),
                 this->pelec->pot->get_fixed_v(),
                 this->solvent);
@@ -226,11 +228,11 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
         // 6) write ELF
         if (PARAM.inp.out_elf[0] > 0)
         {
-            this->pelec->charge->cal_elf = true;
+            this->chr.cal_elf = true;
             Symmetry_rho srho;
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
-                srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
+                srho.begin(is, this->chr, this->pw_rhod, ucell.symm);
             }
 
             std::string out_dir =PARAM.globalv.global_out_dir;
@@ -242,8 +244,8 @@ void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_eso
                 out_dir,
                 istep,
                 PARAM.inp.nspin,
-                this->pelec->charge->rho,
-                this->pelec->charge->kin_r,
+                this->chr.rho,
+                this->chr.kin_r,
                 this->pw_rhod,
                 this->Pgrid,
                 &(ucell),
@@ -260,9 +262,9 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     {
         // only G-vector and K-vector are changed due to the change of lattice
         // vector FFT grids do not change!!
-        pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
-        pw_rho->collect_local_pw();
-        pw_rho->collect_uniqgg();
+        this->pw_rho->initgrids(ucell.lat0, ucell.latvec, pw_rho->nx, pw_rho->ny, pw_rho->nz);
+        this->pw_rho->collect_local_pw();
+        this->pw_rho->collect_uniqgg();
 
         if (PARAM.globalv.double_grid)
         {
@@ -292,7 +294,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->CE.update_all_dis(ucell);
         this->CE.extrapolate_charge(&(this->Pgrid),
                                     ucell,
-                                    this->pelec->charge,
+                                    &this->chr,
                                     &(this->sf),
                                     GlobalV::ofs_running,
                                     GlobalV::ofs_warning);
@@ -327,7 +329,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
             std::stringstream ss;
             ss << PARAM.globalv.global_out_dir << "SPIN" << is + 1 << "_CHG_INI.cube";
             ModuleIO::write_vdata_palgrid(this->Pgrid,
-                                          this->pelec->charge->rho[is],
+                                          this->chr.rho[is],
                                           is,
                                           PARAM.inp.nspin,
                                           istep,
@@ -368,8 +370,8 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
         if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
         {
             std::complex<double>** rhog_tot
-                = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rhog : this->pelec->charge->rhog_save;
-            double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rho : this->pelec->charge->rho_save;
+                = (PARAM.inp.dm_to_rho) ? this->chr.rhog : this->chr.rhog_save;
+            double** rhor_tot = (PARAM.inp.dm_to_rho) ? this->chr.rho : this->chr.rho_save;
             for (int is = 0; is < PARAM.inp.nspin; is++)
             {
                 this->pw_rhod->real2recip(rhor_tot[is], rhog_tot[is]);
@@ -386,12 +388,12 @@ void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool&
 
             if (XC_Functional::get_ked_flag())
             {
-                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->pelec->charge->ngmc, {0.0, 0.0});
+                std::vector<std::complex<double>> kin_g_space(PARAM.inp.nspin * this->chr.ngmc, {0.0, 0.0});
                 std::vector<std::complex<double>*> kin_g;
                 for (int is = 0; is < PARAM.inp.nspin; is++)
                 {
-                    kin_g.push_back(kin_g_space.data() + is * this->pelec->charge->ngmc);
-                    this->pw_rhod->real2recip(this->pelec->charge->kin_r_save[is], kin_g[is]);
+                    kin_g.push_back(kin_g_space.data() + is * this->chr.ngmc);
+                    this->pw_rhod->real2recip(this->chr.kin_r_save[is], kin_g[is]);
                 }
                 ModuleIO::write_rhog(PARAM.globalv.global_out_dir + PARAM.inp.suffix + "-TAU-DENSITY.restart",
                                      PARAM.globalv.gamma_only_pw || PARAM.globalv.gamma_only_local,

source/module_esolver/esolver_fp.h

@@ -2,16 +2,31 @@
 #define ESOLVER_FP_H
 
 #include "esolver.h"
+
+//! plane wave basis
 #include "module_basis/module_pw/pw_basis.h"
+
+//! symmetry analysis
 #include "module_cell/module_symmetry/symmetry.h"
+
+//! electronic states
 #include "module_elecstate/elecstate.h"
+
+//! charge extrapolation
 #include "module_elecstate/module_charge/charge_extra.h"
+
+//! solvation model
 #include "module_hamilt_general/module_surchem/surchem.h"
+
+//! local pseudopotential
 #include "module_hamilt_pw/hamilt_pwdft/VL_in_pw.h"
+
+//! structure factor related to plane wave basis
 #include "module_hamilt_pw/hamilt_pwdft/structure_factor.h"
 
 #include <fstream>
 
+
 //! The First-Principles (FP) Energy Solver Class
 /**
  * This class represents components that needed in
@@ -22,7 +37,7 @@
 
 namespace ModuleESolver
 {
-class ESolver_FP : public ESolver
+class ESolver_FP: public ESolver
 {
   public:
     //! Constructor
@@ -49,39 +64,34 @@ class ESolver_FP : public ESolver
     //! ------------------------------------------------------------------------------
     elecstate::ElecState* pelec = nullptr; ///< Electronic states
 
-    //! ------------------------------------------------------------------------------
+    //! K points in Brillouin zone
+    K_Vectors kv;
 
     //! Electorn charge density
     Charge chr;
 
-    //! Structure factors that used with plane-wave basis set
-    Structure_Factor sf;
-
-    //! K points in Brillouin zone
-    K_Vectors kv;
-
-    //! Plane-wave basis set for charge density
+    //! pw_rho: Plane-wave basis set for charge density
+    //! pw_rhod: same as pw_rho for NCPP. Here 'd' stands for 'dense',
+    //!          dense grid for for uspp, used for ultrasoft augmented charge density.
+    //!          charge density and potential are defined on dense grids,
+    //!          but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
     ModulePW::PW_Basis* pw_rho;
+    ModulePW::PW_Basis* pw_rhod;    //! dense grid for USPP
+    ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class
 
     //! parallel for rho grid
     Parallel_Grid Pgrid;
 
-    //! pointer to local pseudopotential
-    pseudopot_cell_vl locpp;
+    //! Structure factors that used with plane-wave basis set
+    Structure_Factor sf;
 
-    /**
-     * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense'
-     * dense grid for for uspp, used for ultrasoft augmented charge density.
-     * charge density and potential are defined on dense grids,
-     * but effective potential needs to be interpolated on smooth grids in order to compute Veff|psi>
-     */
-    ModulePW::PW_Basis* pw_rhod;
-    ModulePW::PW_Basis_Big* pw_big; ///< [temp] pw_basis_big class
+    //! local pseudopotentials
+    pseudopot_cell_vl locpp;
 
-    //! Charge extrapolation
+    //! charge extrapolation method
     Charge_Extra CE;
 
-    // solvent model
+    //! solvent model
     surchem solvent;
 };
 } // namespace ModuleESolver

source/module_esolver/esolver_ks.cpp

@@ -10,6 +10,7 @@
 #include "module_io/write_istate_info.h"
 #include "module_parameter/parameter.h"
 #include "module_elecstate/elecstate_print.h"
+#include "module_hsolver/hsolver.h"
 
 #include <ctime>
 #include <iostream>
@@ -367,7 +368,7 @@ void ESolver_KS<T, Device>::hamilt2density(UnitCell& ucell, const int istep, con
         // double drho = this->estate.caldr2();
         // EState should be used after it is constructed.
 
-        drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec);
+        drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
         hsolver_error = 0.0;
         if (iter == 1 && PARAM.inp.calculation != "nscf")
         {
@@ -389,7 +390,7 @@ void ESolver_KS<T, Device>::hamilt2density(UnitCell& ucell, const int istep, con
 
                 this->hamilt2density_single(ucell, istep, iter, diag_ethr);
 
-                drho = p_chgmix->get_drho(pelec->charge, PARAM.inp.nelec);
+                drho = p_chgmix->get_drho(&this->chr, PARAM.inp.nelec);
 
                 hsolver_error = hsolver::cal_hsolve_error(PARAM.inp.basis_type,
                                                           PARAM.inp.esolver_type,
@@ -521,7 +522,7 @@ void ESolver_KS<T, Device>::iter_init(UnitCell& ucell, const int istep, const in
     }
 
     // 1) save input rho
-    this->pelec->charge->save_rho_before_sum_band();
+    this->chr.save_rho_before_sum_band();
 }
 
 template <typename T, typename Device>
@@ -551,9 +552,9 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
 
     // compute magnetization, only for LSDA(spin==2)
     ucell.magnet.compute_magnetization(ucell.omega,
-                                       this->pelec->charge->nrxx,
-                                       this->pelec->charge->nxyz,
-                                       this->pelec->charge->rho,
+                                       this->chr.nrxx,
+                                       this->chr.nxyz,
+                                       this->chr.rho,
                                        this->pelec->nelec_spin.data());
 
     if (PARAM.globalv.ks_run)
@@ -628,11 +629,11 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
             }
             else
             {
-                p_chgmix->mix_rho(pelec->charge); // update chr->rho by mixing
+                p_chgmix->mix_rho(&this->chr); // update chr->rho by mixing
             }
             if (PARAM.inp.scf_thr_type == 2)
             {
-                pelec->charge->renormalize_rho(); // renormalize rho in R-space would
+                this->chr.renormalize_rho(); // renormalize rho in R-space would
                                                   // induce a error in K-space
             }
             //----------charge mixing done-----------
@@ -644,7 +645,7 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
 
     // be careful! conv_esolver is bool, not double !! Maybe a bug 20250302 by mohan 
     MPI_Bcast(&conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD);
-    MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD);
+    MPI_Bcast(this->chr.rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD);
 #endif
 
     // update potential
@@ -676,7 +677,7 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
     double dkin = 0.0; // for meta-GGA
     if (XC_Functional::get_ked_flag())
     {
-        dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec);
+        dkin = p_chgmix->get_dkin(&this->chr, PARAM.inp.nelec);
     }