diff --git a/source/module_cell/setup_nonlocal.cpp b/source/module_cell/setup_nonlocal.cpp
index 9fb7280c1f..f62982db49 100644
--- a/source/module_cell/setup_nonlocal.cpp
+++ b/source/module_cell/setup_nonlocal.cpp
@@ -34,7 +34,7 @@ void InfoNonlocal::Set_NonLocal(const int& it,
     ModuleBase::TITLE("InfoNonlocal", "Set_NonLocal");
 
     // set a pointer
-    // Atom* atom = &GlobalC::ucell.atoms[it];
+    // Atom* atom = &ucell.atoms[it];
 
     // get the number of non-local projectors
     n_projectors = atom->ncpp.nbeta;
diff --git a/source/module_elecstate/elecstate.cpp b/source/module_elecstate/elecstate.cpp
index c84ce504d2..d80c80cf1b 100644
--- a/source/module_elecstate/elecstate.cpp
+++ b/source/module_elecstate/elecstate.cpp
@@ -207,6 +207,7 @@ void ElecState::calEBand()
 
 
 void ElecState::init_scf(const int istep, 
+                         const UnitCell& ucell,
                          const ModuleBase::ComplexMatrix& strucfac, 
                          const bool* numeric,
                          ModuleSymmetry::Symmetry& symm, 
@@ -215,7 +216,7 @@ void ElecState::init_scf(const int istep,
     //! core correction potential.
     if (!PARAM.inp.use_paw)
     {
-        this->charge->set_rho_core(strucfac, numeric);
+        this->charge->set_rho_core(ucell,strucfac, numeric);
     }
     else
     {
@@ -226,7 +227,7 @@ void ElecState::init_scf(const int istep,
     // choose charge density from ionic step 0.
     if (istep == 0)
     {
-        this->charge->init_rho(this->eferm, strucfac, symm, (const void*)this->klist, wfcpw);
+        this->charge->init_rho(this->eferm,ucell, strucfac, symm, (const void*)this->klist, wfcpw);
         this->charge->check_rho(); // check the rho
     }
 
diff --git a/source/module_elecstate/elecstate.h b/source/module_elecstate/elecstate.h
index 7640d43d0f..ba85850b3f 100644
--- a/source/module_elecstate/elecstate.h
+++ b/source/module_elecstate/elecstate.h
@@ -104,11 +104,13 @@ class ElecState
      * @brief Init rho_core, init rho, renormalize rho, init pot
      * 
      * @param istep i-th step
+     * @param ucell unit cell
      * @param strucfac structure factor
      * @param symm symmetry
      * @param wfcpw PW basis for wave function if needed
      */
     void init_scf(const int istep,
+                  const UnitCell& ucell,
                   const ModuleBase::ComplexMatrix& strucfac,
                   const bool* numeric,
                   ModuleSymmetry::Symmetry& symm,
@@ -126,7 +128,7 @@ class ElecState
     void cal_bandgap();
     void cal_bandgap_updw();
 
-    double cal_delta_eband() const;
+    double cal_delta_eband(const UnitCell& ucell) const;
     double cal_delta_escf() const;
 
     ModuleBase::matrix vnew;
@@ -171,7 +173,8 @@ class ElecState
     ModuleBase::matrix wg;  ///< occupation weight for each k-point and band
 
   public: // print something. See elecstate_print.cpp
-    void print_etot(const bool converged,
+    void print_etot(const Magnetism& magnet,
+                    const bool converged,
                     const int& iter,
                     const double& scf_thr,
                     const double& scf_thr_kin,
diff --git a/source/module_elecstate/elecstate_energy.cpp b/source/module_elecstate/elecstate_energy.cpp
index f016da85a5..fae932bd27 100644
--- a/source/module_elecstate/elecstate_energy.cpp
+++ b/source/module_elecstate/elecstate_energy.cpp
@@ -90,7 +90,7 @@ void ElecState::cal_bandgap_updw()
 }
 
 /// @brief calculate deband
-double ElecState::cal_delta_eband() const
+double ElecState::cal_delta_eband(const UnitCell& ucell) const
 {
     // out potentials from potential mixing
     // total energy and band energy corrections
@@ -109,7 +109,7 @@ double ElecState::cal_delta_eband() const
     {
         ModuleBase::matrix v_xc;
         const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
-            = XC_Functional::v_xc(this->charge->nrxx, this->charge, &GlobalC::ucell);
+            = XC_Functional::v_xc(this->charge->nrxx, this->charge, &ucell);
         v_xc = std::get<2>(etxc_vtxc_v);
 
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
diff --git a/source/module_elecstate/elecstate_getters.cpp b/source/module_elecstate/elecstate_getters.cpp
index 30bf81f94c..bf39413498 100644
--- a/source/module_elecstate/elecstate_getters.cpp
+++ b/source/module_elecstate/elecstate_getters.cpp
@@ -8,54 +8,10 @@
 namespace elecstate
 {
 
-double get_ucell_omega()
-{
-    return GlobalC::ucell.omega;
-}
-
-double get_ucell_tpiba()
-{
-    return GlobalC::ucell.tpiba;
-}
-
 int get_xc_func_type()
 {
     return XC_Functional::get_func_type();
 }
 
-std::string get_input_vdw_method()
-{
-    return PARAM.inp.vdw_method;
-}
-
-double get_ucell_tot_magnetization()
-{
-    return GlobalC::ucell.magnet.tot_magnetization;
-}
-
-double get_ucell_abs_magnetization()
-{
-    return GlobalC::ucell.magnet.abs_magnetization;
-}
-
-double get_ucell_tot_magnetization_nc_x()
-{
-    return GlobalC::ucell.magnet.tot_magnetization_nc[0];
-}
-
-double get_ucell_tot_magnetization_nc_y()
-{
-    return GlobalC::ucell.magnet.tot_magnetization_nc[1];
-}
-
-double get_ucell_tot_magnetization_nc_z()
-{
-    return GlobalC::ucell.magnet.tot_magnetization_nc[2];
-}
-
-std::string get_ks_solver_type()
-{
-    return PARAM.inp.ks_solver;
-}
 
 } // namespace elecstate
diff --git a/source/module_elecstate/elecstate_getters.h b/source/module_elecstate/elecstate_getters.h
index 4ee4ca57c0..e9a7cdf6c1 100644
--- a/source/module_elecstate/elecstate_getters.h
+++ b/source/module_elecstate/elecstate_getters.h
@@ -7,26 +7,9 @@
 namespace elecstate
 {
 
-/// @brief get the value of GlobalC::ucell.omega
-double get_ucell_omega();
-/// @brief get the value of GlobalC::ucell.tpiba
-double get_ucell_tpiba();
 /// @brief get the value of XC_Functional::func_type
 int get_xc_func_type();
-/// @brief get the value of INPUT.vdw_method
-std::string get_input_vdw_method();
-/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization
-double get_ucell_tot_magnetization();
-/// @brief get the value of GlobalC::ucell.magnet.abs_magnetization
-double get_ucell_abs_magnetization();
-/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[0]
-double get_ucell_tot_magnetization_nc_x();
-/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[1]
-double get_ucell_tot_magnetization_nc_y();
-/// @brief get the value of GlobalC::ucell.magnet.tot_magnetization_nc[2]
-double get_ucell_tot_magnetization_nc_z();
-/// @brief get the type of KS_SOLVER
-std::string get_ks_solver_type();
+
 
 } // namespace elecstate
 
diff --git a/source/module_elecstate/elecstate_print.cpp b/source/module_elecstate/elecstate_print.cpp
index 2d738345b3..fdca966ea6 100644
--- a/source/module_elecstate/elecstate_print.cpp
+++ b/source/module_elecstate/elecstate_print.cpp
@@ -282,6 +282,7 @@ void ElecState::print_band(const int& ik, const int& printe, const int& iter)
 }
 
 /// @brief print total free energy and other energies
+/// @param ucell: unit cell
 /// @param converged: if converged
 /// @param iter_in: iter
 /// @param scf_thr: threshold for scf
@@ -289,7 +290,8 @@ void ElecState::print_band(const int& ik, const int& printe, const int& iter)
 /// @param pw_diag_thr: threshold for diagonalization
 /// @param avg_iter: averaged diagonalization iteration of each scf iteration
 /// @param print: if print to screen
-void ElecState::print_etot(const bool converged,
+void ElecState::print_etot(const Magnetism& magnet,
+                           const bool converged,
                            const int& iter_in,
                            const double& scf_thr,
                            const double& scf_thr_kin,
@@ -339,7 +341,7 @@ void ElecState::print_etot(const bool converged,
         energies_Ry.push_back(this->f_en.demet);
         titles.push_back("E_descf");
         energies_Ry.push_back(this->f_en.descf);
-        std::string vdw_method = get_input_vdw_method();
+        std::string vdw_method = PARAM.inp.vdw_method;
         if (vdw_method == "d2") // Peize Lin add 2014-04, update 2021-03-09
         {
             titles.push_back("E_vdwD2");
@@ -429,13 +431,13 @@ void ElecState::print_etot(const bool converged,
         switch (PARAM.inp.nspin)
         {
         case 2:
-            mag = {get_ucell_tot_magnetization(), get_ucell_abs_magnetization()};
+            mag = {magnet.tot_magnetization, magnet.abs_magnetization};
             break;
         case 4:
-            mag = {get_ucell_tot_magnetization_nc_x(),
-                   get_ucell_tot_magnetization_nc_y(),
-                   get_ucell_tot_magnetization_nc_z(),
-                   get_ucell_abs_magnetization()};
+            mag = {magnet.tot_magnetization_nc[0],
+                   magnet.tot_magnetization_nc[1],
+                   magnet.tot_magnetization_nc[2],
+                   magnet.abs_magnetization};
             break;
         default:
             mag = {};
@@ -446,7 +448,7 @@ void ElecState::print_etot(const bool converged,
         {
             drho.push_back(scf_thr_kin);
         }
-        elecstate::print_scf_iterinfo(get_ks_solver_type(),
+        elecstate::print_scf_iterinfo(PARAM.inp.ks_solver,
                                       iter,
                                       6,
                                       mag,
diff --git a/source/module_elecstate/elecstate_pw.cpp b/source/module_elecstate/elecstate_pw.cpp
index 696a3a75a2..05fb2b57f3 100644
--- a/source/module_elecstate/elecstate_pw.cpp
+++ b/source/module_elecstate/elecstate_pw.cpp
@@ -176,7 +176,7 @@ void ElecStatePW<T, Device>::rhoBandK(const psi::Psi<T, Device>& psi)
 
             this->basis->recip_to_real(this->ctx, &psi(ibnd,npwx), this->wfcr_another_spin, ik);
 
-            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / get_ucell_omega());
+            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / ucell->omega);
 
             if (w1 != 0.0)
             {
@@ -202,7 +202,7 @@ void ElecStatePW<T, Device>::rhoBandK(const psi::Psi<T, Device>& psi)
 
             this->basis->recip_to_real(this->ctx, &psi(ibnd,0), this->wfcr, ik);
 
-            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / get_ucell_omega());
+            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / ucell->omega);
 
             if (w1 != 0.0)
             {
@@ -222,7 +222,7 @@ void ElecStatePW<T, Device>::rhoBandK(const psi::Psi<T, Device>& psi)
                               j,
                               npw,
                               this->basis->npwk_max,
-                              static_cast<Real>(get_ucell_tpiba()),
+                              static_cast<Real>(ucell->tpiba),
                               this->basis->template get_gcar_data<Real>(),
                               this->basis->template get_kvec_c_data<Real>(),
                               &psi(ibnd, 0),
diff --git a/source/module_elecstate/elecstate_pw_cal_tau.cpp b/source/module_elecstate/elecstate_pw_cal_tau.cpp
index 0662bb4425..fd07f834af 100644
--- a/source/module_elecstate/elecstate_pw_cal_tau.cpp
+++ b/source/module_elecstate/elecstate_pw_cal_tau.cpp
@@ -26,7 +26,7 @@ void ElecStatePW<T, Device>::cal_tau(const psi::Psi<T, Device>& psi)
         {
             this->basis->recip_to_real(this->ctx, &psi(ibnd,0), this->wfcr, ik);
 
-            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / get_ucell_omega());
+            const auto w1 = static_cast<Real>(this->wg(ik, ibnd) / ucell->omega);
 
             // kinetic energy density
             for (int j = 0; j < 3; j++)
@@ -38,7 +38,7 @@ void ElecStatePW<T, Device>::cal_tau(const psi::Psi<T, Device>& psi)
                             j,
                             npw,
                             this->basis->npwk_max,
-                            static_cast<Real>(get_ucell_tpiba()),
+                            static_cast<Real>(ucell->tpiba),
                             this->basis->template get_gcar_data<Real>(),
                             this->basis->template get_kvec_c_data<Real>(),
                             &psi(ibnd, 0),
diff --git a/source/module_elecstate/magnetism.cpp b/source/module_elecstate/magnetism.cpp
index 89bf71eb38..918f1c988c 100644
--- a/source/module_elecstate/magnetism.cpp
+++ b/source/module_elecstate/magnetism.cpp
@@ -15,7 +15,11 @@ Magnetism::~Magnetism()
     delete[] this->start_magnetization;
 }
 
-void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const double* const * rho, double* nelec_spin)
+void Magnetism::compute_magnetization(const double& omega,
+                                      const int& nrxx, 
+                                      const int& nxyz, 
+                                      const double* const * rho, 
+                                      double* nelec_spin)
 {
     if (PARAM.inp.nspin==2)
     {
@@ -32,8 +36,8 @@ void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const do
         Parallel_Reduce::reduce_pool(this->tot_magnetization);
         Parallel_Reduce::reduce_pool(this->abs_magnetization);
 #endif
-        this->tot_magnetization *= elecstate::get_ucell_omega() / nxyz;
-        this->abs_magnetization *= elecstate::get_ucell_omega() / nxyz;
+        this->tot_magnetization *= omega / nxyz;
+        this->abs_magnetization *= omega / nxyz;
 
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"total magnetism (Bohr mag/cell)",this->tot_magnetization);
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization);
@@ -65,8 +69,8 @@ void Magnetism::compute_magnetization(const int& nrxx, const int& nxyz, const do
         Parallel_Reduce::reduce_pool(this->tot_magnetization_nc, 3);
         Parallel_Reduce::reduce_pool(this->abs_magnetization);
 #endif
-		for(int i=0;i<3;i++)this->tot_magnetization_nc[i] *= elecstate::get_ucell_omega() / nxyz;
-		this->abs_magnetization *= elecstate::get_ucell_omega() / nxyz;
+		for(int i=0;i<3;i++)this->tot_magnetization_nc[i] *= omega/ nxyz;
+		this->abs_magnetization *= omega/ nxyz;
 		GlobalV::ofs_running<<"total magnetism (Bohr mag/cell)"<<'\t'<<this->tot_magnetization_nc[0]<<'\t'<<this->tot_magnetization_nc[1]<<'\t'<<this->tot_magnetization_nc[2]<<'\n';
 		ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running,"absolute magnetism (Bohr mag/cell)",this->abs_magnetization);
 	}
diff --git a/source/module_elecstate/magnetism.h b/source/module_elecstate/magnetism.h
index fbc750b10c..bb4de6d28d 100644
--- a/source/module_elecstate/magnetism.h
+++ b/source/module_elecstate/magnetism.h
@@ -21,7 +21,12 @@ class Magnetism
     double tot_magnetization_nc[3];
     double abs_magnetization;
 
-    void compute_magnetization(const int& nrxx, const int& nxyz, const double* const * rho, double* nelec_spin = nullptr);
+    void compute_magnetization(const double& omega,
+                               const int& nrxx, 
+                               const int& nxyz, 
+                               const double* const * rho, 
+                               
+                               double* nelec_spin = nullptr);
 
     ModuleBase::Vector3<double> *m_loc_;   //magnetization for each element along c-axis
 	double *angle1_;           //angle between c-axis and real spin std::vector
diff --git a/source/module_elecstate/module_charge/charge.cpp b/source/module_elecstate/module_charge/charge.cpp
index 6a2405b579..c3e23dd6c6 100644
--- a/source/module_elecstate/module_charge/charge.cpp
+++ b/source/module_elecstate/module_charge/charge.cpp
@@ -218,7 +218,7 @@ double Charge::sum_rho() const
     }
 
     // multiply the sum of charge density by a factor
-    sum_rho *= elecstate::get_ucell_omega() / static_cast<double>(this->rhopw->nxyz);
+    sum_rho *= *this->omega_ / static_cast<double>(this->rhopw->nxyz);
 
 #ifdef __MPI
     Parallel_Reduce::reduce_pool(sum_rho);
@@ -722,7 +722,7 @@ double Charge::cal_rho2ne(const double* rho_in) const
 #ifdef __MPI
     Parallel_Reduce::reduce_pool(ne);
 #endif
-    ne = ne * elecstate::get_ucell_omega() / (double)this->rhopw->nxyz;
+    ne = ne * *this->omega_ / (double)this->rhopw->nxyz;
 
     return ne;
 }
diff --git a/source/module_elecstate/module_charge/charge.h b/source/module_elecstate/module_charge/charge.h
index bc8a1aa4ec..5f01773580 100644
--- a/source/module_elecstate/module_charge/charge.h
+++ b/source/module_elecstate/module_charge/charge.h
@@ -65,12 +65,14 @@ class Charge
      * @brief Init charge density from file or atomic pseudo-wave-functions
      *
      * @param eferm_iout [out] fermi energy to be initialized
+     * @param ucell [in] unit cell
      * @param strucFac [in] structure factor
      * @param symm [in] symmetry
      * @param klist [in] k points list if needed
      * @param wfcpw [in] PW basis for wave function if needed
      */
     void init_rho(elecstate::efermi& eferm_iout,
+                  const UnitCell& ucell,
                   const ModuleBase::ComplexMatrix& strucFac,
                   ModuleSymmetry::Symmetry& symm,
                   const void* klist = nullptr,
@@ -84,7 +86,9 @@ class Charge
                     const ModuleBase::ComplexMatrix& strucFac,
                     const UnitCell& ucell) const;
 
-    void set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric);
+    void set_rho_core(const UnitCell& ucell,
+                      const ModuleBase::ComplexMatrix& structure_factor, 
+                      const bool* numeric);
     void set_rho_core_paw();
 
     void renormalize_rho();
@@ -97,6 +101,8 @@ class Charge
     void non_linear_core_correction
     (
         const bool &numeric,
+        const double omega,
+        const double tpiba2,
         const int mesh,
         const double *r,
         const double *rab,
@@ -132,6 +138,8 @@ class Charge
 	   */
 	  void reduce_diff_pools(double* array_rho) const;
 
+    void set_omega(double* omega_in){this->omega_ = omega_in;};
+
     // mohan add 2021-02-20
     int nrxx; // number of r vectors in this processor
     int nxyz; // total nuber of r vectors
@@ -143,6 +151,8 @@ class Charge
 
     void destroy();    // free arrays  liuyu 2023-03-12
 
+    double* omega_ = nullptr; // omega for non-linear core correction
+
     bool allocate_rho;
 
     bool allocate_rho_final_scf; // LiuXh add 20180606
diff --git a/source/module_elecstate/module_charge/charge_init.cpp b/source/module_elecstate/module_charge/charge_init.cpp
index b38c5cdd5c..45769c3161 100644
--- a/source/module_elecstate/module_charge/charge_init.cpp
+++ b/source/module_elecstate/module_charge/charge_init.cpp
@@ -22,6 +22,7 @@
 #endif
 
 void Charge::init_rho(elecstate::efermi& eferm_iout,
+                      const UnitCell& ucell,
                       const ModuleBase::ComplexMatrix& strucFac,
                       ModuleSymmetry::Symmetry& symm,
                       const void* klist,
@@ -30,6 +31,9 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "init_chg", PARAM.inp.init_chg);
 
     std::cout << " START CHARGE      : " << PARAM.inp.init_chg << std::endl;
+    //here we need to set the omega for the charge density
+    set_omega(&ucell.omega);
+    
     bool read_error = false;
     if (PARAM.inp.init_chg == "file" || PARAM.inp.init_chg == "auto")
     {
@@ -58,7 +62,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                     GlobalV::ofs_running,
                     ssc.str(),
                     this->rho[is],
-                    GlobalC::ucell.nat))
+                    ucell.nat))
                 {
                     GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
                 }
@@ -108,7 +112,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                         GlobalV::ofs_running,
                         ssc.str(),
                         this->kin_r[is],
-                        GlobalC::ucell.nat))
+                        ucell.nat))
                     {
                         GlobalV::ofs_running << " Read in the kinetic energy density: " << ssc.str() << std::endl;
                     }
@@ -134,7 +138,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
     if (PARAM.inp.init_chg == "atomic" || 
         (PARAM.inp.init_chg == "auto" && read_error)) // mohan add 2007-10-17
     {
-        this->atomic_rho(PARAM.inp.nspin, GlobalC::ucell.omega, rho, strucFac, GlobalC::ucell);
+        this->atomic_rho(PARAM.inp.nspin, ucell.omega, rho, strucFac, ucell);
 
         // liuyu 2023-06-29 : move here from atomic_rho(), which will be called several times in charge extrapolation
         // wenfei 2021-7-29 : initial tau = 3/5 rho^2/3, Thomas-Fermi
@@ -171,7 +175,7 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
                     GlobalV::ofs_running,
                     ssc.str(),
                     this->rho[is],
-                    GlobalC::ucell.nat))
+                    ucell.nat))
                 {
                     GlobalV::ofs_running << " Read in the charge density: " << ssc.str() << std::endl;
                 }
