Merge branch 'develop' of https://github.com/deepmodeling/abacus-develop into refactor

Author: YuLiu98

source/Makefile.Objects

@@ -259,6 +259,7 @@ OBJS_ESOLVER=esolver.o\
 OBJS_ESOLVER_LCAO=esolver_ks_lcao.o\
       esolver_ks_lcao_tddft.o\
       lcao_before_scf.o\
+      lcao_after_scf.o\
       esolver_gets.o\
       lcao_others.o\
 

source/module_base/tool_title.cpp

@@ -17,10 +17,10 @@ namespace ModuleBase
 
 void TITLE(const std::string &class_name,const std::string &function_name,const bool disable)
 {
-    if (disable)
-    {
-        return; // no output
-    }
+    // if (disable)
+    // {
+    //     return; // no output
+    // }
 #ifdef __NORMAL
     std::cout<<" ==> "<<class_name<<"::"<<function_name<<"\t"
 		   <<ModuleBase::GlobalFunc::MemAvailable()/1024.0/1024<<" GB\t"
@@ -37,10 +37,10 @@ void TITLE(const std::string &class_name,const std::string &function_name,const
 
 void TITLE(std::ofstream &ofs,const std::string &class_name,const std::string &function_name,const bool disable)
 {
-    if (disable)
-    {
-        return; // no output
-    }
+    // if (disable)
+    // {
+    //     return; // no output
+    // }
 #ifdef __NORMAL
     std::cout<<"\n\n ==> "<<class_name<<"::"<<function_name<<"\t"
 		   <<ModuleBase::GlobalFunc::MemAvailable()/1024.0/1024<<" GB\t"

source/module_esolver/CMakeLists.txt

@@ -17,6 +17,7 @@ if(ENABLE_LCAO)
       esolver_ks_lcao.cpp
       esolver_ks_lcao_tddft.cpp
       lcao_before_scf.cpp
+      lcao_after_scf.cpp
       esolver_gets.cpp
       lcao_others.cpp
   )

source/module_esolver/esolver_fp.cpp

@@ -130,15 +130,15 @@ void ESolver_FP::before_all_runners(UnitCell& ucell, const Input_para& inp)
 }
 
 //! 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)
+void ESolver_FP::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)
 {
     ModuleBase::TITLE("ESolver_FP", "after_scf");
 
     // 0) output convergence information
-    ModuleIO::output_convergence_after_scf(this->conv_esolver, this->pelec->f_en.etot);
+    ModuleIO::output_convergence_after_scf(conv_esolver, this->pelec->f_en.etot);
 
     // 1) write fermi energy
-    ModuleIO::output_efermi(this->conv_esolver, this->pelec->eferm.ef);
+    ModuleIO::output_efermi(conv_esolver, this->pelec->eferm.ef);
 
     // 2) update delta rho for charge extrapolation
     CE.update_delta_rho(ucell, &(this->chr), &(this->sf));
@@ -313,23 +313,12 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
         this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rhod, this->sf.strucFac);
     }
 
-    //=========================================================
-    // cal_ux should be called before init_scf because
-    // the direction of ux is used in noncoline_rho
-    //=========================================================
+    //----------------------------------------------------------
+    //! cal_ux should be called before init_scf because
+    //! the direction of ux is used in noncoline_rho
+    //----------------------------------------------------------
     elecstate::cal_ux(ucell);
 
-    this->pelec->init_scf(istep, ucell, this->Pgrid, this->sf.strucFac, this->locpp.numeric, ucell.symm);
-
-    //! Symmetry_rho should behind init_scf, because charge should be
-    //! initialized first. liuyu comment: Symmetry_rho should be located between
-    //! init_rho and v_of_rho?
-    Symmetry_rho srho;
-    for (int is = 0; is < PARAM.inp.nspin; is++)
-    {
-        srho.begin(is, *(this->pelec->charge), this->pw_rhod, ucell.symm);
-    }
-
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
     {
@@ -371,12 +360,12 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
     return;
 }
 
-void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter)
+void ESolver_FP::iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver)
 {
     //! output charge density
     if (PARAM.inp.out_chg[0] != -1)
     {
-        if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || this->conv_esolver)
+        if (iter % PARAM.inp.out_freq_elec == 0 || iter == PARAM.inp.scf_nmax || conv_esolver)
         {
             std::complex<double>** rhog_tot
                 = (PARAM.inp.dm_to_rho) ? this->pelec->charge->rhog : this->pelec->charge->rhog_save;

source/module_esolver/esolver_fp.h

@@ -39,10 +39,10 @@ class ESolver_FP : public ESolver
     virtual void before_scf(UnitCell& ucell, const int istep);
 
     //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
-    virtual void after_scf(UnitCell& ucell, const int istep);
+    virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver);
 
     //! Something to do after hamilt2density function in each iter loop.
-    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter);
+    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver);
 
     //! ------------------------------------------------------------------------------
     //! These pointers will be deleted in the free_pointers() function every ion step.

source/module_esolver/esolver_ks.cpp

@@ -437,7 +437,7 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "INIT SCF");
 
     // 4) SCF iterations
-    this->conv_esolver = false;
+    bool conv_esolver = false;
     this->niter = this->maxniter;
     this->diag_ethr = PARAM.inp.pw_diag_thr;
     for (int iter = 1; iter <= this->maxniter; ++iter)
@@ -449,10 +449,10 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
         this->hamilt2density(ucell, istep, iter, diag_ethr);
 
         // 7) finish scf iterations
-        this->iter_finish(ucell, istep, iter);
+        this->iter_finish(ucell, istep, iter, conv_esolver);
 
         // 8) check convergence
