Refactor: runner() of esolver_ks (#5445)

* Refactor: runner() of esolver_ks

* rename hamilt2density and diag as hamilt2density_single and hamilt2density

Author: YuLiu98

source/module_esolver/esolver_ks.cpp

@@ -46,12 +46,15 @@ class ESolver_KS : public ESolver_FP
     //! Something to do after hamilt2density function in each iter loop.
     virtual void iter_finish(const int istep, int& iter);
 
-    // calculate electron density from a specific Hamiltonian
-    virtual void hamilt2density(const int istep, const int iter, const double ethr);
+    // calculate electron density from a specific Hamiltonian with ethr
+    virtual void hamilt2density_single(const int istep, const int iter, const double ethr);
 
     // calculate electron states from a specific Hamiltonian
     virtual void hamilt2estates(const double ethr) {};
 
+    // calculate electron density from a specific Hamiltonian
+    void hamilt2density(const int istep, const int iter, const double ethr);
+
     //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
     virtual void after_scf(const int istep) override;
 
@@ -82,6 +85,7 @@ class ESolver_KS : public ESolver_FP
     double scf_thr;                 // scf density threshold
     double scf_ene_thr;             // scf energy threshold
     double drho;                    // the difference between rho_in (before HSolver) and rho_out (After HSolver)
+    double hsolver_error;           // the error of HSolver
     int maxniter;                   // maximum iter steps for scf
     int niter;                      // iter steps actually used in scf
     int out_freq_elec;              // frequency for output

source/module_esolver/esolver_ks.h

@@ -681,10 +681,17 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
         SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
         sc.run_lambda_loop(iter - 1);
     }
+
+    // save density matrix DMR for mixing
+    if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
+    {
+        elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
+        dm->save_DMR();
+    }
 }
 
 //------------------------------------------------------------------------------
-//! the 11th function of ESolver_KS_LCAO: hamilt2density
+//! 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
@@ -700,47 +707,16 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(const int istep, const int iter)
 //! 12) calculate delta energy
 //------------------------------------------------------------------------------
 template <typename TK, typename TR>
-void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
+void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(int istep, int iter, double ethr)
 {
-    ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density");
-
-    // 1) save input rho
-    this->pelec->charge->save_rho_before_sum_band();
-
-    // 2) save density matrix DMR for mixing
-    if (PARAM.inp.mixing_restart > 0 && PARAM.inp.mixing_dmr && this->p_chgmix->mixing_restart_count > 0)
-    {
-        elecstate::DensityMatrix<TK, double>* dm = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM();
-        dm->save_DMR();
-    }
+    ModuleBase::TITLE("ESolver_KS_LCAO", "hamilt2density_single");
 
-    // 3) solve the Hamiltonian and output band gap
-    {
-        // reset energy
-        this->pelec->f_en.eband = 0.0;
-        this->pelec->f_en.demet = 0.0;
-
-        hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
-        hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false);
+    // reset energy
+    this->pelec->f_en.eband = 0.0;
+    this->pelec->f_en.demet = 0.0;
 
-        if (PARAM.inp.out_bandgap)
-        {
-            if (!PARAM.globalv.two_fermi)
-            {
-                this->pelec->cal_bandgap();
-            }
-            else
-            {
-                this->pelec->cal_bandgap_updw();
-            }
-        }
-    }
-
-    // 4) print bands for each k-point and each band
-    for (int ik = 0; ik < this->kv.get_nks(); ++ik)
-    {
-        this->pelec->print_band(ik, PARAM.inp.printe, iter);
-    }
+    hsolver::HSolverLCAO<TK> hsolver_lcao_obj(&(this->pv), PARAM.inp.ks_solver);
+    hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec, false);
 
     // 5) what's the exd used for?
 #ifdef __EXX
@@ -754,59 +730,13 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density(int istep, int iter, double ethr)
     }
 #endif
 
-    // 6) 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,
-        // new DFT+U method will calculate energy in calculating Hamiltonian
-        if (PARAM.inp.dft_plus_u == 2)
-        {
-            if (GlobalC::dftu.omc != 2)
-            {
-                const std::vector<std::vector<TK>>& tmp_dm
-                    = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
-                this->dftu_cal_occup_m(iter, tmp_dm);
-            }
-            GlobalC::dftu.cal_energy_correction(istep);
-        }
-        GlobalC::dftu.output();
-    }
-
-    // (7) for deepks, calculate delta_e
-#ifdef __DEEPKS
-    if (PARAM.inp.deepks_scf)
-    {
-        const std::vector<std::vector<TK>>& dm
-            = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
-
-        this->dpks_cal_e_delta_band(dm);
-    }
-#endif
-
-    // 8) for delta spin
-    if (PARAM.inp.sc_mag_switch)
-    {
-        SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
-        sc.cal_MW(iter, this->p_hamilt);
-    }
-
-    // 9) use new charge density to calculate energy
-    this->pelec->cal_energies(1);
-
     // 10) symmetrize the charge density
     Symmetry_rho srho;
     for (int is = 0; is < PARAM.inp.nspin; is++)
     {
         srho.begin(is, *(this->pelec->charge), this->pw_rho, GlobalC::ucell.symm);
     }
 
-    // 11) compute magnetization, only for spin==2
-    GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
-                                                this->pelec->charge->nxyz,
-                                                this->pelec->charge->rho,
-                                                this->pelec->nelec_spin.data());
-
     // 12) calculate delta energy
     this->pelec->f_en.deband = this->pelec->cal_delta_eband();
 }
