Feature: move scf_ene_thr before mix_rho (#5375)

* add the energy threshold for SCF before mix_rho

* delete the energy threshold in iter_finish()

* calculate the energy of output density every iteration and output it

* only update_pot for scf_ene_thr>0

Author: WHUweiqingzhou

source/module_esolver/esolver_ks.cpp

@@ -547,6 +547,24 @@ void ESolver_KS<T, Device>::runner(const int istep, UnitCell& ucell)
 
             this->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(istep, iter);
+                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
+                {
+                    // update the convergence flag
+                    this->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)
@@ -669,13 +687,6 @@ void ESolver_KS<T, Device>::iter_finish(int& iter)
     }
     this->pelec->f_en.etot_delta = this->pelec->f_en.etot - this->pelec->f_en.etot_old;
     this->pelec->f_en.etot_old = this->pelec->f_en.etot;
-
-    // add a energy threshold for SCF convergence
-    if (this->scf_ene_thr > 0.0 && this->conv_esolver == 1) // only check when density is converged
-    {
-        this->conv_esolver
-            = (std::abs(this->pelec->f_en.etot_delta * ModuleBase::Ry_to_eV) < this->scf_ene_thr);
-    }
 }
 
 //! Something to do after SCF iterations when SCF is converged or comes to the max iter step.