update before_all_runners in esolver

mohanchen · mohanchen · commit 74b3500ed1ca · 2025-10-08T11:47:44.000+08:00
diff --git a/source/source_esolver/esolver_fp.cpp b/source/source_esolver/esolver_fp.cpp
@@ -48,21 +48,21 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
       this->pw_rho, this->pw_rhod, this->pw_big, 
       this->classname, inp);
 
-    // setup the structure factors
+    // setup structure factors
     this->sf.set(this->pw_rhod, inp.nbspline);
 
+    // write geometry file
     ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif",
                                ucell,
                                "# Generated by ABACUS ModuleIO::CifParser",
                                "data_?");
 
-    //! initialize the charge extrapolation method if necessary
+    // init charge extrapolation
     this->CE.Init_CE(inp.nspin, ucell.nat, this->pw_rhod->nrxx, inp.chg_extrap);
 
     return;
 }
 
-//! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
 void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)
 {
     ModuleBase::TITLE("ESolver_FP", "after_scf");
@@ -125,12 +125,8 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     if (ucell.ionic_position_updated)
     {
         this->CE.update_all_dis(ucell);
-        this->CE.extrapolate_charge(&this->Pgrid,
-                                    ucell,
-                                    &this->chr,
-                                    &this->sf,
-                                    GlobalV::ofs_running,
-                                    GlobalV::ofs_warning);
+        this->CE.extrapolate_charge(&this->Pgrid, ucell, &this->chr, &this->sf,
+                                    GlobalV::ofs_running, GlobalV::ofs_warning);
     }
 
     //! calculate D2 or D3 vdW
diff --git a/source/source_esolver/esolver_ks.cpp b/source/source_esolver/esolver_ks.cpp
@@ -51,9 +51,11 @@ template <typename T, typename Device>
 void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
     ModuleBase::TITLE("ESolver_KS", "before_all_runners");
-    //! 1) initialize "before_all_runniers" in ESolver_FP
+
+    //! 1) init "before_all_runniers" in ESolver_FP
     ESolver_FP::before_all_runners(ucell, inp);
     
+    //! 2) setup some parameters
     classname = "ESolver_KS";
     basisname = "";
 
@@ -65,8 +67,8 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
 
     std::string fft_device = inp.device;
 
-    // Fast Fourier Transform
-    // LCAO basis doesn't support GPU acceleration on FFT currently
+    //! 3) setup pw_wfc
+    // currently LCAO doesn't support GPU acceleration of FFT
     if(inp.basis_type == "lcao")
     {
         fft_device = "cpu";
@@ -82,71 +84,49 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
     pw_wfc = new ModulePW::PW_Basis_K_Big(fft_device, fft_precision);
     ModulePW::PW_Basis_K_Big* tmp = static_cast<ModulePW::PW_Basis_K_Big*>(pw_wfc);
 
-    // should not use INPUT here, mohan 2024-05-12
     tmp->setbxyz(inp.bx, inp.by, inp.bz);
 
-    ///----------------------------------------------------------
-    /// charge mixing
-    ///----------------------------------------------------------
+    //! 4) setup charge mixing
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
 
     // cell_factor
     this->ppcell.cell_factor = inp.cell_factor;
 
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
 
-    //! 3) it has been established that
-    // xc_func is same for all elements, therefore
-    // only the first one if used
+    //! 5) setup Exc for the first element '0', because all elements have same exc 
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func);
     
-    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "SETUP UNITCELL");
-
-    //! 4) 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);
+    //! 6) 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();
 
-    //! 5) ESolver depends on the Symmetry module
-    // symmetry analysis should be performed every time the cell is changed
+    //! 7) 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");
     }
 
-    //! 6) Setup the k points according to symmetry.
+    //! 8) Setup the k points 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");
 
-    //! 7) print information
+    //! 9) print information
     ModuleIO::setup_parameters(ucell, this->kv);
 
-    //! 8) setup plane wave for electronic wave functions
+    //! 10) setup plane wave for electronic wave functions
     ModuleESolver::pw_setup(inp, ucell, *this->pw_rho, this->kv, *this->pw_wfc);
 
-    //! 9) initialize the real-space uniform grid for FFT and parallel
-    //! distribution of plane waves
-	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
+    //! 11) 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);
+
+    //! 12) calculate the structure factor
     this->sf.setup_structure_factor(&ucell, Pgrid, this->pw_rhod);
 }
 
diff --git a/source/source_esolver/esolver_ks_lcao.cpp b/source/source_esolver/esolver_ks_lcao.cpp
@@ -1,53 +1,24 @@
 #include "esolver_ks_lcao.h"
-
-#include "source_base/global_variable.h"
-#include "source_base/tool_title.h"
 #include "source_estate/elecstate_tools.h"
-
-#include "source_estate/module_dm/cal_dm_psi.h"
 #include "source_lcao/module_deltaspin/spin_constrain.h"
-#include "source_lcao/module_dftu/dftu.h"
-#include "source_io/cube_io.h"
-#include "source_io/output_log.h"
 #include "source_io/read_wfc_nao.h"
-#include "source_io/write_elecstat_pot.h"
-#include "source_io/module_parameter/parameter.h"
-
-// be careful of hpp, there may be multiple definitions of functions, 20250302, mohan
-#include "source_lcao/hs_matrix_k.hpp"
-
-#include "source_base/global_function.h"
-#include "source_cell/module_neighbor/sltk_grid_driver.h"
+#include "source_lcao/hs_matrix_k.hpp" // there may be multiple definitions if using hpp
 #include "source_estate/cal_ux.h"
 #include "source_estate/module_charge/symmetry_rho.h"