@@ -199,7 +203,9 @@ void Charge::init_rho(elecstate::efermi& eferm_iout,
 //==========================================================
 // computes the core charge on the real space 3D mesh.
 //==========================================================
-void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric)
+void Charge::set_rho_core(const UnitCell& ucell,
+                          const ModuleBase::ComplexMatrix& structure_factor, 
+                          const bool* numeric)
 {
     ModuleBase::TITLE("Charge","set_rho_core");
     ModuleBase::timer::tick("Charge","set_rho_core");
@@ -216,9 +222,9 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con
     // int ig = 0;
 
     bool bl = false;
-    for (int it = 0; it<GlobalC::ucell.ntype; it++)
+    for (int it = 0; it<ucell.ntype; it++)
     {
-        if (GlobalC::ucell.atoms[it].ncpp.nlcc)
+        if (ucell.atoms[it].ncpp.nlcc)
         {
             bl = true;
             break;
@@ -238,9 +244,9 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con
 	// three dimension.
     std::complex<double> *vg = new std::complex<double>[this->rhopw->npw];	
 
-    for (int it = 0; it < GlobalC::ucell.ntype;it++)
+    for (int it = 0; it < ucell.ntype;it++)
     {
-        if (GlobalC::ucell.atoms[it].ncpp.nlcc)
+        if (ucell.atoms[it].ncpp.nlcc)
         {
 //----------------------------------------------------------
 // EXPLAIN : drhoc compute the radial fourier transform for
@@ -248,10 +254,12 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con
 //----------------------------------------------------------
             this->non_linear_core_correction(
                 numeric,
-                GlobalC::ucell.atoms[it].ncpp.msh,
-                GlobalC::ucell.atoms[it].ncpp.r.data(),
-                GlobalC::ucell.atoms[it].ncpp.rab.data(),
-                GlobalC::ucell.atoms[it].ncpp.rho_atc.data(),
+                ucell.omega,
+                ucell.tpiba2,
+                ucell.atoms[it].ncpp.msh,
+                ucell.atoms[it].ncpp.r.data(),
+                ucell.atoms[it].ncpp.rab.data(),
+                ucell.atoms[it].ncpp.rho_atc.data(),
                 rhocg);
 //----------------------------------------------------------
 // EXPLAIN : multiply by the structure factor and sum
@@ -296,8 +304,8 @@ void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, con
 
 	// mohan changed 2010-2-2, make this same as in atomic_rho.
 	// still lack something......
-    rhoneg /= this->rhopw->nxyz * GlobalC::ucell.omega;
-    rhoima /= this->rhopw->nxyz * GlobalC::ucell.omega;
+    rhoneg /= this->rhopw->nxyz * ucell.omega;
+    rhoima /= this->rhopw->nxyz * ucell.omega;
 
     // calculate core_only exch-corr energy etxcc=E_xc[rho_core] if required
     // The term was present in previous versions of the code but it shouldn't
@@ -322,6 +330,8 @@ void Charge::set_rho_core_paw()
 void Charge::non_linear_core_correction
 (
     const bool &numeric,
+    const double omega,
+    const double tpiba2,
     const int mesh,
     const double *r,
     const double *rab,
@@ -356,7 +366,7 @@ void Charge::non_linear_core_correction
 				}
 				ModuleBase::Integral::Simpson_Integral(mesh, aux, rab, rhocg1);
 				//rhocg [1] = fpi * rhocg1 / omega;
-				rhocg [0] = ModuleBase::FOUR_PI * rhocg1 / GlobalC::ucell.omega;//mohan modify 2008-01-19
+				rhocg [0] = ModuleBase::FOUR_PI * rhocg1 / omega;//mohan modify 2008-01-19
 			}
             igl0 = 1;
         }
@@ -370,14 +380,14 @@ void Charge::non_linear_core_correction
         // G <> 0 term
         for (int igl = igl_beg; igl < igl_end;igl++) 
         {
-            gx = sqrt(this->rhopw->gg_uniq[igl] * GlobalC::ucell.tpiba2);
+            gx = sqrt(this->rhopw->gg_uniq[igl] * tpiba2);
             ModuleBase::Sphbes::Spherical_Bessel(mesh, r, gx, 0, aux);
             for (int ir = 0;ir < mesh; ir++) 
             {
                 aux [ir] = r[ir] * r[ir] * rhoc [ir] * aux [ir];
             } //  enddo
             ModuleBase::Integral::Simpson_Integral(mesh, aux, rab, rhocg1);
-            rhocg [igl] = ModuleBase::FOUR_PI * rhocg1 / GlobalC::ucell.omega;
+            rhocg [igl] = ModuleBase::FOUR_PI * rhocg1 / omega;
         } //  enddo
         delete [] aux;
     }
diff --git a/source/module_elecstate/module_charge/charge_mixing.cpp b/source/module_elecstate/module_charge/charge_mixing.cpp
index 37860691b7..5bc3188190 100644
--- a/source/module_elecstate/module_charge/charge_mixing.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing.cpp
@@ -25,7 +25,9 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
                                const double& mixing_gg0_mag_in,
                                const double& mixing_gg0_min_in,
                                const double& mixing_angle_in,
-                               const bool& mixing_dmr_in)
+                               const bool& mixing_dmr_in,
+                               double& omega_in,
+                               double& tpiba_in)
 {
     // get private mixing parameters
     this->mixing_mode = mixing_mode_in;
@@ -38,7 +40,8 @@ void Charge_Mixing::set_mixing(const std::string& mixing_mode_in,
     this->mixing_gg0_min = mixing_gg0_min_in;
     this->mixing_angle = mixing_angle_in;
     this->mixing_dmr = mixing_dmr_in;
-
+    this->omega = &omega_in;
+    this->tpiba = &tpiba_in;
     // check the paramters
     if (this->mixing_beta > 1.0 || this->mixing_beta < 0.0)
     {
diff --git a/source/module_elecstate/module_charge/charge_mixing.h b/source/module_elecstate/module_charge/charge_mixing.h
index 9e85908d12..d6fdbc2586 100644
--- a/source/module_elecstate/module_charge/charge_mixing.h
+++ b/source/module_elecstate/module_charge/charge_mixing.h
@@ -33,6 +33,8 @@ class Charge_Mixing
      * @param mixing_gg0_min_in minimum kerker coefficient
      * @param mixing_angle_in mixing angle for nspin=4
      * @param mixing_dmr_in whether to mixing real space density matrix
+     * @param omega_in omega for non-linear core correction
+     * @param tpiba_in 2*pi/beta for non-linear core correction
      */
     void set_mixing(const std::string& mixing_mode_in,
                     const double& mixing_beta_in,
@@ -43,7 +45,9 @@ class Charge_Mixing
                     const double& mixing_gg0_mag_in,
                     const double& mixing_gg0_min_in,
                     const double& mixing_angle_in,
-                    const bool& mixing_dmr_in);
+                    const bool& mixing_dmr_in,
+                    double& omega_in,
+                    double& tpiba_in);
 
     void close_kerker_gg0() { mixing_gg0 = 0.0; mixing_gg0_mag = 0.0; }
     /**
@@ -129,7 +133,9 @@ class Charge_Mixing
     double mixing_gg0_min = 0.1;         ///< minimum kerker coefficient
     double mixing_angle = 0.0;           ///< mixing angle for nspin=4
     bool mixing_dmr = false;             ///< whether to mixing real space density matrix
-
+    double* omega = nullptr;                  ///< omega for non-linear core correction
+    double* tpiba = nullptr;                  ///< 2*pi/beta for non-linear core correction
+    double* tpiba2 = nullptr;                 ///< 2*pi/beta^2 for non-linear core correction
     std::vector<double> _drho_history; ///< history of drho used to determine the oscillation, size is scf_nmax
     
     bool new_e_iteration = true;
diff --git a/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp b/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp
index 3564b86f98..ad8ea3801b 100644
--- a/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing_preconditioner.cpp
@@ -45,7 +45,7 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
             amin = this->mixing_beta;
         }
 
-        gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
+        gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 512)
 #endif
@@ -110,7 +110,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
             amin = this->mixing_beta;
         }
         
-        gg0 = std::pow(fac * 0.529177 / GlobalC::ucell.tpiba, 2);
+        gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 512)
 #endif
diff --git a/source/module_elecstate/module_charge/charge_mixing_residual.cpp b/source/module_elecstate/module_charge/charge_mixing_residual.cpp
index d2683629d5..b2d1425597 100644
--- a/source/module_elecstate/module_charge/charge_mixing_residual.cpp
+++ b/source/module_elecstate/module_charge/charge_mixing_residual.cpp
@@ -60,9 +60,9 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
         Parallel_Reduce::reduce_pool(drho);
 #endif
         assert(nelec != 0);
-        assert(GlobalC::ucell.omega > 0);
+        assert(*this->omega > 0);
         assert(this->rhopw->nxyz > 0);
-        drho *= GlobalC::ucell.omega / static_cast<double>(this->rhopw->nxyz);
+        drho *= *this->omega / static_cast<double>(this->rhopw->nxyz);
         drho /= nelec;
     }
 
@@ -99,9 +99,9 @@ double Charge_Mixing::get_dkin(Charge* chr, const double nelec)
     Parallel_Reduce::reduce_pool(dkin);
 #endif
     assert(nelec != 0);
-    assert(GlobalC::ucell.omega > 0);
+    assert(*this->omega > 0);
     assert(this->rhopw->nxyz > 0);
-    dkin *= GlobalC::ucell.omega / static_cast<double>(this->rhopw->nxyz);
+    dkin *= *this->omega / static_cast<double>(this->rhopw->nxyz);
     dkin /= nelec;
 
     ModuleBase::timer::tick("Charge_Mixing", "get_dkin");
@@ -121,7 +121,7 @@ double Charge_Mixing::inner_product_recip_rho(std::complex<double>* rho1, std::c
         rhog2[is] = rho2 + is * this->rhopw->npw;
     }
 
-    static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2;
+    static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / ((*this->tpiba) * (*this->tpiba));
     static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI);
 
     double sum = 0.0;
@@ -245,7 +245,7 @@ double Charge_Mixing::inner_product_recip_rho(std::complex<double>* rho1, std::c
 #endif
     ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_rho");
 
-    sum *= GlobalC::ucell.omega * 0.5;
+    sum *= *this->omega * 0.5;
 
     delete[] rhog1;
     delete[] rhog2;
@@ -285,7 +285,7 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex<double>* rhog1, s
     ModuleBase::TITLE("Charge_Mixing", "inner_product_recip_hartree");
     ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree");
 
-    static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / GlobalC::ucell.tpiba2;
+    static const double fac = ModuleBase::e2 * ModuleBase::FOUR_PI / ((*this->tpiba) * (*this->tpiba));
     static const double fac2 = ModuleBase::e2 * ModuleBase::FOUR_PI / (ModuleBase::TWO_PI * ModuleBase::TWO_PI);
 
     double sum = 0.0;
@@ -446,7 +446,7 @@ double Charge_Mixing::inner_product_recip_hartree(std::complex<double>* rhog1, s
 
     ModuleBase::timer::tick("Charge_Mixing", "inner_product_recip_hartree");
 
-    sum *= GlobalC::ucell.omega * 0.5;
+    sum *= *this->omega * 0.5;
 
     return sum;
 }
diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.cpp b/source/module_elecstate/potentials/H_TDDFT_pw.cpp
index 7d92e3959a..a0533fda51 100644
--- a/source/module_elecstate/potentials/H_TDDFT_pw.cpp
+++ b/source/module_elecstate/potentials/H_TDDFT_pw.cpp
@@ -156,7 +156,7 @@ void H_TDDFT_pw::cal_v_space_length(std::vector<double>& vext_space, int direc)
     ModuleBase::TITLE("H_TDDFT_pw", "cal_v_space_length");
     ModuleBase::timer::tick("H_TDDFT_pw", "cal_v_space_length");
 
-    prepare(GlobalC::ucell, direc);
+    prepare(ucell_->G, direc);
 
     for (int ir = 0; ir < this->rho_basis_->nrxx; ++ir)
     {
@@ -436,25 +436,25 @@ double H_TDDFT_pw::cal_v_time_heaviside()
     return vext_time;
 }
 
-void H_TDDFT_pw::prepare(const UnitCell& cell, int& dir)
+void H_TDDFT_pw::prepare(const ModuleBase::Matrix3& G, int& dir)
 {
     if (dir == 1)
     {
-        bvec[0] = cell.G.e11;
-        bvec[1] = cell.G.e12;
-        bvec[2] = cell.G.e13;
+        bvec[0] = G.e11;
+        bvec[1] = G.e12;
+        bvec[2] = G.e13;
     }
     else if (dir == 2)
     {
-        bvec[0] = cell.G.e21;
-        bvec[1] = cell.G.e22;
-        bvec[2] = cell.G.e23;
+        bvec[0] = G.e21;
+        bvec[1] = G.e22;
+        bvec[2] = G.e23;
     }
     else if (dir == 3)
     {
-        bvec[0] = cell.G.e31;
-        bvec[1] = cell.G.e32;
-        bvec[2] = cell.G.e33;
+        bvec[0] = G.e31;
+        bvec[1] = G.e32;
+        bvec[2] = G.e33;
     }
     else
     {
diff --git a/source/module_elecstate/potentials/H_TDDFT_pw.h b/source/module_elecstate/potentials/H_TDDFT_pw.h
index 4925da63bd..8d7c42f9e5 100644
--- a/source/module_elecstate/potentials/H_TDDFT_pw.h
+++ b/source/module_elecstate/potentials/H_TDDFT_pw.h
@@ -95,7 +95,7 @@ class H_TDDFT_pw : public PotBase
     static std::vector<double> heavi_amp; // Ry/bohr
 
     //update At for velocity gauge by intergral of E(t)dt
-    static void update_At(void);
+    static void update_At();
 
   private:
     // internal time-step,
@@ -127,7 +127,7 @@ class H_TDDFT_pw : public PotBase
     //get ncut number for At integral
     static int check_ncut(int t_type);
 
-    void prepare(const UnitCell& cell, int& dir);
+    void prepare(const ModuleBase::Matrix3& G, int& dir);
 };
 
 } // namespace elecstate
diff --git a/source/module_elecstate/test/charge_mixing_test.cpp b/source/module_elecstate/test/charge_mixing_test.cpp
index b86f5dd64c..049ef74fc2 100644
--- a/source/module_elecstate/test/charge_mixing_test.cpp
+++ b/source/module_elecstate/test/charge_mixing_test.cpp
@@ -48,10 +48,6 @@ InfoNonlocal::~InfoNonlocal()
 }
 #endif
 // mock class cell
-namespace GlobalC
-{
-UnitCell ucell;
-};
 /************************************************
  *  unit test of charge_mixing.cpp
  ***********************************************/
@@ -89,6 +85,7 @@ UnitCell ucell;
 class ChargeMixingTest : public ::testing::Test
 {
   public:
+    UnitCell ucell;
     ChargeMixingTest()
     {
         // Init pw_basis
@@ -111,6 +108,8 @@ class ChargeMixingTest : public ::testing::Test
         PARAM.input.mixing_gg0_min = 0.1;
         PARAM.input.mixing_angle = -10.0;
         PARAM.input.mixing_dmr = false;
+        ucell.omega = 1.0;
+        ucell.tpiba = 1.0;
     }
     ModulePW::PW_Basis pw_basis;
     ModulePW::PW_Basis_Sup pw_dbasis;
@@ -138,7 +137,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     EXPECT_EQ(CMtest.get_mixing_mode(), "broyden");
     EXPECT_EQ(CMtest.get_mixing_beta(), 1.0);
     EXPECT_EQ(CMtest.get_mixing_ndim(), 1);
@@ -161,7 +162,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     EXPECT_EQ(CMtest.mixing_mode, "plain");
     EXPECT_EQ(CMtest.mixing_tau, true);
 
@@ -177,7 +180,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
                                 PARAM.input.mixing_gg0_mag,
                                 PARAM.input.mixing_gg0_min,
                                 PARAM.input.mixing_angle,
-                                PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), "");
+                                PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);, ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("You'd better set mixing_beta to [0.0, 1.0]!"));
 
