Refactor:Modify Weights and Ebands Functions in the elecstate Module

source/Makefile.Objects

@@ -222,6 +222,7 @@ OBJS_ELECSTAT=elecstate.o\
     elecstate_pw_sdft.o\
     elecstate_pw_cal_tau.o\
     elecstate_op.o\
+    elecstate_tools.o\
     efield.o\
     gatefield.o\
     potential_new.o\

source/module_elecstate/CMakeLists.txt

@@ -4,6 +4,7 @@ list(APPEND objects
     elecstate_energy.cpp
     elecstate_exx.cpp
     elecstate_print.cpp
+    elecstate_tools.cpp
     elecstate_pw.cpp
     elecstate_pw_sdft.cpp
     elecstate_pw_cal_tau.cpp

source/module_elecstate/elecstate.cpp

@@ -15,43 +15,6 @@ const double* ElecState::getRho(int spin) const
     return &(this->charge->rho[spin][0]);
 }
 
-void ElecState::fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec)
-{
-    assert(nbands > 0);
-    assert(nelec > 0.0);
-
-    const double ne_thr = 1.0e-5;
-
-    const int num = this->klist->get_nks() * nbands;
-    if (num != ocp_kb.size())
-    {
-        ModuleBase::WARNING_QUIT("ElecState::fixed_weights",
-                                 "size of occupation array is wrong , please check ocp_set");
-    }
-
-    double num_elec = 0.0;
-    for (int i = 0; i < ocp_kb.size(); ++i)
-    {
-        num_elec += ocp_kb[i];
-    }
-
-    if (std::abs(num_elec - nelec) > ne_thr)
-    {
-        ModuleBase::WARNING_QUIT("ElecState::fixed_weights",
-                                 "total number of occupations is wrong , please check ocp_set");
-    }
-
-    for (int ik = 0; ik < this->wg.nr; ++ik)
-    {
-        for (int ib = 0; ib < this->wg.nc; ++ib)
-        {
-            this->wg(ik, ib) = ocp_kb[ik * this->wg.nc + ib];
-        }
-    }
-    this->skip_weights = true;
-
-    return;
-}
 
 
 void ElecState::init_nelec_spin()
@@ -64,143 +27,6 @@ void ElecState::init_nelec_spin()
     }
 }
 
-
-void ElecState::calculate_weights()
-{
-    ModuleBase::TITLE("ElecState", "calculate_weights");
-    if (this->skip_weights)
-    {
-        return;
-    }
-
-    const int nbands = this->ekb.nc;
-    const int nks = this->ekb.nr;
-
-    if (!Occupy::use_gaussian_broadening && !Occupy::fixed_occupations)
-    {
-        if (PARAM.globalv.two_fermi)
-        {
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[0],
-                             this->ekb,
-                             this->eferm.ef_up,
-                             this->wg,
-                             0,
-                             this->klist->isk);
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[1],
-                             this->ekb,
-                             this->eferm.ef_dw,
-                             this->wg,
-                             1,
-                             this->klist->isk);
-            // ef = ( ef_up + ef_dw ) / 2.0_dp need??? mohan add 2012-04-16
-        }
-        else
-        {
-            // -1 means don't need to consider spin.
-            Occupy::iweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             PARAM.inp.nelec,
-                             this->ekb,
-                             this->eferm.ef,
-                             this->wg,
-                             -1,
-                             this->klist->isk);
-        }
-    }
-    else if (Occupy::use_gaussian_broadening)
-    {
-        if (PARAM.globalv.two_fermi)
-        {
-            double demet_up = 0.0;
-            double demet_dw = 0.0;
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[0],
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef_up,
-                             demet_up,
-                             this->wg,
-                             0,
-                             this->klist->isk);
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             this->nelec_spin[1],
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef_dw,
-                             demet_dw,
-                             this->wg,
-                             1,
-                             this->klist->isk);
-            this->f_en.demet = demet_up + demet_dw;
-        }
-        else
-        {
-            // -1 means is no related to spin.
-            Occupy::gweights(nks,
-                             this->klist->wk,
-                             nbands,
-                             PARAM.inp.nelec,
-                             Occupy::gaussian_parameter,
-                             Occupy::gaussian_type,
-                             this->ekb,
-                             this->eferm.ef,
-                             this->f_en.demet,
-                             this->wg,
-                             -1,
-                             this->klist->isk);
-        }
-#ifdef __MPI
-        const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-        Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, this->f_en.demet);
-#endif
-    }
-    else if (Occupy::fixed_occupations)
-    {
-        ModuleBase::WARNING_QUIT("calculate_weights", "other occupations, not implemented");
-    }
-
-    return;
-}
-
-
-void ElecState::calEBand()
-{
-    ModuleBase::TITLE("ElecState", "calEBand");
-    // calculate ebands using wg and ekb
-    double eband = 0.0;
-#ifdef _OPENMP
-#pragma omp parallel for collapse(2) reduction(+ : eband)
-#endif
-    for (int ik = 0; ik < this->ekb.nr; ++ik)
-    {
-        for (int ibnd = 0; ibnd < this->ekb.nc; ibnd++)
-        {
-            eband += this->ekb(ik, ibnd) * this->wg(ik, ibnd);
-        }
-    }
-    this->f_en.eband = eband;
-
-#ifdef __MPI
-    const int npool = GlobalV::KPAR * PARAM.inp.bndpar;
-    Parallel_Reduce::reduce_double_allpool(npool, GlobalV::NPROC_IN_POOL, this->f_en.eband);
-#endif
-    return;
-}
-
-
 void ElecState::init_scf(const int istep, 
                          const UnitCell& ucell,
                          const Parallel_Grid& pgrid,

source/module_elecstate/elecstate.h

@@ -72,12 +72,11 @@ class ElecState
         return;
     }
 
-    // calculate wg from ekb
-    virtual void calculate_weights();
+
 
     // use occupied weights from INPUT and skip calculate_weights
     // mohan updated on 2024-06-08
-    void fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec);
+
 
     // if nupdown is not 0(TWO_EFERMI case),
     // nelec_spin will be fixed and weights will be constrained
@@ -167,11 +166,11 @@ class ElecState
     ModuleBase::matrix wg;  ///< occupation weight for each k-point and band
 
   public:
-    // calculate ebands for all k points and all occupied bands
-    void calEBand();
 
     bool skip_weights = false;
 };
+
+
 
 } // namespace elecstate
 #endif