-#include "source_estate/occupy.h"
 #include "source_lcao/LCAO_domain.h" // need DeePKS_init
 #include "source_lcao/module_dftu/dftu.h"
-#include "source_pw/module_pwdft/global.h"
-#include "source_io/print_info.h"
-
-#include <memory>
-
 #ifdef __MLALGO
 #include "source_lcao/module_deepks/LCAO_deepks.h"
 #include "source_lcao/module_deepks/LCAO_deepks_interface.h"
 #endif
-//-----force& stress-------------------
 #include "source_lcao/FORCE_STRESS.h"
-
-//-----HSolver ElecState Hamilt--------
 #include "source_estate/elecstate_lcao.h"
 #include "source_lcao/hamilt_lcao.h"
 #include "source_hsolver/hsolver_lcao.h"
-
 #ifdef __EXX
 #include "../source_lcao/module_ri/exx_opt_orb.h"
 #endif
-
-// test RDMFT
 #include "source_lcao/module_rdmft/rdmft.h"
-#include "source_lcao/module_gint/temp_gint/gint_info.h"
-
 #include "source_estate/module_charge/chgmixing.h" // use charge mixing, mohan add 20251006 
 #include "source_estate/module_dm/setup_dm.h" // setup dm from electronic wave functions
 #include "source_io/ctrl_runner_lcao.h" // use ctrl_runner_lcao() 
@@ -96,32 +67,31 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     // 1) before_all_runners in ESolver_KS
     ESolver_KS<TK>::before_all_runners(ucell, inp);
 
-    // 2) init ElecState
-    // autoset nbands in ElecState before init_basis (for Psi 2d division)
+    // 2) autoset nbands in ElecState before init_basis (for Psi 2d division)
     if (this->pelec == nullptr)
     {
         // TK stands for double and std::complex<double>?
         this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), &(this->kv),
           this->kv.get_nks(), &(this->GG), &(this->GK), this->pw_rho, this->pw_big);
     }
 
-    // 3) read the LCAO orbitals/projectors and construct the interpolation tables.
+    // 3) read LCAO orbitals/projectors and construct the interpolation tables.
     LCAO_domain::init_basis_lcao(this->pv, inp.onsite_radius, inp.lcao_ecut,
       inp.lcao_dk, inp.lcao_dr, inp.lcao_rmax, ucell, two_center_bundle_, orb_);
 
     // 4) setup EXX calculations
     if (PARAM.inp.calculation == "gen_opt_abfs")
     {
-      #ifdef __EXX
+#ifdef __EXX
         Exx_Opt_Orb exx_opt_orb;
         exx_opt_orb.generate_matrix(GlobalC::exx_info.info_opt_abfs, this->kv, ucell, this->orb_);
-      #else
+#else
         ModuleBase::WARNING_QUIT("ESolver_KS_LCAO::before_all_runners", "calculation=gen_opt_abfs must compile __EXX");
-      #endif
+#endif
         return;
     }
 
-    // 5) initialize electronic wave function psi
+    // 5) init electronic wave function psi
     if (this->psi == nullptr)
     {
         int nsk = 0;
@@ -160,11 +130,10 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         }
     }
 
-    // 7) initialize the density matrix
-    // DMK are allocated here, but DMR is constructed in before_scf()
+    // 7) init DMK, but DMR is constructed in before_scf()
     dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), inp.nspin);
 
-    // 8) initialize exact exchange calculations
+    // 8) init exact exchange calculations
 #ifdef __EXX
     if (inp.calculation == "scf" || inp.calculation == "relax" || inp.calculation == "cell-relax"
         || inp.calculation == "md")
@@ -198,23 +167,23 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
     }
 
-    // 10) initialize local pseudopotentials
+    // 10) init local pseudopotentials
     this->locpp.init_vloc(ucell, this->pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
-    // 11) inititlize the charge density
+    // 11) init charge density
     this->chr.allocate(inp.nspin);
     this->pelec->omega = ucell.omega;
 
-    // 12) initialize the potential
+    // 12) init potentials
     if (this->pelec->pot == nullptr)
     {
         this->pelec->pot = new elecstate::Potential(this->pw_rhod, this->pw_rho,
           &ucell, &(this->locpp.vloc), &(this->sf), &(this->solvent),
           &(this->pelec->f_en.etxc), &(this->pelec->f_en.vtxc));
     }
 
-    // 13) initialize deepks
+    // 13) init deepks
 #ifdef __MLALGO
     LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld, GlobalV::ofs_running);
     if (inp.deepks_scf)
@@ -228,8 +197,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 #endif
 
-    // 14) set occupations
-    // tddft does not need to set occupations in the first scf
+    // 14) set occupations, tddft does not need to set occupations in the first scf
     if (inp.ocp && inp.esolver_type != "tddft")
     {
         elecstate::fixed_weights(inp.ocp_kb, inp.nbands, inp.nelec,
@@ -257,7 +225,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         }
     }
 
-    // 16) initialize rdmft, added by jghan
+    // 16) init rdmft, added by jghan
     if (inp.rdmft == true)
     {
         rdmft_solver.init(this->GG, this->GK, this->pv, ucell,
diff --git a/source/source_esolver/lcao_before_scf.cpp b/source/source_esolver/lcao_before_scf.cpp