@@ -194,7 +199,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
                                 PARAM.input.mixing_gg0_mag,
                                 PARAM.input.mixing_gg0_min,
                                 PARAM.input.mixing_angle,
-                                PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), "");
+                                PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);, ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("You'd better set mixing_beta_mag >= 0.0!"));
 
@@ -212,7 +219,9 @@ TEST_F(ChargeMixingTest, SetMixingTest)
                                 PARAM.input.mixing_gg0_mag,
                                 PARAM.input.mixing_gg0_min,
                                 PARAM.input.mixing_angle,
-                                PARAM.input.mixing_dmr);, ::testing::ExitedWithCode(1), "");
+                                PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);, ::testing::ExitedWithCode(1), "");
     output = testing::internal::GetCapturedStdout();
     EXPECT_THAT(output, testing::HasSubstr("This Mixing mode is not implemended yet,coming soon."));
 }
@@ -236,7 +245,9 @@ TEST_F(ChargeMixingTest, InitMixingTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     
     PARAM.input.scf_thr_type= 1;
     CMtest.init_mixing();
@@ -261,7 +272,9 @@ TEST_F(ChargeMixingTest, InitMixingTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     FUNC_TYPE = 3;
     CMtest.init_mixing();
     EXPECT_EQ(CMtest.tau_mdata.length, pw_basis.nrxx);
@@ -277,7 +290,9 @@ TEST_F(ChargeMixingTest, InitMixingTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     CMtest.init_mixing();
     EXPECT_EQ(CMtest.rho_mdata.length, 2 * pw_basis.nrxx);
 }
@@ -295,7 +310,9 @@ TEST_F(ChargeMixingTest, InnerDotRealTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     CMtest.set_rhopw(&pw_basis, &pw_basis);
     PARAM.input.nspin = 4;
 
@@ -321,7 +338,9 @@ TEST_F(ChargeMixingTest, InnerDotRealTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     PARAM.input.nspin = 4;
 
     // a simple sum for inner product
@@ -349,7 +368,9 @@ TEST_F(ChargeMixingTest, InnerDotRecipSimpleTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     CMtest.set_rhopw(&pw_basis, &pw_basis);
     PARAM.input.nspin = 2;
 
@@ -384,9 +405,20 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest)
     EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8);
 
     // RECIPROCAL NSPIN=1
-    GlobalC::ucell.tpiba2 = 1.0;
-    GlobalC::ucell.omega = 2.0;
-
+    ucell.tpiba2 = 1.0;
+    ucell.omega = 2.0;
+    CMtest.set_mixing(PARAM.input.mixing_mode,
+                    PARAM.input.mixing_beta,
+                    PARAM.input.mixing_ndim,
+                    PARAM.input.mixing_gg0,
+                    PARAM.input.mixing_tau,
+                    PARAM.input.mixing_beta_mag,
+                    PARAM.input.mixing_gg0_mag,
+                    PARAM.input.mixing_gg0_min,
+                    PARAM.input.mixing_angle,
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     PARAM.input.nspin = 1;
     std::vector<std::complex<double>> drhog1(pw_basis.npw);
     std::vector<std::complex<double>> drhog2(pw_basis.npw);
@@ -459,7 +491,9 @@ TEST_F(ChargeMixingTest, InnerDotRecipHartreeTest)
                     PARAM.input.mixing_gg0_mag,
                     PARAM.input.mixing_gg0_min,
                     PARAM.input.mixing_angle,
-                    PARAM.input.mixing_dmr);
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     drhog1.resize(pw_basis.npw * 2);
     drhog2.resize(pw_basis.npw * 2);
     for (int i = 0; i < pw_basis.npw * 2; ++i)
@@ -492,9 +526,20 @@ TEST_F(ChargeMixingTest, InnerDotRecipRhoTest)
     EXPECT_NEAR(inner, 0.5 * pw_basis.nrxx * (pw_basis.nrxx - 1), 1e-8);
 
     // RECIPROCAL
-    GlobalC::ucell.tpiba2 = 1.0;
-    GlobalC::ucell.omega = 2.0;
-
+    ucell.tpiba2 = 1.0;
+    ucell.omega = 2.0;
+    CMtest.set_mixing(PARAM.input.mixing_mode,
+                    PARAM.input.mixing_beta,
+                    PARAM.input.mixing_ndim,
+                    PARAM.input.mixing_gg0,
+                    PARAM.input.mixing_tau,
+                    PARAM.input.mixing_beta_mag,
+                    PARAM.input.mixing_gg0_mag,
+                    PARAM.input.mixing_gg0_min,
+                    PARAM.input.mixing_angle,
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     PARAM.input.nspin = 1;
     std::vector<std::complex<double>> drhog1(pw_basis.npw);
     std::vector<std::complex<double>> drhog2(pw_basis.npw);
@@ -548,7 +593,19 @@ TEST_F(ChargeMixingTest, KerkerScreenRecipTest)
 {
     Charge_Mixing CMtest;
     CMtest.set_rhopw(&pw_basis, &pw_basis);
-    GlobalC::ucell.tpiba = 1.0;
+    ucell.tpiba = 1.0;
+    CMtest.set_mixing(PARAM.input.mixing_mode,
+                    PARAM.input.mixing_beta,
+                    PARAM.input.mixing_ndim,
+                    PARAM.input.mixing_gg0,
+                    PARAM.input.mixing_tau,
+                    PARAM.input.mixing_beta_mag,
+                    PARAM.input.mixing_gg0_mag,
+                    PARAM.input.mixing_gg0_min,
+                    PARAM.input.mixing_angle,
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     // nspin = 1
     PARAM.input.nspin = 1;
     std::complex<double>* drhog = new std::complex<double>[PARAM.input.nspin*pw_basis.npw];
@@ -678,8 +735,19 @@ TEST_F(ChargeMixingTest, KerkerScreenRealTest)
 {
     Charge_Mixing CMtest;
     CMtest.set_rhopw(&pw_basis, &pw_basis);
-    GlobalC::ucell.tpiba = 1.0;
-
+    ucell.tpiba = 1.0;
+    CMtest.set_mixing(PARAM.input.mixing_mode,
+                    PARAM.input.mixing_beta,
+                    PARAM.input.mixing_ndim,
+                    PARAM.input.mixing_gg0,
+                    PARAM.input.mixing_tau,
+                    PARAM.input.mixing_beta_mag,
+                    PARAM.input.mixing_gg0_mag,
+                    PARAM.input.mixing_gg0_min,
+                    PARAM.input.mixing_angle,
+                    PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     // nspin = 1
     PARAM.input.nspin = 1;
     double* drhor = new double[PARAM.input.nspin*pw_basis.nrxx];
@@ -815,7 +883,9 @@ TEST_F(ChargeMixingTest, MixRhoTest)
                             PARAM.input.mixing_gg0_mag,
                             PARAM.input.mixing_gg0_min,
                             PARAM.input.mixing_angle,
-                            PARAM.input.mixing_dmr);
+                            PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     CMtest_recip.init_mixing();
     for(int i = 0 ; i < nspin * npw; ++i)
     {
@@ -854,7 +924,9 @@ TEST_F(ChargeMixingTest, MixRhoTest)
                         PARAM.input.mixing_gg0_mag,
                         PARAM.input.mixing_gg0_min,
                         PARAM.input.mixing_angle,
-                        PARAM.input.mixing_dmr);
+                        PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
     CMtest_real.init_mixing();
     for(int i = 0 ; i < nspin * nrxx; ++i)
     {
@@ -950,7 +1022,9 @@ TEST_F(ChargeMixingTest, MixDoubleGridRhoTest)
                             PARAM.input.mixing_gg0_mag,
                             PARAM.input.mixing_gg0_min,
                             PARAM.input.mixing_angle,
-                            PARAM.input.mixing_dmr);
+                            PARAM.input.mixing_dmr,
+                    ucell.omega,
+                    ucell.tpiba);
 
     CMtest_recip.init_mixing();
     for (int i = 0; i < nspin * npw; ++i)
diff --git a/source/module_elecstate/test/charge_test.cpp b/source/module_elecstate/test/charge_test.cpp
index f780285427..b6c98f85a2 100644
--- a/source/module_elecstate/test/charge_test.cpp
+++ b/source/module_elecstate/test/charge_test.cpp
@@ -39,10 +39,6 @@ int get_xc_func_type()
     return tmp_xc_func_type;
 }
 double tmp_ucell_omega = 500.0;
-double get_ucell_omega()
-{
-    return tmp_ucell_omega;
-}
 double tmp_gridecut = 80.0;
 void Set_GlobalV_Default()
 {
@@ -146,7 +142,7 @@ TEST_F(ChargeTest, SumRho)
             charge->rho[is][ir] = 0.1;
         }
     }
-    elecstate::tmp_ucell_omega = ucell->omega;
+    charge->set_omega(&ucell->omega);;
     EXPECT_NEAR(charge->sum_rho(), 0.1 * nspin * rhopw->nrxx * ucell->omega / rhopw->nxyz, 1E-10);
 }
 
@@ -165,7 +161,7 @@ TEST_F(ChargeTest, RenormalizeRho)
         }
     }
     EXPECT_EQ(PARAM.input.nelec, 8);
-    elecstate::tmp_ucell_omega = ucell->omega;
+    charge->set_omega(&ucell->omega);;
     charge->renormalize_rho();
     EXPECT_NEAR(charge->sum_rho(), 8.0, 1e-10);
 }
