Refactoring LCAO codes.

source/module_esolver/esolver_ks.cpp

@@ -392,7 +392,9 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 		ModuleBase::timer::tick(this->classname, "run");
 
 		this->before_scf(istep); //Something else to do before the iter loop
+
 		ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
+
 		if(this->maxniter > 0)  
 		{
 			this->print_head(); //print the headline on the screen.
@@ -499,9 +501,6 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
                (std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::system_clock::now() 
                 - iterstart)).count() / static_cast<double>(1e6);
 #endif
-			/*
-			   SCF print: G1    -3.435545e+03  0.000000e+00   3.607e-01  2.862e-01
-			 */
 
 			double dkin = 0.0; // for meta-GGA
 			if (XC_Functional::get_func_type() == 3 || XC_Functional::get_func_type() == 5)
@@ -525,7 +524,10 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 			{
 				this->niter = iter;
 				bool stop = this->do_after_converge(iter);
-				if(stop) break;
+				if(stop) 
+				{
+					break;
+				}
 			}
 
 			// notice for restart
@@ -544,7 +546,9 @@ void ESolver_KS<T, Device>::run(const int istep, UnitCell& ucell)
 				this->conv_elec
 				);
 #endif //__RAPIDJSON 
-		after_scf(istep);
+
+		this->after_scf(istep);
+
 		ModuleBase::timer::tick(this->classname, "run");
 	} 
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -1,4 +1,5 @@
 #include "esolver_ks_lcao.h"
+#include "module_base/tool_title.h"
 
 #include "module_base/global_variable.h"
 #include "module_io/dos_nao.h"
@@ -21,7 +22,6 @@
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
 #include "module_io/print_info.h"
 #ifdef __EXX
-// #include "module_rpa/rpa.h"
 #include "module_ri/RPA_LRI.h"
 #endif
 
@@ -81,6 +81,8 @@ template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "init");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
+
     // if we are only calculating S, then there is no need
     // to prepare for potentials and so on
 
@@ -123,7 +125,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     // Init Basis should be put outside of Ensolver.
     // * reading the localized orbitals/projectors
     // * construct the interpolation tables.
-    this->Init_Basis_lcao(this->orb_con, inp, ucell);
+    this->init_basis_lcao(this->orb_con, inp, ucell);
     //------------------init Basis_lcao----------------------
 
     //! pass Hamilt-pointer to Operator
@@ -137,6 +139,7 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
 
     if (GlobalV::CALCULATION == "get_S")
     {
+        ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
         return;
     }
 
@@ -230,13 +233,19 @@ void ESolver_KS_LCAO<TK, TR>::init(Input& inp, UnitCell& ucell)
     if (GlobalV::ocp)
     {
         this->pelec->fixed_weights(GlobalV::ocp_kb);
-    }
+	}
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init");
+	return;
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "init_after_vc");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc");
+
 	ESolver_KS<TK>::init_after_vc(inp, ucell);
 
     if (GlobalV::md_prec_level == 2)
@@ -269,26 +278,38 @@ void ESolver_KS_LCAO<TK, TR>::init_after_vc(Input& inp, UnitCell& ucell)
                                                         &(this->pelec->f_en.etxc),
                                                         &(this->pelec->f_en.vtxc));
         }
-    }
+	}
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "init_after_vc");
+	return;
 }
 
 
 template <typename TK, typename TR>
 double ESolver_KS_LCAO<TK, TR>::cal_energy()
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_energy");
+
     return this->pelec->f_en.etot;
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 {
-	Force_Stress_LCAO<TK> FSL(this->RA, GlobalC::ucell.nat);
-	FSL.getForceStress(GlobalV::CAL_FORCE,
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_force");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force");
+
+	Force_Stress_LCAO<TK> fsl(this->RA, GlobalC::ucell.nat);
+
+	fsl.getForceStress(
+            GlobalV::CAL_FORCE,
 			GlobalV::CAL_STRESS,
 			GlobalV::TEST_FORCE,
 			GlobalV::TEST_STRESS,
 			this->LOC,
+            this->orb_con.ParaV, 
+            this->LM,
 			this->pelec,
 			this->psi,
 			this->UHM,
@@ -308,25 +329,35 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(ModuleBase::matrix& force)
 	this->RA.delete_grid();
 
 	this->have_force = true;
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_force");
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::cal_stress(ModuleBase::matrix& stress)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "cal_stress");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress");
+
     if (!this->have_force)
     {
         ModuleBase::matrix fcs;
         this->cal_force(fcs);
     }
     stress = this->scs; // copy the stress
     this->have_force = false;
+
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "cal_stress");
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::post_process(void)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "post_process");
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process");
+
     GlobalV::ofs_running << "\n\n --------------------------------------------" << std::endl;
     GlobalV::ofs_running << std::setprecision(16);
     GlobalV::ofs_running << " !FINAL_ETOT_IS " << this->pelec->f_en.etot * ModuleBase::Ry_to_eV << " eV" << std::endl;
@@ -413,15 +444,18 @@ void ESolver_KS_LCAO<TK, TR>::post_process(void)
 				GlobalV::NBANDS,
 				this->p_hamilt);
 	}
+    ModuleBase::timer::tick("ESolver_KS_LCAO", "post_process");
 }
 
 
 template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::Init_Basis_lcao(
+void ESolver_KS_LCAO<TK, TR>::init_basis_lcao(
 		ORB_control& orb_con, 
 		Input& inp, 
 		UnitCell& ucell)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "init_basis_lcao");
