Refactor : Remove GlobalC::ppcell

source/driver.cpp

@@ -19,9 +19,6 @@ Driver::Driver()
 
 Driver::~Driver()
 {
-    // Release the device memory within singleton object GlobalC::ppcell
-    // before the main function exits.
-    GlobalC::ppcell.release_memory();
 }
 
 void Driver::init()

source/module_elecstate/elecstate.cpp

@@ -208,13 +208,14 @@ void ElecState::calEBand()
 
 void ElecState::init_scf(const int istep, 
                          const ModuleBase::ComplexMatrix& strucfac, 
+                         const bool* numeric,
                          ModuleSymmetry::Symmetry& symm, 
                          const void* wfcpw)
 {
     //! core correction potential.
     if (!PARAM.inp.use_paw)
     {
-        this->charge->set_rho_core(strucfac);
+        this->charge->set_rho_core(strucfac, numeric);
     }
     else
     {

source/module_elecstate/elecstate.h

@@ -110,6 +110,7 @@ class ElecState
      */
     void init_scf(const int istep,
                   const ModuleBase::ComplexMatrix& strucfac,
+                  const bool* numeric,
                   ModuleSymmetry::Symmetry& symm,
                   const void* wfcpw = nullptr);
     std::string classname = "elecstate";

source/module_elecstate/module_charge/charge.h

@@ -84,7 +84,7 @@ class Charge
                     const ModuleBase::ComplexMatrix& strucFac,
                     const UnitCell& ucell) const;
 
-    void set_rho_core(const ModuleBase::ComplexMatrix &structure_factor);
+    void set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric);
     void set_rho_core_paw();
 
     void renormalize_rho();

source/module_elecstate/module_charge/charge_init.cpp

@@ -199,9 +199,7 @@ 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
-)
+void Charge::set_rho_core(const ModuleBase::ComplexMatrix& structure_factor, const bool* numeric)
 {
     ModuleBase::TITLE("Charge","set_rho_core");
     ModuleBase::timer::tick("Charge","set_rho_core");
@@ -249,7 +247,7 @@ void Charge::set_rho_core(
 // each shell of g vec
 //----------------------------------------------------------
             this->non_linear_core_correction(
-                GlobalC::ppcell.numeric,
+                numeric,
                 GlobalC::ucell.atoms[it].ncpp.msh,
                 GlobalC::ucell.atoms[it].ncpp.r.data(),
                 GlobalC::ucell.atoms[it].ncpp.rab.data(),

source/module_elecstate/test/elecstate_base_test.cpp

@@ -60,7 +60,7 @@ 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&)
+void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
 void Charge::set_rho_core_paw()
@@ -249,7 +249,7 @@ TEST_F(ElecStateTest, InitSCF)
     ModuleBase::ComplexMatrix strucfac;
     elecstate->eferm = efermi;
     ModuleSymmetry::Symmetry symm;
-    EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, symm));
+    EXPECT_NO_THROW(elecstate->init_scf(istep, strucfac, nullptr, symm));
     // delete elecstate->pot is done in the destructor of elecstate
     delete charge;
 }

source/module_elecstate/test/elecstate_pw_test.cpp

@@ -139,7 +139,7 @@ K_Vectors::~K_Vectors()
 {
 }
 
-void Charge::set_rho_core(ModuleBase::ComplexMatrix const&)
+void Charge::set_rho_core(ModuleBase::ComplexMatrix const&, const bool*)
 {
 }
 void Charge::set_rho_core_paw()

source/module_esolver/esolver_fp.cpp

@@ -281,7 +281,7 @@ void ESolver_FP::before_scf(UnitCell& ucell, const int istep)
             this->pw_rhod->collect_uniqgg();
         }
 