@@ -185,7 +181,7 @@ TEST_F(ChargeTest, CheckNe)
         }
     }
     EXPECT_EQ(PARAM.input.nelec, 8);
-    elecstate::tmp_ucell_omega = ucell->omega;
+    charge->set_omega(&ucell->omega);;
     charge->renormalize_rho();
     EXPECT_NEAR(charge->sum_rho(), 8.0, 1e-10);
     EXPECT_NEAR(charge->cal_rho2ne(charge->rho[0]), 8.0, 1e-10);
@@ -207,7 +203,7 @@ TEST_F(ChargeTest, SaveRhoBeforeSumBand)
     }
     EXPECT_EQ(PARAM.input.nelec, 8);
     elecstate::tmp_xc_func_type = 3;
-    elecstate::tmp_ucell_omega = ucell->omega;
+    charge->set_omega(&ucell->omega);;
     charge->renormalize_rho();
     charge->save_rho_before_sum_band();
     EXPECT_NEAR(charge->cal_rho2ne(charge->rho_save[0]), 8.0, 1e-10);
diff --git a/source/module_elecstate/test/elecstate_base_test.cpp b/source/module_elecstate/test/elecstate_base_test.cpp
index 050a50b74a..ae63d16508 100644
--- a/source/module_elecstate/test/elecstate_base_test.cpp
+++ b/source/module_elecstate/test/elecstate_base_test.cpp
@@ -32,7 +32,12 @@ Charge::Charge()
 Charge::~Charge()
 {
 }