+
     // autoset NB2D first
     if (GlobalV::NB2D == 0)
     {
@@ -500,12 +534,16 @@ void ESolver_KS_LCAO<TK, TR>::Init_Basis_lcao(
         this->orb_con.setup_2d_division(GlobalV::ofs_running, GlobalV::ofs_warning);
         this->orb_con.ParaV.set_atomic_trace(GlobalC::ucell.get_iat2iwt(), GlobalC::ucell.nat, GlobalV::NLOCAL);
     }
+
+    return;
 }
 
 
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_init");
+
     if (iter == 1)
     {
         this->p_chgmix->init_mixing(); // init mixing
@@ -645,6 +683,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density");
+
     // save input rho
     this->pelec->charge->save_rho_before_sum_band();
     // save density matrix for mixing
@@ -756,6 +796,8 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "update_pot");
+
     // print Hamiltonian and Overlap matrix
     if (this->conv_elec)
     {
@@ -852,6 +894,8 @@ void ESolver_KS_LCAO<TK, TR>::update_pot(const int istep, const int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::iter_finish(int iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
+
     // mix density matrix
     if (GlobalV::MIXING_RESTART > 0 && iter >= this->p_chgmix->mixing_restart && GlobalV::MIXING_DMR )
     {
@@ -930,6 +974,8 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(int iter)
 template <typename TK, typename TR>
 void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO", "after_scf");
+
     // save charge difference into files for charge extrapolation
     if (GlobalV::CALCULATION != "scf")
     {
@@ -1015,6 +1061,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
             this->LOWF.ParaV,
             *(this->psi),
             dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM());
+
     ModuleBase::timer::tick("ESolver_KS_LCAO", "out_deepks_labels");
 
 #endif
@@ -1042,6 +1089,7 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
         } // qifeng add 2019/9/10, jiyy modify 2023/2/27, liuyu move here 2023-04-18
     }
 
+    // spin constrain calculations, added by Tianqi Zhao.
     if (GlobalV::sc_mag_switch)
     {
         SpinConstrain<TK, psi::DEVICE_CPU>& sc = SpinConstrain<TK, psi::DEVICE_CPU>::getScInstance();
@@ -1052,6 +1100,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
     {
         RA.delete_grid();
     }
+
+    // quasi-orbitals, added by Yike Huang.
     if(GlobalV::qo_switch)
     {
         toQO tqo(GlobalV::qo_basis, GlobalV::qo_strategy, GlobalV::qo_thr, GlobalV::qo_screening_coeff);
@@ -1069,6 +1119,8 @@ void ESolver_KS_LCAO<TK, TR>::after_scf(const int istep)
 template <typename TK, typename TR>
 bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter)
 {
+    ModuleBase::TITLE("ESolver_KS_LCAO","do_after_converge");
+
 #ifdef __EXX
     if (GlobalC::exx_info.info_ri.real_number)
     {
@@ -1088,11 +1140,13 @@ bool ESolver_KS_LCAO<TK, TR>::do_after_converge(int& iter)
 				iter);
 	}
 #endif // __EXX
+
     if(GlobalV::dft_plus_u)
     {
         // use the converged occupation matrix for next MD/Relax SCF calculation
         GlobalC::dftu.initialed_locale = true;
     }
+
     return true;
 }
 
@@ -1127,17 +1181,17 @@ ModuleIO::Output_DM1 ESolver_KS_LCAO<TK, TR>::create_Output_DM1(int istep)
 template <typename TK, typename TR>
 ModuleIO::Output_Mat_Sparse<TK> ESolver_KS_LCAO<TK, TR>::create_Output_Mat_Sparse(int istep)
 {
-        return ModuleIO::Output_Mat_Sparse<TK>(hsolver::HSolverLCAO<TK>::out_mat_hsR,
-            hsolver::HSolverLCAO<TK>::out_mat_dh,
-            hsolver::HSolverLCAO<TK>::out_mat_t,
-            INPUT.out_mat_r,
-            istep,
-            this->pelec->pot->get_effective_v(),
-            *this->LOWF.ParaV,
-            this->UHM,
-            this->LM,
-            this->kv,
-            this->p_hamilt);
+	return ModuleIO::Output_Mat_Sparse<TK>(hsolver::HSolverLCAO<TK>::out_mat_hsR,
+			hsolver::HSolverLCAO<TK>::out_mat_dh,
+			hsolver::HSolverLCAO<TK>::out_mat_t,
+			INPUT.out_mat_r,
+			istep,
+			this->pelec->pot->get_effective_v(),
+			*this->LOWF.ParaV,
+			this->UHM,
+			this->LM,
+			this->kv,
+			this->p_hamilt);
 }
 
 

source/module_esolver/esolver_ks_lcao.h

@@ -32,9 +32,9 @@ namespace ModuleESolver
 
         double cal_energy() override;
 
-        void cal_force(ModuleBase::matrix& force) override;
+        void cal_force(ModuleBase::matrix &force) override;
 
-        void cal_stress(ModuleBase::matrix& stress) override;
+        void cal_stress(ModuleBase::matrix &stress) override;
 
         void post_process() override;
 
@@ -89,7 +89,7 @@ namespace ModuleESolver
         ModuleBase::matrix scs;
         bool have_force = false;
 
-        void Init_Basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell);
+        void init_basis_lcao(ORB_control& orb_con, Input& inp, UnitCell& ucell);
 
         //--------------common for all calculation, not only scf-------------
         // set matrix and grid integral

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -90,7 +90,7 @@ void ESolver_KS_LCAO_TDDFT::init(Input& inp, UnitCell& ucell)
     // Init Basis should be put outside of Ensolver.
     // * reading the localized orbitals/projectors
     // * construct the interpolation tables.
-    this->Init_Basis_lcao(this->orb_con, inp, ucell);
+    this->init_basis_lcao(this->orb_con, inp, ucell);
     //------------------init Basis_lcao----------------------
 
     //------------------init Hamilt_lcao----------------------