-        if (this->conv_esolver || this->oscillate_esolver)
+        if (conv_esolver || this->oscillate_esolver)
         {
             this->niter = iter;
             if (this->oscillate_esolver)
@@ -464,7 +464,7 @@ void ESolver_KS<T, Device>::runner(UnitCell& ucell, const int istep)
     } // end scf iterations
 
     // 9) after scf
-    this->after_scf(ucell, istep);
+    this->after_scf(ucell, istep, conv_esolver);
 
     ModuleBase::timer::tick(this->classname, "runner");
     return;
@@ -524,7 +524,7 @@ void ESolver_KS<T, Device>::iter_init(UnitCell& ucell, const int istep, const in
 }
 
 template <typename T, typename Device>
-void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& iter)
+void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& iter, bool &conv_esolver)
 {
     if (PARAM.inp.out_bandgap)
     {
@@ -578,30 +578,30 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
         }
 #endif
 
-        this->conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged);
+        conv_esolver = (drho < this->scf_thr && not_restart_step && is_U_converged);
 
         // add energy threshold for SCF convergence
         if (this->scf_ene_thr > 0.0)
         {
             // calculate energy of output charge density
-            this->update_pot(ucell, istep, iter);
+            this->update_pot(ucell, istep, iter, conv_esolver);
             this->pelec->cal_energies(2); // 2 means Kohn-Sham functional
             // now, etot_old is the energy of input density, while etot is the energy of output density
             this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old;
             // output etot_delta
             GlobalV::ofs_running << " DeltaE_womix = " << this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV << " eV"
                                  << std::endl;
-            if (iter > 1 && this->conv_esolver == 1) // only check when density is converged
+            if (iter > 1 && conv_esolver == 1) // only check when density is converged
             {
                 // update the convergence flag
-                this->conv_esolver
+                conv_esolver
                     = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr);
             }
         }
 
         // If drho < hsolver_error in the first iter or drho < scf_thr, we
         // do not change rho.
-        if (drho < hsolver_error || this->conv_esolver || PARAM.inp.calculation == "nscf")
+        if (drho < hsolver_error || conv_esolver || PARAM.inp.calculation == "nscf")
         {
             if (drho < hsolver_error)
             {
@@ -635,14 +635,16 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
 
 #ifdef __MPI
     MPI_Bcast(&drho, 1, MPI_DOUBLE, 0, BP_WORLD);
-    MPI_Bcast(&this->conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD);
+
+    // be careful! conv_esolver is bool, not double !! Maybe a bug 20250302 by mohan 
+    MPI_Bcast(&conv_esolver, 1, MPI_DOUBLE, 0, BP_WORLD);
     MPI_Bcast(pelec->charge->rho[0], this->pw_rhod->nrxx, MPI_DOUBLE, 0, BP_WORLD);
 #endif
 
     // update potential
     // Hamilt should be used after it is constructed.
     // this->phamilt->update(conv_esolver);
-    this->update_pot(ucell, istep, iter);
+    this->update_pot(ucell, istep, iter, conv_esolver);
 
     // 1 means Harris-Foulkes functional
     // 2 means Kohn-Sham functional
@@ -671,7 +673,7 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
         dkin = p_chgmix->get_dkin(pelec->charge, PARAM.inp.nelec);
     }
 
-    this->pelec->print_etot(ucell.magnet,this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr);
+    this->pelec->print_etot(ucell.magnet,conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr);
 
     // Json, need to be moved to somewhere else
 #ifdef __RAPIDJSON
@@ -691,17 +693,17 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
         std::cout << " SCF restart after this step!" << std::endl;
     }
 
-    ESolver_FP::iter_finish(ucell, istep, iter);
+    ESolver_FP::iter_finish(ucell, istep, iter, conv_esolver);
 }
 
 //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.
 template <typename T, typename Device>
-void ESolver_KS<T, Device>::after_scf(UnitCell& ucell, const int istep)
+void ESolver_KS<T, Device>::after_scf(UnitCell& ucell, const int istep, const bool conv_esolver)
 {
     ModuleBase::TITLE("ESolver_KS", "after_scf");
 
     // 1) call after_scf() of ESolver_FP
-    ESolver_FP::after_scf(ucell, istep);
+    ESolver_FP::after_scf(ucell, istep, conv_esolver);
 
     // 2) write eigenvalues
     if (istep % PARAM.inp.out_interval == 0)

source/module_esolver/esolver_ks.h

@@ -43,7 +43,7 @@ class ESolver_KS : public ESolver_FP
     virtual void iter_init(UnitCell& ucell, const int istep, const int iter);
 
     //! Something to do after hamilt2density function in each iter loop.
-    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter) override;
+    virtual void iter_finish(UnitCell& ucell, const int istep, int& iter, bool& conv_esolver) override;
 
     // calculate electron density from a specific Hamiltonian with ethr
     virtual void hamilt2density_single(UnitCell& ucell, const int istep, const int iter, const double ethr);
@@ -52,10 +52,10 @@ class ESolver_KS : public ESolver_FP
     void hamilt2density(UnitCell& ucell, 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(UnitCell& ucell, const int istep) override;
+    virtual void after_scf(UnitCell& ucell, const int istep, const bool conv_esolver) override;
 
     //! <Temporary> It should be replaced by a function in Hamilt Class
-    virtual void update_pot(UnitCell& ucell, const int istep, const int iter){};
+    virtual void update_pot(UnitCell& ucell, const int istep, const int iter, const bool conv_esolver){};
 
     //! Hamiltonian
     hamilt::Hamilt<T, Device>* p_hamilt = nullptr;