-
+UnitCell::UnitCell(){}
+UnitCell::~UnitCell(){}
+Magnetism::Magnetism(){}
+Magnetism::~Magnetism(){}
+InfoNonlocal::InfoNonlocal(){}
+InfoNonlocal::~InfoNonlocal(){}
 #include "module_cell/klist.h"
 K_Vectors::K_Vectors()
 {
@@ -60,13 +65,14 @@ void ModulePW::PW_Basis::initgrids(double, ModuleBase::Matrix3, int, int, int)
 void ModulePW::PW_Basis::distribute_r()
 {
 }
-void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*)
+void Charge::set_rho_core(const UnitCell& ucell, ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
 void Charge::set_rho_core_paw()
 {
 }
 void Charge::init_rho(elecstate::efermi&,
+                      const UnitCell&,
                       ModuleBase::ComplexMatrix const&,
                       ModuleSymmetry::Symmetry& symm,
                       const void*,
@@ -139,6 +145,7 @@ class ElecStateTest : public ::testing::Test
 {
   protected:
     elecstate::MockElecState* elecstate;
+    UnitCell ucell;
     std::string output;
     void SetUp()
     {
@@ -249,7 +256,7 @@ TEST_F(ElecStateTest, InitSCF)
     ModuleBase::ComplexMatrix strucfac;
     elecstate->eferm = efermi;
     ModuleSymmetry::Symmetry symm;
-    EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, nullptr, symm));
+    EXPECT_NO_THROW(elecstate->init_scf(istep, ucell,strucfac, nullptr, symm));
     // delete elecstate->pot is done in the destructor of elecstate
     delete charge;
 }
diff --git a/source/module_elecstate/test/elecstate_magnetism_test.cpp b/source/module_elecstate/test/elecstate_magnetism_test.cpp
index b9266c7244..a527de29d5 100644
--- a/source/module_elecstate/test/elecstate_magnetism_test.cpp
+++ b/source/module_elecstate/test/elecstate_magnetism_test.cpp
@@ -31,10 +31,6 @@ Charge::~Charge()
 {
 }
 
-double elecstate::get_ucell_omega()
-{
-    return 500.0;
-}
 
 class MagnetismTest : public ::testing::Test
 {
@@ -57,11 +53,6 @@ TEST_F(MagnetismTest, Magnetism)
     EXPECT_EQ(nullptr, magnetism->start_magnetization);
 }
 
-TEST_F(MagnetismTest, GlobalInfo)
-{
-    EXPECT_EQ(500.0, elecstate::get_ucell_omega());
-}
-
 TEST_F(MagnetismTest, JudgeParallel)
 {
     double a[3] = {1.0, 0.0, 0.0};
@@ -91,7 +82,7 @@ TEST_F(MagnetismTest, ComputeMagnetizationS2)
                                     chr->rho[1][ir] = 1.01;
                   }
                   double* nelec_spin = new double[2];
-                  magnetism->compute_magnetization(chr->nrxx, chr->nxyz, chr->rho, nelec_spin);
+                  magnetism->compute_magnetization(500.0,chr->nrxx, chr->nxyz, chr->rho, nelec_spin);
                   EXPECT_DOUBLE_EQ(-0.5, magnetism->tot_magnetization);
                   EXPECT_DOUBLE_EQ(0.5, magnetism->abs_magnetization);
                   EXPECT_DOUBLE_EQ(4.75, nelec_spin[0]);
@@ -125,7 +116,7 @@ TEST_F(MagnetismTest, ComputeMagnetizationS4)
                                         chr->rho[3][ir] = 1.00;
                     }
                     double* nelec_spin = new double[4];
