Refactor: Replace Esolver_OF::mu_ with efermi. (#4849)

* Refactor: Replace `Esolver_OF::mu_` with `efermi`

* [pre-commit.ci lite] apply automatic fixes

---------

Co-authored-by: pre-commit-ci-lite[bot] <117423508+pre-commit-ci-lite[bot]@users.noreply.github.com>

source/module_esolver/esolver_of.cpp

@@ -42,7 +42,6 @@ ESolver_OF::~ESolver_OF()
 
     delete[] this->nelec_;
     delete[] this->theta_;
-    delete[] this->mu_;
     delete[] this->task_;
     delete this->ptemp_rho_;
 
@@ -311,7 +310,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
 
     for (int is = 0; is < GlobalV::NSPIN; ++is)
     {
-        this->mu_[is] = 0;
+        this->pelec->eferm.get_ef(is) = 0.;
         this->theta_[is] = 0.;
         ModuleBase::GlobalFunc::ZEROS(this->pdLdphi_[is], this->pw_rho->nrxx);
         ModuleBase::GlobalFunc::ZEROS(this->pdEdphi_[is], this->pw_rho->nrxx);
@@ -348,11 +347,11 @@ void ESolver_OF::update_potential(UnitCell& ucell)
         {
             this->pdEdphi_[is][ir] = vr_eff[ir];
         }
-        this->mu_[is] = this->cal_mu(this->pphi_[is], this->pdEdphi_[is], this->nelec_[is]);
+        this->pelec->eferm.get_ef(is) = this->cal_mu(this->pphi_[is], this->pdEdphi_[is], this->nelec_[is]);
 
         for (int ir = 0; ir < this->pw_rho->nrxx; ++ir)
         {
-            this->pdLdphi_[is][ir] = this->pdEdphi_[is][ir] - 2. * this->mu_[is] * this->pphi_[is][ir];
+            this->pdLdphi_[is][ir] = this->pdEdphi_[is][ir] - 2. * this->pelec->eferm.get_efval(is) * this->pphi_[is][ir];
         }
     }
 
@@ -458,15 +457,18 @@ bool ESolver_OF::check_exit()
     bool potHold = false; // if normdLdphi nearly remains unchanged
     bool energyConv = false;
 
-    if (this->normdLdphi_ < this->of_tolp_)
+    if (this->normdLdphi_ < this->of_tolp_) {
         potConv = true;
+}
     if (this->iter_ >= 3 && std::abs(this->normdLdphi_ - this->normdLdphi_last_) < 1e-10
-        && std::abs(this->normdLdphi_ - this->normdLdphi_llast_) < 1e-10)
+        && std::abs(this->normdLdphi_ - this->normdLdphi_llast_) < 1e-10) {
         potHold = true;
+}
 
     if (this->iter_ >= 3 && std::abs(this->energy_current_ - this->energy_last_) < this->of_tole_
-        && std::abs(this->energy_current_ - this->energy_llast_) < this->of_tole_)
+        && std::abs(this->energy_current_ - this->energy_llast_) < this->of_tole_) {
         energyConv = true;
+}
 
     this->conv_ = (this->of_conv_ == "energy" && energyConv) || (this->of_conv_ == "potential" && potConv)
                   || (this->of_conv_ == "both" && potConv && energyConv);
@@ -527,7 +529,7 @@ void ESolver_OF::after_opt(const int istep, UnitCell& ucell)
                 this->pw_rhod->nx,
                 this->pw_rhod->ny,
                 this->pw_rhod->nz,
-                this->mu_[is],
+                this->pelec->eferm.get_efval(is),
                 &(ucell),
                 3,
                 1);

source/module_esolver/esolver_of.h

@@ -73,7 +73,6 @@ class ESolver_OF : public ESolver_FP
     double** pdLdphi_ = nullptr;                  // dL/dphi
     double** pphi_ = nullptr;                     // pphi[i] = ppsi.get_pointer(i), which will be freed in ~Psi().
     char* task_ = nullptr;                        // used in line search
-    double* mu_ = nullptr;                        // chemical potential
     int tn_spin_flag_ = -1;                       // spin flag used in cal_potential, which will be called by opt_tn
     int max_dcsrch_ = 200;                        // max no. of line search
     int flag_ = -1;                               // flag of TN

source/module_esolver/esolver_of_tool.cpp

@@ -83,7 +83,6 @@ void ESolver_OF::allocate_array()
     this->ptemp_rho_->set_rhopw(this->pw_rho);
     this->ptemp_rho_->allocate(GlobalV::NSPIN);
 
-    this->mu_ = new double[GlobalV::NSPIN];
     this->theta_ = new double[GlobalV::NSPIN];
     this->pdLdphi_ = new double*[GlobalV::NSPIN];
     this->pdEdphi_ = new double*[GlobalV::NSPIN];
@@ -132,9 +131,10 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi)
         }
     }
 
-    if (GlobalV::NSPIN == 4) {
+    if (GlobalV::NSPIN == 4) 
+    {
         GlobalC::ucell.cal_ux();
-}
+    }
     this->pelec->pot->update_from_charge(this->ptemp_rho_, &GlobalC::ucell);
     ModuleBase::matrix& vr_eff = this->pelec->pot->get_effective_v();
 
@@ -413,7 +413,7 @@ void ESolver_OF::print_info()
     //     maxPot = this->pdEdphi_[0][i];
     // }
     std::cout << std::setw(6) << this->iter_ << std::setw(22) << std::scientific << std::setprecision(12)
-              << this->energy_current_ / 2. << std::setw(12) << std::setprecision(3) << this->mu_[0] / 2.
+              << this->energy_current_ / 2. << std::setw(12) << std::setprecision(3) << this->pelec->eferm.get_efval(0) / 2.
               << std::setw(12) << this->theta_[0] << std::setw(12) << this->normdLdphi_ << std::setw(12)
               << (this->energy_current_ - this->energy_last_) / 2. << std::endl;
     // ============ used to compare with PROFESS3.0 ================
@@ -502,14 +502,14 @@ void ESolver_OF::print_info()
     if (GlobalV::TWO_EFERMI)
     {
         titles.push_back("E_Fermi_up");
-        energies_Ry.push_back(this->mu_[0]);
+        energies_Ry.push_back(this->pelec->eferm.get_efval(0));
         titles.push_back("E_Fermi_dw");
-        energies_Ry.push_back(this->mu_[1]);
+        energies_Ry.push_back(this->pelec->eferm.get_efval(1));
     }
     else
     {
         titles.push_back("E_Fermi");
-        energies_Ry.push_back(this->mu_[0]);
+        energies_Ry.push_back(this->pelec->eferm.get_efval(0));
     }
     energies_eV.resize(energies_Ry.size());
     std::transform(energies_Ry.begin(), energies_Ry.end(), energies_eV.begin(), [](double energy) {