Update some formats in LCAO_operator and ESolver (deepmodeling#5987)

* small format updates in operator_lcao

* update esolver information

Author: mohanchen

source/module_esolver/esolver_ks_lcao.cpp

@@ -67,14 +67,12 @@
 namespace ModuleESolver
 {
 
-//------------------------------------------------------------------------------
-//! the 1st function of ESolver_KS_LCAO: constructor
-//------------------------------------------------------------------------------
 template <typename TK, typename TR>
 ESolver_KS_LCAO<TK, TR>::ESolver_KS_LCAO()
 {
     this->classname = "ESolver_KS_LCAO";
     this->basisname = "LCAO";
+
 #ifdef __EXX
     // 1. currently this initialization must be put in constructor rather than `before_all_runners()`
     //  because the latter is not reused by ESolver_LCAO_TDDFT,
@@ -94,47 +92,30 @@ ESolver_KS_LCAO<TK, TR>::ESolver_KS_LCAO()
 #endif
 }
 
-//------------------------------------------------------------------------------
-//! the 2nd function of ESolver_KS_LCAO: deconstructor
-//------------------------------------------------------------------------------
 template <typename TK, typename TR>
 ESolver_KS_LCAO<TK, TR>::~ESolver_KS_LCAO()
 {
 }
 
-//------------------------------------------------------------------------------
-//! the 3rd function of ESolver_KS_LCAO: init
-//! 1) calculate overlap matrix S or initialize
-//! 2) init ElecState
-//! 3) init LCAO basis
-//! 4) initialize the density matrix
-//! 5) initialize exx
-//! 6) initialize DFT+U
-//! 7) ppcell
-//! 8) inititlize the charge density
-//! 9) initialize the potential.
-//! 10) initialize deepks
-//! 11) set occupations
-//! 12) print a warning if needed
-//------------------------------------------------------------------------------
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_para& inp)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "before_all_runners");
     ModuleBase::timer::tick("ESolver_KS_LCAO", "before_all_runners");
 
+    // 1) before_all_runners in ESolver_KS
     ESolver_KS<TK>::before_all_runners(ucell, inp);
 
     // 2) init ElecState
-    // autoset nbands in ElecState, it should before basis_init (for Psi 2d division)
+    // autoset nbands in ElecState before basis_init (for Psi 2d division)
     if (this->pelec == nullptr)
     {
         // TK stands for double and complex<double>?
         this->pelec = new elecstate::ElecStateLCAO<TK>(&(this->chr), // use which parameter?
                                                        &(this->kv),
                                                        this->kv.get_nks(),
-                                                       &(this->GG), // mohan add 2024-04-01
-                                                       &(this->GK), // mohan add 2024-04-01
+                                                       &(this->GG),
+                                                       &(this->GK),
                                                        this->pw_rho,
                                                        this->pw_big);
     }
@@ -157,9 +138,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     // DMR is not initialized here, it will be constructed in each before_scf
     dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->init_DM(&this->kv, &(this->pv), PARAM.inp.nspin);
 
+    // 5) initialize exact exchange calculations
 #ifdef __EXX
-    // 5) initialize exx
-    // PLEASE simplify the Exx_Global interface
     if (PARAM.inp.calculation == "scf" 
         || PARAM.inp.calculation == "relax" 
         || PARAM.inp.calculation == "cell-relax"
@@ -190,7 +170,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         dftu->init(ucell, &this->pv, this->kv.get_nks(), &orb_);
     }
 
-    // 7) initialize ppcell
+    // 7) initialize local pseudopotentials
     this->locpp.init_vloc(ucell, this->pw_rho);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
@@ -211,8 +191,9 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
                                                     &(this->pelec->f_en.vtxc));
     }
 
-#ifdef __DEEPKS
+
     // 10) initialize deepks
+#ifdef __DEEPKS
     LCAO_domain::DeePKS_init(ucell, pv, this->kv.get_nks(), orb_, this->ld, GlobalV::ofs_running);
     if (PARAM.inp.deepks_scf)
     {
@@ -279,10 +260,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     return;
 }
 
-//------------------------------------------------------------------------------
-//! the 5th function of ESolver_KS_LCAO: cal_energy
-//! mohan add 2024-05-11
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 double ESolver_KS_LCAO<TK, TR>::cal_energy()
 {
@@ -291,10 +269,7 @@ double ESolver_KS_LCAO<TK, TR>::cal_energy()
     return this->pelec->f_en.etot;
 }
 