-                    magnetism->compute_magnetization(chr->nrxx, chr->nxyz, chr->rho, nelec_spin);
+                    magnetism->compute_magnetization(500.0,chr->nrxx, chr->nxyz, chr->rho, nelec_spin);
                     EXPECT_DOUBLE_EQ(100.0, magnetism->abs_magnetization);
                     EXPECT_DOUBLE_EQ(50.0*std::sqrt(2.0), magnetism->tot_magnetization_nc[0]);
                     EXPECT_DOUBLE_EQ(50.0, magnetism->tot_magnetization_nc[1]);
diff --git a/source/module_elecstate/test/elecstate_print_test.cpp b/source/module_elecstate/test/elecstate_print_test.cpp
index 657ac999b8..0fb0248e64 100644
--- a/source/module_elecstate/test/elecstate_print_test.cpp
+++ b/source/module_elecstate/test/elecstate_print_test.cpp
@@ -33,38 +33,14 @@ void ElecState::calculate_weights()
 } // just for mock
 double Efield::etotefield = 1.1;
 double elecstate::Gatefield::etotgatefield = 2.2;
-std::string tmp_vdw_method = "d2";
-std::string get_input_vdw_method()
-{
-    return tmp_vdw_method;
-}
-double get_ucell_tot_magnetization()
-{
-    return 1.1;
-}
-double get_ucell_abs_magnetization()
-{
-    return 2.2;
-}
-double get_ucell_tot_magnetization_nc_x()
-{
-    return 3.3;
-}
-double get_ucell_tot_magnetization_nc_y()
-{
-    return 4.4;
-}
-double get_ucell_tot_magnetization_nc_z()
-{
-    return 5.5;
-}
-std::string tmp_ks_solver = "dav";
-std::string get_ks_solver_type()
-{
-    return tmp_ks_solver;
-}
-} // namespace elecstate
 
