Refactor::Remove GlobalC::ucell in module_elecstate

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;

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(ucell,this->eferm, strucfac, symm, (const void*)this->klist, wfcpw);
         this->charge->check_rho(); // check the rho
     }
 

source/module_elecstate/elecstate.h

@@ -109,6 +109,7 @@ class ElecState
      * @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,
@@ -124,9 +125,9 @@ class ElecState
 
   public: // something aboud energies. See elecstate_energy.cpp
     void cal_bandgap();
-    void cal_bandgap_updw();
+    void cal_bandgap_updw(const UnitCell& ucell);
 
-    double cal_delta_eband() const;
+    double cal_delta_eband(const UnitCell& ucell) const;
     double cal_delta_escf() const;
 
     ModuleBase::matrix vnew;
@@ -173,7 +174,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 UnitCell& ucell,
+                    const bool converged,
                     const int& iter,
                     const double& scf_thr,
                     const double& scf_thr_kin,

source/module_elecstate/elecstate_energy.cpp

@@ -43,7 +43,7 @@ void ElecState::cal_bandgap()
 
 /// @brief calculate spin up & down band gap
 /// @todo add isk[ik] so as to discriminate different spins
-void ElecState::cal_bandgap_updw()
+void ElecState::cal_bandgap_updw(const UnitCell& ucell)
 {
     if (this->ekb.nr == 0 || this->ekb.nc == 0)
     { // which means no homo and no lumo
@@ -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++)

source/module_elecstate/elecstate_getters.cpp

@@ -8,54 +8,13 @@
 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

source/module_elecstate/elecstate_getters.h

@@ -7,26 +7,8 @@
 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
 

source/module_elecstate/elecstate_print.cpp

@@ -289,7 +289,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 UnitCell& ucell,
+                           const bool converged,
                            const int& iter_in,
                            const double& scf_thr,
                            const double& scf_thr_kin,
@@ -339,7 +340,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 +430,13 @@ void ElecState::print_etot(const bool converged,
         switch (PARAM.inp.nspin)
         {
         case 2:
-            mag = {get_ucell_tot_magnetization(), get_ucell_abs_magnetization()};
+            mag = {ucell.magnet.tot_magnetization, ucell.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 = {ucell.magnet.tot_magnetization_nc[0],
+                   ucell.magnet.tot_magnetization_nc[1],
+                   ucell.magnet.tot_magnetization_nc[2],
+                   ucell.magnet.abs_magnetization};
             break;
         default:
             mag = {};
@@ -446,7 +447,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,

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),

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),

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);
 	}

source/module_elecstate/magnetism.h

@@ -21,7 +21,7 @@ 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

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;
 }

source/module_elecstate/module_charge/charge.h

@@ -70,7 +70,8 @@ class Charge
      * @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,
+    void init_rho(const UnitCell& ucell,
+                  elecstate::efermi& eferm_iout,
                   const ModuleBase::ComplexMatrix& strucFac,
                   ModuleSymmetry::Symmetry& symm,
                   const void* klist = nullptr,
@@ -84,7 +85,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 +100,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,
@@ -144,7 +149,8 @@ class Charge
     void destroy();    // free arrays  liuyu 2023-03-12
 
     bool allocate_rho;
-
+    double omega; // volume of the unit cell
+    double tpiba2; // 2*pi/lat0
     bool allocate_rho_final_scf; // LiuXh add 20180606
 #ifdef __MPI
   private: