Refactor class Charge to template

Author: PeizeLin

source/source_esolver/esolver_double_xc.h

@@ -36,7 +36,7 @@ class ESolver_DoubleXC : public ESolver_KS_LCAO<TK, TR>
     LCAO_domain::Setup_DM<TK> dmat_base;
 
     //! Electorn charge density
-    Charge chr_base;
+    Charge<double> chr_base;
 };
 } // namespace ModuleESolver
 #endif

source/source_esolver/esolver_fp.h

@@ -54,7 +54,7 @@ class ESolver_FP: public ESolver
     K_Vectors kv;
 
     //! Electorn charge density
-    Charge chr;
+    Charge<double> chr;
 
     //! pw_rho: Plane-wave basis set for charge density
     //! pw_rhod: same as pw_rho for NCPP. Here 'd' stands for 'dense',

source/source_esolver/esolver_of.cpp

@@ -217,7 +217,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
         }
 
         delete this->ptemp_rho_;
-        this->ptemp_rho_ = new Charge();
+        this->ptemp_rho_ = new Charge<double>();
 		this->ptemp_rho_->set_rhopw(this->pw_rho);
 		const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
 		this->ptemp_rho_->allocate(PARAM.inp.nspin, kin_den);
@@ -250,7 +250,7 @@ void ESolver_OF::before_opt(const int istep, UnitCell& ucell)
             for (int ibs = 0; ibs < this->pw_rho->nrxx; ++ibs)
             {
                 // Here we initialize rho to be uniform,
-                // because the rho got by pot.init_pot -> Charge::atomic_rho may contain minus elements.
+                // because the rho got by pot.init_pot -> Charge<double>::atomic_rho may contain minus elements.
                 this->chr.rho[is][ibs] = this->nelec_[is] / ucell.omega;
                 this->pphi_[is][ibs] = sqrt(this->chr.rho[is][ibs]);
             }

source/source_esolver/esolver_of.h

@@ -62,7 +62,7 @@ class ESolver_OF : public ESolver_FP
     int tn_spin_flag_ = -1;                       // spin flag used in cal_potential, which will be called by opt_tn
     int max_dcsrch_ = 200;                        // max no. of line search
     int flag_ = -1;                               // flag of TN
-    Charge* ptemp_rho_ = nullptr;                 // used in line search
+    Charge<double>* ptemp_rho_ = nullptr;                 // used in line search
     psi::Psi<double>* psi_ = nullptr;             // sqrt(rho)
 
     // ----------------- used for convergence check -------------------
@@ -90,7 +90,7 @@ class ESolver_OF : public ESolver_FP
     std::function<void(double*, double*)> bound_cal_potential_;
     void cal_potential_wrapper(double* ptemp_phi, double* rdLdphi);
     void cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& ucell);
-    void cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta);
+    void cal_dEdtheta(double** ptemp_phi, Charge<double>* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta);
     double cal_mu(double* pphi, double* pdEdphi, double nelec);
 
     // --------------------- determine the optimization direction -------

source/source_esolver/esolver_of_tool.cpp

@@ -19,7 +19,7 @@ void ESolver_OF::init_elecstate(UnitCell& ucell)
 {
     if (this->pelec == nullptr)
     {
-        this->pelec = new elecstate::ElecState((Charge*)(&chr), this->pw_rho, pw_big);
+        this->pelec = new elecstate::ElecState((Charge<double>*)(&chr), this->pw_rho, pw_big);
     }
 
     delete this->pelec->pot;
@@ -84,7 +84,7 @@ void ESolver_OF::allocate_array()
 
     // initialize chemical potential, step length, ...
     delete this->ptemp_rho_;
-    this->ptemp_rho_ = new Charge();
+    this->ptemp_rho_ = new Charge<double>();
     this->ptemp_rho_->set_rhopw(this->pw_rho);
     const bool kin_den = this->ptemp_rho_->kin_density(); // mohan add 20251202
     this->ptemp_rho_->allocate(PARAM.inp.nspin, kin_den);
@@ -173,7 +173,7 @@ void ESolver_OF::cal_potential(double* ptemp_phi, double* rdLdphi, UnitCell& uce
  * @param [in] ptheta
  * @param [out] rdEdtheta dE/dTheta
  */
-void ESolver_OF::cal_dEdtheta(double** ptemp_phi, Charge* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta)
+void ESolver_OF::cal_dEdtheta(double** ptemp_phi, Charge<double>* temp_rho, UnitCell& ucell, double* ptheta, double* rdEdtheta)
 {
     double* dphi_dtheta = new double[this->pw_rho->nrxx];
 

source/source_estate/elecstate.cpp

@@ -45,7 +45,7 @@ void ElecState::init_scf(const UnitCell& ucell,
 }
 
 
-void ElecState::init_ks(Charge* chr_in, // pointer for class Charge
+void ElecState::init_ks(Charge<double>* chr_in, // pointer for class Charge<double>
                         const K_Vectors* klist_in,
                         int nk_in,
                         const ModulePW::PW_Basis_Big* bigpw_in)

source/source_estate/elecstate.h

@@ -17,7 +17,7 @@ class ElecState
     ElecState()
     {
     }
-    ElecState(Charge* chr_in, ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis_Big* bigpw_in)
+    ElecState(Charge<double>* chr_in, ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis_Big* bigpw_in)
     {
         this->charge = chr_in;
         this->charge->set_rhopw(rhopw_in);
@@ -32,7 +32,7 @@ class ElecState
             this->pot = nullptr;
         }
     }
-    void init_ks(Charge* chr_in, // pointer for class Charge
+    void init_ks(Charge<double>* chr_in, // pointer for class Charge<double>
                  const K_Vectors* klist_in,
                  int nk_in, // number of k points
                  const ModulePW::PW_Basis_Big* bigpw_in);
@@ -112,7 +112,7 @@ class ElecState
 
     int iter = 0;                                  ///< scf iteration
     Potential* pot = nullptr;                      ///< pointer to potential
-    Charge* charge = nullptr;                      ///< pointer to charge density
+    Charge<double>* charge = nullptr;                      ///< pointer to charge density
     const K_Vectors* klist = nullptr;              ///< pointer to k points lists
     const ModulePW::PW_Basis_Big* bigpw = nullptr; ///< bigpw will be removed later
 

source/source_estate/elecstate_lcao.h

@@ -15,7 +15,7 @@ class ElecStateLCAO : public ElecState
     ElecStateLCAO()
     {
     } // will be called by ElecStateLCAO_TDDFT
-    ElecStateLCAO(Charge* chr_in,
+    ElecStateLCAO(Charge<double>* chr_in,
                   const K_Vectors* klist_in,
                   int nks_in,
                   ModulePW::PW_Basis_Big* bigpw_in)

source/source_estate/elecstate_pw.cpp

@@ -13,7 +13,7 @@ namespace elecstate {
 
 template <typename T, typename Device>
 ElecStatePW<T, Device>::ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
-                                    Charge* chr_in,
+                                    Charge<double>* chr_in,
                                     K_Vectors* pkv_in,
                                     UnitCell* ucell_in,
                                     pseudopot_cell_vnl* ppcell_in,

source/source_estate/elecstate_pw.h

@@ -20,7 +20,7 @@ class ElecStatePW : public ElecState
 
   public:
     ElecStatePW(ModulePW::PW_Basis_K* wfc_basis_in,
-                Charge* chr_in,
+                Charge<double>* chr_in,
                 K_Vectors* pkv_in,
                 UnitCell* ucell_in,
                 pseudopot_cell_vnl* ppcell_in,