+} // namespace elecstate
+UnitCell::UnitCell(){}
+UnitCell::~UnitCell(){}
+Magnetism::Magnetism(){}
+Magnetism::~Magnetism(){}
+InfoNonlocal::InfoNonlocal(){}
+InfoNonlocal::~InfoNonlocal(){}
 Charge::Charge()
 {
 }
@@ -90,6 +66,7 @@ class ElecStatePrintTest : public ::testing::Test
 {
   protected:
     elecstate::ElecState elecstate;
+    UnitCell ucell;
     std::string output;
     std::ifstream ifs;
     std::ofstream ofs;
@@ -117,6 +94,12 @@ class ElecStatePrintTest : public ::testing::Test
         elecstate.wg(0, 1) = 0.2;
         elecstate.wg(1, 0) = 0.3;
         elecstate.wg(1, 1) = 0.4;
+        ucell.magnet.tot_magnetization = 1.1;
+        ucell.magnet.abs_magnetization = 2.2;
+        ucell.magnet.tot_magnetization_nc[0] = 3.3;
+        ucell.magnet.tot_magnetization_nc[1] = 4.4;
+        ucell.magnet.tot_magnetization_nc[2] = 5.5;
+        PARAM.input.ks_solver = "dav";
     }
     void TearDown()
     {
@@ -251,38 +234,38 @@ TEST_F(ElecStatePrintTest, PrintEtot)
     std::vector<std::string> vdw_methods = {"d2", "d3_0", "d3_bj"};
     for (int i = 0; i < vdw_methods.size(); i++)
     {
-        elecstate::tmp_vdw_method = vdw_methods[i];
-        elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, false);
+        PARAM.input.vdw_method = vdw_methods[i];
+        elecstate.print_etot(ucell.magnet,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, false);
     }
     // iteration of different ks_solver
     std::vector<std::string> ks_solvers = {"cg", "lapack", "genelpa", "dav", "scalapack_gvx", "cusolver"};
     for (int i = 0; i < ks_solvers.size(); i++)
     {
-        elecstate::tmp_ks_solver = ks_solvers[i];
+        PARAM.input.ks_solver = ks_solvers[i];
         testing::internal::CaptureStdout();
-        elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
+        elecstate.print_etot(ucell.magnet,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
         output = testing::internal::GetCapturedStdout();
-        if (elecstate::tmp_ks_solver == "cg")
+        if (PARAM.input.ks_solver == "cg")
         {
             EXPECT_THAT(output, testing::HasSubstr("CG"));
         }
-        else if (elecstate::tmp_ks_solver == "lapack")
+        else if (PARAM.input.ks_solver == "lapack")
         {
             EXPECT_THAT(output, testing::HasSubstr("LA"));
         }
-        else if (elecstate::tmp_ks_solver == "genelpa")
+        else if (PARAM.input.ks_solver == "genelpa")
         {
             EXPECT_THAT(output, testing::HasSubstr("GE"));
         }
-        else if (elecstate::tmp_ks_solver == "dav")
+        else if (PARAM.input.ks_solver == "dav")
         {
             EXPECT_THAT(output, testing::HasSubstr("DA"));
         }
-        else if (elecstate::tmp_ks_solver == "scalapack_gvx")
+        else if (PARAM.input.ks_solver == "scalapack_gvx")
         {
             EXPECT_THAT(output, testing::HasSubstr("GV"));
         }
-        else if (elecstate::tmp_ks_solver == "cusolver")
+        else if (PARAM.input.ks_solver == "cusolver")
         {
             EXPECT_THAT(output, testing::HasSubstr("CU"));
         }
@@ -328,7 +311,7 @@ TEST_F(ElecStatePrintTest, PrintEtot2)
     PARAM.input.basis_type = "pw";
     PARAM.input.scf_nmax = 100;
 
-    elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
+    elecstate.print_etot(ucell.magnet,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
     GlobalV::ofs_running.close();
     ifs.open("test.dat", std::ios::in);
     std::string str((std::istreambuf_iterator<char>(ifs)), std::istreambuf_iterator<char>());
@@ -364,7 +347,7 @@ TEST_F(ElecStatePrintTest, PrintEtotColorS2)
     PARAM.input.out_bandgap = true;
     PARAM.input.nspin = 2;
     GlobalV::MY_RANK = 0;
-    elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
+    elecstate.print_etot(ucell.magnet,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
 }
 
 TEST_F(ElecStatePrintTest, PrintEtotColorS4)
@@ -389,5 +372,5 @@ TEST_F(ElecStatePrintTest, PrintEtotColorS4)
     PARAM.input.nspin = 4;
     PARAM.input.noncolin = true;
     GlobalV::MY_RANK = 0;
-    elecstate.print_etot(converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
+    elecstate.print_etot(ucell.magnet,converged, iter, scf_thr, scf_thr_kin, duration, printe, pw_diag_thr, avg_iter, print);
 }
\ No newline at end of file
diff --git a/source/module_elecstate/test/elecstate_pw_test.cpp b/source/module_elecstate/test/elecstate_pw_test.cpp
index 85ba4d1b24..dae7ab46a4 100644
--- a/source/module_elecstate/test/elecstate_pw_test.cpp
+++ b/source/module_elecstate/test/elecstate_pw_test.cpp
@@ -11,14 +11,6 @@
 // mock functions for testing
 namespace elecstate
 {
-double get_ucell_omega()
-{
-    return 500.0;
-}
-double get_ucell_tpiba()
-{
-    return 2.0;
-}
 int tmp_xc_func_type = 1;
 int get_xc_func_type()
 {
@@ -141,13 +133,14 @@ K_Vectors::~K_Vectors()
 {
 }
 
-void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*)
+void Charge::set_rho_core(const UnitCell& ucell, ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
 void Charge::set_rho_core_paw()
 {
 }
 void Charge::init_rho(elecstate::efermi&,
+                      const UnitCell&,
                       ModuleBase::ComplexMatrix const&,
                       ModuleSymmetry::Symmetry& symm,
                       const void*,
@@ -221,6 +214,8 @@ class ElecStatePWTest : public ::testing::Test
         klist = new K_Vectors;
         klist->set_nks(5);
         ucell = new UnitCell;
+        ucell->omega = 500.0;
+        ucell->tpiba = 2.0;
         ppcell = new pseudopot_cell_vnl;
         rhodpw = new ModulePW::PW_Basis;
         rhopw = new ModulePW::PW_Basis;
diff --git a/source/module_esolver/esolver_ks.cpp b/source/module_esolver/esolver_ks.cpp
index e644d19e0a..d5d596712b 100644
--- a/source/module_esolver/esolver_ks.cpp
+++ b/source/module_esolver/esolver_ks.cpp
@@ -108,7 +108,9 @@ void ESolver_KS<T, Device>::before_all_runners(UnitCell& ucell, const Input_para
                          PARAM.inp.mixing_gg0_mag,
                          PARAM.inp.mixing_gg0_min,
                          PARAM.inp.mixing_angle,
-                         PARAM.inp.mixing_dmr);
+                         PARAM.inp.mixing_dmr,
+                         ucell.omega,
+                         ucell.tpiba);
     p_chgmix->init_mixing();
 
     /// PAW Section
@@ -546,7 +548,8 @@ void ESolver_KS<T, Device>::iter_finish(UnitCell& ucell, const int istep, int& i
     }
 
     // compute magnetization, only for LSDA(spin==2)
-    ucell.magnet.compute_magnetization(this->pelec->charge->nrxx,
+    ucell.magnet.compute_magnetization(ucell.omega,
+                                       this->pelec->charge->nrxx,
                                        this->pelec->charge->nxyz,
                                        this->pelec->charge->rho,
                                        this->pelec->nelec_spin.data());
@@ -671,7 +674,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(this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr);
+    this->pelec->print_etot(ucell.magnet,this->conv_esolver, iter, drho, dkin, duration, PARAM.inp.printe, diag_ethr);
 
     // Json, need to be moved to somewhere else
 #ifdef __RAPIDJSON
diff --git a/source/module_esolver/esolver_ks_lcao.cpp b/source/module_esolver/esolver_ks_lcao.cpp
index 44e8105dc9..65c4c774b1 100644
--- a/source/module_esolver/esolver_ks_lcao.cpp
+++ b/source/module_esolver/esolver_ks_lcao.cpp
@@ -570,7 +570,7 @@ void ESolver_KS_LCAO<TK, TR>::iter_init(UnitCell& ucell, const int istep, const
     }
 
     // mohan update 2012-06-05
-    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();
+    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
 
     // mohan move it outside 2011-01-13
     // first need to calculate the weight according to
@@ -760,7 +760,7 @@ void ESolver_KS_LCAO<TK, TR>::hamilt2density_single(UnitCell& ucell, int istep,
     }
 
     // 12) calculate delta energy
-    this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+    this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
 }
 
 //------------------------------------------------------------------------------
diff --git a/source/module_esolver/esolver_ks_lcao_tddft.cpp b/source/module_esolver/esolver_ks_lcao_tddft.cpp
index fc1ba66bc6..6b50347716 100644
--- a/source/module_esolver/esolver_ks_lcao_tddft.cpp
+++ b/source/module_esolver/esolver_ks_lcao_tddft.cpp
@@ -139,7 +139,7 @@ void ESolver_KS_LCAO_TDDFT::hamilt2density_single(UnitCell& ucell, const int ist
     }
 
     // (7) calculate delta energy
-    this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+    this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
 }
 
 void ESolver_KS_LCAO_TDDFT::iter_finish(UnitCell& ucell, const int istep, int& iter)
diff --git a/source/module_esolver/esolver_ks_lcaopw.cpp b/source/module_esolver/esolver_ks_lcaopw.cpp
index 1060a80d29..e79346cd6d 100644
--- a/source/module_esolver/esolver_ks_lcaopw.cpp
+++ b/source/module_esolver/esolver_ks_lcaopw.cpp
@@ -166,7 +166,7 @@ namespace ModuleESolver
         // deband is calculated from "output" charge density calculated
         // in sum_band
         // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))'
-        this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+        this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
 
         ModuleBase::timer::tick("ESolver_KS_LIP", "hamilt2density_single");
     }
diff --git a/source/module_esolver/esolver_ks_pw.cpp b/source/module_esolver/esolver_ks_pw.cpp
index 385cf27621..5248bcce0b 100644
--- a/source/module_esolver/esolver_ks_pw.cpp
+++ b/source/module_esolver/esolver_ks_pw.cpp
@@ -307,7 +307,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     //! calculate the total local pseudopotential in real space
-    this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc);
+    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc);
 
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
@@ -481,7 +481,7 @@ void ESolver_KS_PW<T, Device>::iter_init(UnitCell& ucell, const int istep, const
     }
     // mohan move harris functional to here, 2012-06-05
     // use 'rho(in)' and 'v_h and v_xc'(in)
-    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband();
+    this->pelec->f_en.deband_harris = this->pelec->cal_delta_eband(ucell);
 
     // update local occupations for DFT+U
     // should before lambda loop in DeltaSpin
@@ -576,7 +576,7 @@ void ESolver_KS_PW<T, Device>::hamilt2density_single(UnitCell& ucell,
     // deband is calculated from "output" charge density calculated
     // in sum_band
     // need 'rho(out)' and 'vr (v_h(in) and v_xc(in))'
-    this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+    this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
 
     ModuleBase::timer::tick("ESolver_KS_PW", "hamilt2density_single");
 }
diff --git a/source/module_esolver/esolver_of.cpp b/source/module_esolver/esolver_of.cpp
index bdb24829e4..7abb63553b 100644
--- a/source/module_esolver/esolver_of.cpp
+++ b/source/module_esolver/esolver_of.cpp
@@ -120,7 +120,7 @@ void ESolver_OF::before_all_runners(UnitCell& ucell, const Input_para& inp)
     this->init_elecstate(ucell);
 
     // calculate the total local pseudopotential in real space
-    this->pelec->init_scf(0, sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs
+    this->pelec->init_scf(0, ucell,sf.strucFac, locpp.numeric, ucell.symm); // atomic_rho, v_of_rho, set_vrs
 
     // liuyu move here 2023-10-09
     // D in uspp need vloc, thus behind init_scf()
@@ -258,7 +258,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
             GlobalV::ofs_warning);
     }
 