-        GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, pw_rhod);
+        this->p_locpp->init_vloc(this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
         this->pelec->omega = ucell.omega;

source/module_esolver/esolver_fp.h

@@ -53,6 +53,9 @@ namespace ModuleESolver
 
         ModulePW::PW_Basis* pw_rho;
 
+        //! pointer to pseudopotential
+        pseudopot_cell_vl* p_locpp = nullptr;
+
         /**
          * @brief same as pw_rho for ncpp. Here 'd' stands for 'dense'
          * dense grid for for uspp, used for ultrasoft augmented charge density.

source/module_esolver/esolver_ks.cpp

@@ -68,6 +68,8 @@ ESolver_KS<T, Device>::ESolver_KS()
     ///----------------------------------------------------------
     p_chgmix = new Charge_Mixing();
     p_chgmix->set_rhopw(this->pw_rho, this->pw_rhod);
+    this->ppcell.cell_factor = PARAM.inp.cell_factor;
+    this->p_locpp = &this->ppcell;
 }
 
 //------------------------------------------------------------------------------
@@ -81,6 +83,7 @@ ESolver_KS<T, Device>::~ESolver_KS()
     delete this->pw_wfc;
     delete this->p_hamilt;
     delete this->p_chgmix;
+    this->ppcell.release_memory();
 }
 
 //------------------------------------------------------------------------------

source/module_esolver/esolver_ks.h

@@ -66,6 +66,9 @@ class ESolver_KS : public ESolver_FP
     //! Charge mixing method, only used in KDSFT, not in OFDFT
     Charge_Mixing* p_chgmix = nullptr;
 
+    //! nonlocal pseudo potential
+    pseudopot_cell_vnl ppcell;
+
     //! Electronic wavefunctions
     psi::Psi<T>* psi = nullptr;
 

source/module_esolver/esolver_ks_lcao.cpp

@@ -197,7 +197,8 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
     }
 
     // 8) initialize ppcell
-    GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rho);
+    this->ppcell.init_vloc(this->pw_rho);
+    ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     // 9) inititlize the charge density
     this->pelec->charge->allocate(PARAM.inp.nspin);
@@ -209,7 +210,7 @@ void ESolver_KS_LCAO<TK, TR>::before_all_runners(UnitCell& ucell, const Input_pa
         this->pelec->pot = new elecstate::Potential(this->pw_rhod,
                                                     this->pw_rho,
                                                     &ucell,
-                                                    &(GlobalC::ppcell.vloc),
+                                                    &(this->ppcell.vloc),
                                                     &(this->sf),
                                                     &(this->pelec->f_en.etxc),
                                                     &(this->pelec->f_en.vtxc));
@@ -302,6 +303,7 @@ void ESolver_KS_LCAO<TK, TR>::cal_force(UnitCell& ucell, ModuleBase::matrix& for
                        orb_,
                        force,
                        this->scs,
+                       this->ppcell,
                        this->sf,
                        this->kv,
                        this->pw_rho,
@@ -447,7 +449,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                     this->sf,
                                     *this->pw_rho,
                                     *this->pw_rhod,
-                                    GlobalC::ppcell.vloc,
+                                    this->ppcell.vloc,
                                     *this->pelec->charge,
                                     this->GG,
                                     this->GK,
@@ -474,7 +476,7 @@ void ESolver_KS_LCAO<TK, TR>::after_all_runners(UnitCell& ucell)
                                             this->sf,
                                             *this->pw_rho,
                                             *this->pw_rhod,
-                                            GlobalC::ppcell.vloc,
+                                            this->ppcell.vloc,
                                             *this->pelec->charge,
                                             this->GG,
                                             this->GK,

source/module_esolver/esolver_ks_lcaopw.cpp

@@ -64,7 +64,7 @@ namespace ModuleESolver
     template <typename T>
     void ESolver_KS_LIP<T>::allocate_hamilt()
     {
-        this->p_hamilt = new hamilt::HamiltLIP<T>(this->pelec->pot, this->pw_wfc, &this->kv
+        this->p_hamilt = new hamilt::HamiltLIP<T>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell
 #ifdef __EXX
             , *this->exx_lip
 #endif
@@ -250,7 +250,7 @@ namespace ModuleESolver
                                 *this->pw_wfc,
                                 *this->pw_rho,
                                 *this->pw_rhod,
-                                GlobalC::ppcell.vloc,
+                                this->ppcell.vloc,
                                 *this->pelec->charge,
                                 this->kv,
                                 this->pelec->wg

source/module_esolver/esolver_ks_pw.cpp

@@ -117,7 +117,7 @@ ESolver_KS_PW<T, Device>::~ESolver_KS_PW()
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::allocate_hamilt()
 {
-    this->p_hamilt = new hamilt::HamiltPW<T, Device>(this->pelec->pot, this->pw_wfc, &this->kv);
+    this->p_hamilt = new hamilt::HamiltPW<T, Device>(this->pelec->pot, this->pw_wfc, &this->kv, &this->ppcell);
 }
 template <typename T, typename Device>
 void ESolver_KS_PW<T, Device>::deallocate_hamilt()
@@ -145,7 +145,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                                                                         &(this->chr),
                                                                                         &(this->kv),
                                                                                         &ucell,
-                                                                                        &(GlobalC::ppcell),
+                                                                                        &(this->ppcell),
                                                                                         this->pw_rhod,
                                                                                         this->pw_rho,
                                                                                         this->pw_big);
@@ -156,7 +156,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                                                 &(this->chr),
                                                                 &(this->kv),
                                                                 &ucell,
-                                                                &GlobalC::ppcell,
+                                                                &this->ppcell,
                                                                 this->pw_rhod,
                                                                 this->pw_rho,
                                                                 this->pw_big);
@@ -175,7 +175,7 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
         this->pelec->pot = new elecstate::Potential(this->pw_rhod,
                                                     this->pw_rho,
                                                     &ucell,
-                                                    &GlobalC::ppcell.vloc,
+                                                    &this->ppcell.vloc,
                                                     &(this->sf),
                                                     &(this->pelec->f_en.etxc),
                                                     &(this->pelec->f_en.vtxc));
@@ -194,22 +194,22 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                   &GlobalC::Pkpoints,
                                   GlobalV::MY_RANK,
 #endif
-                                  &GlobalC::ppcell);
+                                  &this->ppcell);
 
     if (this->psi != nullptr)
     {
         delete this->psi;
     }
 
     //! init pseudopotential
-    GlobalC::ppcell.init(ucell.ntype, &this->sf, this->pw_wfc);
+    this->ppcell.init(ucell.ntype, &this->sf, this->pw_wfc);
 
     //! initalize local pseudopotential
-    GlobalC::ppcell.init_vloc(GlobalC::ppcell.vloc, this->pw_rhod);
+    this->ppcell.init_vloc(this->pw_rhod);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "LOCAL POTENTIAL");
 
     //! Initalize non-local pseudopotential
-    GlobalC::ppcell.init_vnl(ucell, this->pw_rhod);
+    this->ppcell.init_vnl(ucell, this->pw_rhod);
     ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL");
 
     //! Allocate psi
@@ -218,7 +218,8 @@ void ESolver_KS_PW<T, Device>::before_all_runners(UnitCell& ucell, const Input_p
                                   this->kv.get_nks(),
                                   this->kv.ngk.data(),
                                   this->pw_wfc->npwk_max,
-                                  &(this->sf));
+                                  &this->sf,
+                                  &this->ppcell);
 
     this->kspw_psi = PARAM.inp.device == "gpu" || PARAM.inp.precision == "single"
                          ? new psi::Psi<T, Device>(this->psi[0])
@@ -247,7 +248,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
     if (ucell.cell_parameter_updated)
     {
-        GlobalC::ppcell.init_vnl(ucell, this->pw_rhod);
+        this->ppcell.init_vnl(ucell, this->pw_rhod);
         ModuleBase::GlobalFunc::DONE(GlobalV::ofs_running, "NON-LOCAL POTENTIAL");
 
         this->pw_wfc->initgrids(ucell.lat0, ucell.latvec, this->pw_wfc->nx, this->pw_wfc->ny, this->pw_wfc->nz);
@@ -256,7 +257,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
 
         this->pw_wfc->collect_local_pw(PARAM.inp.erf_ecut, PARAM.inp.erf_height, PARAM.inp.erf_sigma);
 
-        this->p_wf_init->make_table(this->kv.get_nks(), &this->sf);
+        this->p_wf_init->make_table(this->kv.get_nks(), &this->sf, &this->ppcell);
     }
     if (ucell.ionic_position_updated)
     {
@@ -303,7 +304,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     }
 
     //! calculate the total local pseudopotential in real space
-    this->pelec->init_scf(istep, this->sf.strucFac, ucell.symm, (void*)this->pw_wfc);
+    this->pelec->init_scf(istep, this->sf.strucFac, this->ppcell.numeric, ucell.symm, (void*)this->pw_wfc);
 
     //! output the initial charge density
     if (PARAM.inp.out_chg[0] == 2)
@@ -358,7 +359,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
     // projectors
     ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
 
-    GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
+    this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
 
     // after init_rho (in pelec->init_scf), we have rho now.
     // before hamilt2density, we update Hk and initialize psi
@@ -374,6 +375,7 @@ void ESolver_KS_PW<T, Device>::before_scf(UnitCell& ucell, const int istep)
         this->p_wf_init->initialize_psi(this->psi,
                                         this->kspw_psi,
                                         this->p_hamilt,
+                                        this->ppcell,
                                         GlobalV::ofs_running,
                                         this->already_initpsi);
 
@@ -500,7 +502,7 @@ void ESolver_KS_PW<T, Device>::iter_finish(UnitCell& ucell, const int istep, int
     if (PARAM.globalv.use_uspp)
     {
         ModuleBase::matrix veff = this->pelec->pot->get_effective_v();
-        GlobalC::ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
+        this->ppcell.cal_effective_D(veff, this->pw_rhod, ucell);
     }
 
     if (this->out_freq_elec && iter % this->out_freq_elec == 0)
@@ -623,7 +625,16 @@ void ESolver_KS_PW<T, Device>::cal_force(UnitCell& ucell, ModuleBase::matrix& fo
                            : reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->kspw_psi);
 
     // Calculate forces
-    ff.cal_force(force, *this->pelec, this->pw_rhod, &ucell.symm, &this->sf, &this->kv, this->pw_wfc, this->__kspw_psi);
+    ff.cal_force(force,
+                 *this->pelec,
+                 this->pw_rhod,
+                 &ucell.symm,
+                 &this->sf,
+                 &this->ppcell,
+                 &this->ppcell,
+                 &this->kv,
+                 this->pw_wfc,
+                 this->__kspw_psi);
 }
 
 template <typename T, typename Device>
@@ -641,7 +652,7 @@ void ESolver_KS_PW<T, Device>::cal_stress(UnitCell& ucell, ModuleBase::matrix& s
                            : reinterpret_cast<psi::Psi<std::complex<double>, Device>*>(this->kspw_psi);
     ss.cal_stress(stress,
                   ucell,
-                  &GlobalC::ppcell,
+                  this->ppcell,
                   this->pw_rhod,
                   &ucell.symm,
                   &this->sf,
@@ -781,7 +792,7 @@ void ESolver_KS_PW<T, Device>::after_all_runners(UnitCell& ucell)
     //! 7) Use Kubo-Greenwood method to compute conductivities
     if (PARAM.inp.cal_cond)
     {
-        EleCond elec_cond(&ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &GlobalC::ppcell);
+        EleCond elec_cond(&ucell, &this->kv, this->pelec, this->pw_wfc, this->psi, &this->ppcell);
         elec_cond.KG(PARAM.inp.cond_smear,
                      PARAM.inp.cond_fwhm,
                      PARAM.inp.cond_wcut,