Clean and Reconstruct ESolver --- 3rd time (#6567)

* add setup_dm in module_dm

* update before_all_runners in esolver

* change setup_dm back

* add setup_exx

* fix bug

* add setup_exx in normal LCAO codes

* fix bug about p_hamilt and fix setup_exx

* fix bug

* fix bug

* eliminate exd and exc from ctrl_iter_lcao input parameters

* add setup_deepks

* update the input variables of ctrl_scf_lcao

* reconstruct the constructor of p_hamilt

* update ld

* fix a few bugs

* DeePKS has been used in operator, fix it

Author: mohanchen

examples/relax/lcao_output/INPUT

@@ -40,4 +40,3 @@ out_band        0
 out_stru        0
 out_app_flag    0
 
-out_interval    1

source/Makefile.Objects

@@ -253,6 +253,7 @@ OBJS_ELECSTAT=elecstate.o\
 
 OBJS_ELECSTAT_LCAO=elecstate_lcao.o\
       elecstate_lcao_cal_tau.o\
+      setup_dm.o\
       density_matrix.o\
       density_matrix_io.o\
       cal_dm_psi.o\
@@ -668,6 +669,8 @@ OBJS_LCAO=evolve_elec.o\
       LCAO_allocate.o\
       LCAO_set_mat2d.o\
       LCAO_init_basis.o\
+      setup_exx.o\
+      setup_deepks.o\
       center2_orb.o\
       center2_orb-orb11.o\
       center2_orb-orb21.o\
@@ -676,13 +679,13 @@ OBJS_LCAO=evolve_elec.o\
       wavefunc_in_pw.o\
 
 OBJS_MODULE_RI=conv_coulomb_pot_k.o\
-      exx_abfs-abfs_index.o \
-      exx_abfs-jle.o \
-      exx_abfs-io.o \
-      exx_abfs-construct_orbs.o \
-      ABFs_Construct-PCA.o \
-      exx_opt_orb.o \
-      exx_opt_orb-print.o \
+      exx_abfs-abfs_index.o\
+      exx_abfs-jle.o\
+      exx_abfs-io.o\
+      exx_abfs-construct_orbs.o\
+      ABFs_Construct-PCA.o\
+      exx_opt_orb-print.o\
+      exx_opt_orb.o\
       Matrix_Orbs11.o\
       Matrix_Orbs21.o\
       Matrix_Orbs22.o\

source/source_esolver/esolver_double_xc.cpp

@@ -154,17 +154,14 @@ void ESolver_DoubleXC<TK, TR>::before_scf(UnitCell& ucell, const int istep)
             this->kv,
             this->two_center_bundle_,
             this->orb_,
-            DM
-#ifdef __MLALGO
-            ,
-            &this->ld
-#endif
+            DM,
+            this->deepks
 #ifdef __EXX
             ,
             istep,
-            GlobalC::exx_info.info_ri.real_number ? &this->exd->two_level_step : &this->exc->two_level_step,
-            GlobalC::exx_info.info_ri.real_number ? &this->exd->get_Hexxs() : nullptr,
-            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exc->get_Hexxs()
+            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->two_level_step : &this->exx_nao.exc->two_level_step,
+            GlobalC::exx_info.info_ri.real_number ? &this->exx_nao.exd->get_Hexxs() : nullptr,
+            GlobalC::exx_info.info_ri.real_number ? nullptr : &this->exx_nao.exc->get_Hexxs()
 #endif
         );
     }
@@ -238,7 +235,7 @@ void ESolver_DoubleXC<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int
 #ifdef __MLALGO
         // ---------- output tot and precalc ----------
         hamilt::HamiltLCAO<TK, TR>* p_ham_deepks = dynamic_cast<hamilt::HamiltLCAO<TK, TR>*>(this->p_hamilt);
-        std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->ld, [](LCAO_Deepks<TK>*) {});
+        std::shared_ptr<LCAO_Deepks<TK>> ld_shared_ptr(&this->deepks.ld, [](LCAO_Deepks<TK>*) {});
         LCAO_Deepks_Interface<TK, TR> deepks_interface(ld_shared_ptr);
 
         deepks_interface.out_deepks_labels(this->pelec->f_en.etot,
@@ -412,14 +409,11 @@ void ESolver_DoubleXC<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                        this->pw_rho,
                        this->solvent,
 #ifdef __MLALGO
-                       this->ld,
+                       this->deepks.ld,
                        "base",
 #endif
-#ifdef __EXX
-                       *this->exd,
-                       *this->exc,
-#endif
-                       &ucell.symm);
+					   this->exx_nao,
+					   &ucell.symm);
     // restore to original xc
     XC_Functional::set_xc_type(ucell.atoms[0].ncpp.xc_func); 
 

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

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);
 }