-    this->pelec->init_scf(istep, sf.strucFac, locpp.numeric, ucell.symm);
+    this->pelec->init_scf(istep,ucell, sf.strucFac, locpp.numeric, ucell.symm);
 
     // calculate ewald energy
     this->pelec->f_en.ewald_energy = H_Ewald_pw::compute_ewald(ucell, this->pw_rho, sf.strucFac);
diff --git a/source/module_esolver/esolver_sdft_pw.cpp b/source/module_esolver/esolver_sdft_pw.cpp
index abeda884b1..ab9141fcb5 100644
--- a/source/module_esolver/esolver_sdft_pw.cpp
+++ b/source/module_esolver/esolver_sdft_pw.cpp
@@ -198,7 +198,7 @@ void ESolver_SDFT_PW<T, Device>::hamilt2density_single(UnitCell& ucell, int iste
         {
             srho.begin(is, *(this->pelec->charge), this->pw_rho, ucell.symm);
         }
-        this->pelec->f_en.deband = this->pelec->cal_delta_eband();
+        this->pelec->f_en.deband = this->pelec->cal_delta_eband(ucell);
     }
     else
     {
diff --git a/source/module_esolver/lcao_before_scf.cpp b/source/module_esolver/lcao_before_scf.cpp
index 9eef0aafce..6ce0994d44 100644
--- a/source/module_esolver/lcao_before_scf.cpp
+++ b/source/module_esolver/lcao_before_scf.cpp
@@ -260,7 +260,7 @@ void ESolver_KS_LCAO<TK, TR>::before_scf(UnitCell& ucell, const int istep)
     }
 #endif // __EXX
 
-    this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm);
+    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm);
 
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
diff --git a/source/module_esolver/lcao_others.cpp b/source/module_esolver/lcao_others.cpp
index 55ee33c593..a14e9e203d 100644
--- a/source/module_esolver/lcao_others.cpp
+++ b/source/module_esolver/lcao_others.cpp
@@ -249,7 +249,7 @@ void ESolver_KS_LCAO<TK, TR>::others(UnitCell& ucell, const int istep)
     elecstate::cal_ux(ucell);
 
     // pelec should be initialized before these calculations
-    this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm);
+    this->pelec->init_scf(istep, ucell,this->sf.strucFac, this->ppcell.numeric, ucell.symm);
     // self consistent calculations for electronic ground state
     if (cal_type == "get_pchg")
     {
diff --git a/source/module_hsolver/test/hsolver_supplementary_mock.h b/source/module_hsolver/test/hsolver_supplementary_mock.h
index 2ee381ad48..3ec0817825 100644
--- a/source/module_hsolver/test/hsolver_supplementary_mock.h
+++ b/source/module_hsolver/test/hsolver_supplementary_mock.h
@@ -43,6 +43,7 @@ void ElecState::print_eigenvalue(std::ofstream& ofs)
 }
 
 void ElecState::init_scf(const int istep,
+                         const UnitCell& ucell,
                          const ModuleBase::ComplexMatrix& strucfac,
                          const bool*,
                          ModuleSymmetry::Symmetry&,
diff --git a/source/module_rdmft/rdmft.cpp b/source/module_rdmft/rdmft.cpp
index a37be99626..3649c0fa26 100644
--- a/source/module_rdmft/rdmft.cpp
+++ b/source/module_rdmft/rdmft.cpp
@@ -343,7 +343,7 @@ void RDMFT<TK, TR>::cal_Energy(const int cal_type)
     }
     else
     {
-        this->pelec->f_en.deband  = this->pelec->cal_delta_eband();
+        this->pelec->f_en.deband  = this->pelec->cal_delta_eband(*ucell);
         E_descf = pelec->f_en.descf = 0.0;
         this->pelec->cal_energies(2);
         Etotal = this->pelec->f_en.etot;