-//------------------------------------------------------------------------------
-//! the 6th function of ESolver_KS_LCAO: cal_force
-//! mohan add 2024-05-11
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& force)
 {
@@ -380,6 +355,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
     GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl;
 
+    // 1) write information
     if (PARAM.inp.out_dos != 0 || PARAM.inp.out_band[0] != 0 || PARAM.inp.out_proj_band != 0)
     {
         GlobalV::ofs_running << "\n\n\n\n";
@@ -412,14 +388,16 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                              << std::endl;
         GlobalV::ofs_running << "\n\n\n\n";
     }
-    // qianrui modify 2020-10-18
+
+    // 2) write information
     if (PARAM.inp.calculation == "scf" || PARAM.inp.calculation == "md" || PARAM.inp.calculation == "relax")
     {
         ModuleIO::write_istate_info(this->pelec->ekb, this->pelec->wg, this->kv);
     }
 
     const int nspin0 = (PARAM.inp.nspin == 2) ? 2 : 1;
 
+    // 3) print out band information
     if (PARAM.inp.out_band[0])
     {
         for (int is = 0; is < nspin0; is++)
@@ -435,13 +413,15 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                 this->pelec->ekb,
                                 this->kv);
         }
-    } // out_band
+    } 
 
-    if (PARAM.inp.out_proj_band) // Projeced band structure added by jiyy-2022-4-20
+    // 4) write projected band structure by jiyy-2022-4-20
+    if (PARAM.inp.out_proj_band)
     {
         ModuleIO::write_proj_band_lcao(this->psi, this->pv, this->pelec, this->kv, ucell, this->p_hamilt);
     }
 
+    // 5) print out density of states (DOS)
     if (PARAM.inp.out_dos)
     {
         ModuleIO::out_dos_nao(this->psi,
@@ -458,6 +438,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                               this->p_hamilt);
     }
 
+    // 6) print out exchange-correlation potential
     if (PARAM.inp.out_mat_xc)
     {
         ModuleIO::write_Vxc<TK, TR>(PARAM.inp.nspin,
@@ -486,6 +467,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
         );
     }
 
+    // 7) write eband terms
     if (PARAM.inp.out_eband_terms)
     {
         ModuleIO::write_eband_terms<TK, TR>(PARAM.inp.nspin,
@@ -518,10 +500,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
     ModuleBase::timer::tick("ESolver_KS_LCAO", "after_all_runners");
 }
 
-//------------------------------------------------------------------------------
-//! the 10th function of ESolver_KS_LCAO: iter_init
-//! mohan add 2024-05-11
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const int iter)
 {
@@ -546,6 +525,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
             std::cout << " eV " << std::endl;
         }
     }
+
     // for mixing restart
     if (iter == this->p_chgmix->mixing_restart_step && PARAM.inp.mixing_restart > 0.0)
     {
@@ -581,7 +561,6 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     // mohan update 2012-06-05
     this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
 
-    // mohan move it outside 2011-01-13
     // first need to calculate the weight according to
     // electrons number.
     if (istep == 0 && PARAM.inp.init_wfc == "file")
@@ -591,7 +570,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
             std::cout << " WAVEFUN -> CHARGE " << std::endl;
 
             // calculate the density matrix using read in wave functions
-            // and the ncalculate the charge density on grid.
+            // and then calculate the charge density on grid.
 
             this->pelec->skip_weights = true;
             this->pelec->calculate_weights();
@@ -695,33 +674,18 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     }
 }
 
-//------------------------------------------------------------------------------
-//! the 11th function of ESolver_KS_LCAO: hamilt2density_single
-//! mohan add 2024-05-11
-//! 1) save input rho
-//! 2) save density matrix DMR for mixing
-//! 3) solve the Hamiltonian and output band gap
-//! 4) print bands for each k-point and each band
-//! 5) EXX:
-//! 6) DFT+U: compute local occupation number matrix and energy correction
-//! 7) DeePKS: compute delta_e
-//! 8) DeltaSpin:
-//! 9) use new charge density to calculate energy
-//! 10) symmetrize the charge density
-//! 11) compute magnetization, only for spin==2
-//! 12) calculate delta energy
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(UnitCell& ucell, int istep, int iter, double ethr)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density_single");
 