@@ -923,6 +853,43 @@ void ESolver_KS_LCAO<TK, TR>::iter_finish(const int istep, int& iter)
 {
     ModuleBase::TITLE("ESolver_KS_LCAO", "iter_finish");
 
+    // 6) 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,
+        // new DFT+U method will calculate energy in calculating Hamiltonian
+        if (PARAM.inp.dft_plus_u == 2)
+        {
+            if (GlobalC::dftu.omc != 2)
+            {
+                const std::vector<std::vector<TK>>& tmp_dm
+                    = dynamic_cast<elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
+                this->dftu_cal_occup_m(iter, tmp_dm);
+            }
+            GlobalC::dftu.cal_energy_correction(istep);
+        }
+        GlobalC::dftu.output();
+    }
+
+    // (7) for deepks, calculate delta_e
+#ifdef __DEEPKS
+    if (PARAM.inp.deepks_scf)
+    {
+        const std::vector<std::vector<TK>>& dm
+            = dynamic_cast<const elecstate::ElecStateLCAO<TK>*>(this->pelec)->get_DM()->get_DMK_vector();
+
+        this->dpks_cal_e_delta_band(dm);
+    }
+#endif
+
+    // 8) for delta spin
+    if (PARAM.inp.sc_mag_switch)
+    {
+        SpinConstrain<TK, base_device::DEVICE_CPU>& sc = SpinConstrain<TK, base_device::DEVICE_CPU>::getScInstance();
+        sc.cal_MW(iter, this->p_hamilt);
+    }
+
     // call iter_finish() of ESolver_KS
     ESolver_KS<TK>::iter_finish(istep, iter);
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -50,9 +50,7 @@ class ESolver_KS_LCAO : public ESolver_KS<TK> {
 
     virtual void iter_init(const int istep, const int iter) override;
 
-    virtual void hamilt2density(const int istep,
-                                const int iter,
-                                const double ethr) override;
+    virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override;
 
     virtual void update_pot(const int istep, const int iter) override;
 

source/module_esolver/esolver_ks_lcao.h

@@ -118,10 +118,8 @@ void ESolver_KS_LCAO_TDDFT::before_all_runners(const Input_para& inp, UnitCell&
     this->pelec_td = dynamic_cast<elecstate::ElecStateLCAO_TDDFT*>(this->pelec);
 }
 
-void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, const double ethr)
+void ESolver_KS_LCAO_TDDFT::hamilt2density_single(const int istep, const int iter, const double ethr)
 {
-    pelec->charge->save_rho_before_sum_band();
-
     if (wf.init_wfc == "file")
     {
         if (istep >= 1)
@@ -171,11 +169,23 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons
             hsolver_lcao_obj.solve(this->p_hamilt, this->psi[0], this->pelec_td, false);
         }
     }
-    // else
-    // {
-    //     ModuleBase::WARNING_QUIT("ESolver_KS_LCAO", "HSolver has not been initialed!");
-    // }
 
+    // symmetrize the charge density only for ground state
+    if (istep <= 1)
+    {
+        Symmetry_rho srho;
+        for (int is = 0; is < PARAM.inp.nspin; is++)
+        {
+            srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm);
+        }
+    }
+
+    // (7) calculate delta energy
+    this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+}
+
+void ESolver_KS_LCAO_TDDFT::iter_finish(const int istep, int& iter)
+{
     // print occupation of each band
     if (iter == 1 && istep <= 2)
     {
@@ -201,32 +211,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density(const int istep, const int iter, cons
                              << std::endl;
     }
 
-    for (int ik = 0; ik < kv.get_nks(); ++ik)
-    {
-        this->pelec_td->print_band(ik, PARAM.inp.printe, iter);
-    }
-
-    // using new charge density.
-    this->pelec->cal_energies(1);
-
-    // symmetrize the charge density only for ground state
-    if (istep <= 1)
-    {
-        Symmetry_rho srho;
-        for (int is = 0; is < PARAM.inp.nspin; is++)
-        {
-            srho.begin(is, *(pelec->charge), pw_rho, GlobalC::ucell.symm);
-        }
-    }
-
-    // (6) compute magnetization, only for spin==2
-    GlobalC::ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
-                                                this->pelec->charge->nxyz,
-                                                this->pelec->charge->rho,
-                                                pelec->nelec_spin.data());
-
-    // (7) calculate delta energy
-    this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+    ESolver_KS_LCAO<std::complex<double>, double>::iter_finish(istep, iter);
 }
 
 void ESolver_KS_LCAO_TDDFT::update_pot(const int istep, const int iter)

source/module_esolver/esolver_ks_lcao_tddft.cpp

@@ -30,10 +30,12 @@ class ESolver_KS_LCAO_TDDFT : public ESolver_KS_LCAO<std::complex<double>, doubl
     int td_htype = 1;
 
   protected:
-    virtual void hamilt2density(const int istep, const int iter, const double ethr) override;
+    virtual void hamilt2density_single(const int istep, const int iter, const double ethr) override;
 
     virtual void update_pot(const int istep, const int iter) override;
 
+    virtual void iter_finish(const int istep, int& iter) override;
+
     virtual void after_scf(const int istep) override;
 
     void cal_edm_tddft();