-    // reset energy
+    // i1) reset energy
     this->pelec->f_en.eband = 0.0;
     this->pelec->f_en.demet = 0.0;
     bool skip_charge = PARAM.inp.calculation == "nscf" ? true : false;
 
-    // run the inner lambda loop to contrain atomic moments with the DeltaSpin method
+    // 2) run the inner lambda loop to contrain atomic moments with the DeltaSpin method
     bool skip_solve = false;
     if (PARAM.inp.sc_mag_switch)
     {
@@ -740,13 +704,15 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(UnitCell& ucell, int istep,
             skip_solve = true;
         }
     }
+
+    // 3) run Hsolver
     if (!skip_solve)
     {
         hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
         hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, skip_charge);
     }
 
-    // 5) what's the exd used for?
+    // 4) EXX
 #ifdef __EXX
     if (PARAM.inp.calculation != "nscf")
     {
@@ -761,22 +727,18 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(UnitCell& ucell, int istep,
     }
 #endif
 
-    // 10) symmetrize the charge density
+    // 5) symmetrize the charge density
     Symmetry_rho srho;
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         srho.begin(is, this->chr, this->pw_rho, ucell.symm);
     }
 
-    // 12) calculate delta energy
+    // 6) calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
 }
 
-//------------------------------------------------------------------------------
-//! the 12th function of ESolver_KS_LCAO: update_pot
-//! mohan add 2024-05-11
-//! 1) print potential
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver)
 {
@@ -794,21 +756,14 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(UnitCell& ucell, const int istep, const
     }
 }
 
-//------------------------------------------------------------------------------
-//! the 13th function of ESolver_KS_LCAO: iter_finish
-//! mohan add 2024-05-11
-//! 1) mix density matrix
-//! 2) output charge density
-//! 3) output exx matrix
-//! 4) output charge density and density matrix
-//------------------------------------------------------------------------------
+
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
 
-    // 6) calculate the local occupation number matrix and energy correction in
-    // DFT+U
+    // 1) calculate the local occupation number matrix and energy correction
+    // in DFT+U
     if (PARAM.inp.dft_plus_u)
     {
         // only old DFT+U method should calculated energy correction in esolver,
@@ -831,7 +786,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
         GlobalC::dftu.output(ucell);
     }
 
-    // (7) for deepks, calculate delta_e
+    // 2) for deepks, calculate delta_e
 #ifdef __DEEPKS
     if (PARAM.inp.deepks_scf)
     {
@@ -844,25 +799,25 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 #endif
 
-    // 8) for delta spin
+    // 3) for delta spin
     if (PARAM.inp.sc_mag_switch)
     {
         spinconstrain::SpinConstrain<TK>& sc = spinconstrain::SpinConstrain<TK>::getScInstance();
         sc.cal_mi_lcao(iter);
     }
 
-    // call iter_finish() of ESolver_KS
+    // 4) call iter_finish() of ESolver_KS
     ESolver_KS<TK>::iter_finish(ucell, istep, iter, conv_esolver);
 
-    // 1) mix density matrix if mixing_restart + mixing_dmr + not first
+    // 5) mix density matrix if mixing_restart + mixing_dmr + not first
     // mixing_restart at every iter
     if (PARAM.inp.mixing_restart > 0 && this->p_chgmix->mixing_restart_count > 0 && PARAM.inp.mixing_dmr)
     {
         elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
         this->p_chgmix->mix_dmr(dm);
     }
 
-    // 2) save charge density
+    // 6) save charge density
     // Peize Lin add 2020.04.04
     if (GlobalC::restart.info_save.save_charge)
     {
@@ -873,7 +828,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 
 #ifdef __EXX
-    // 3) save exx matrix
+    // 7) save exx matrix
     if (PARAM.inp.calculation != "nscf")
     {
         if (GlobalC::exx_info.info_global.cal_exx)
@@ -900,7 +855,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(UnitCell& ucell, const int istep, int&
     }
 #endif
 
-    // 6) use the converged occupation matrix for next MD/Relax SCF calculation
+    // 8) use the converged occupation matrix for next MD/Relax SCF calculation
     if (PARAM.inp.dft_plus_u && conv_esolver)
     {
         GlobalC::dftu.initialed_locale = true;