Refactor: Remove LCAO_Deepks from GlobalC. (#5844)

* Remove GlobalC::ld.

* Fix bug for compling without DeePKS.

Author: ErjieWu

source/module_elecstate/elecstate.h

@@ -1,10 +1,9 @@
 #ifndef ELECSTATE_H
 #define ELECSTATE_H
-#include "module_parameter/parameter.h"
-
 #include "fp_energy.h"
 #include "module_cell/klist.h"
 #include "module_elecstate/module_charge/charge.h"
+#include "module_parameter/parameter.h"
 #include "module_psi/psi.h"
 #include "potentials/potential_new.h"
 
@@ -14,10 +13,10 @@ namespace elecstate
 class ElecState
 {
   public:
-    ElecState(){}
-    ElecState(Charge* charge_in,
-              ModulePW::PW_Basis* rhopw_in,
-              ModulePW::PW_Basis_Big* bigpw_in)
+    ElecState()
+    {
+    }
+    ElecState(Charge* charge_in, ModulePW::PW_Basis* rhopw_in, ModulePW::PW_Basis_Big* bigpw_in)
     {
         this->charge = charge_in;
         this->charge->set_rhopw(rhopw_in);
@@ -26,20 +25,20 @@ class ElecState
     }
     virtual ~ElecState()
     {
-        if(this->pot != nullptr) 
+        if (this->pot != nullptr)
         {
             delete this->pot;
             this->pot = nullptr;
         }
     }
-    void init_ks(Charge *chg_in, // pointer for class Charge
-                      const K_Vectors *klist_in,
-                      int nk_in, // number of k points
-                      ModulePW::PW_Basis* rhopw_in,
-                      const ModulePW::PW_Basis_Big* bigpw_in); 
+    void init_ks(Charge* chg_in, // pointer for class Charge
+                 const K_Vectors* klist_in,
+                 int nk_in, // number of k points
+                 ModulePW::PW_Basis* rhopw_in,
+                 const ModulePW::PW_Basis_Big* bigpw_in);
 
     // return current electronic density rho, as a input for constructing Hamiltonian
-    virtual const double *getRho(int spin) const;
+    virtual const double* getRho(int spin) const;
 
     // calculate electronic charge density on grid points or density matrix in real space
     // the consequence charge density rho saved into rho_out, preparing for charge mixing.
@@ -78,17 +77,14 @@ class ElecState
 
     // use occupied weights from INPUT and skip calculate_weights
     // mohan updated on 2024-06-08
-	void fixed_weights(
-			const std::vector<double>& ocp_kb,
-			const int &nbands,
-			const double &nelec);
+    void fixed_weights(const std::vector<double>& ocp_kb, const int& nbands, const double& nelec);
 
-    // if nupdown is not 0(TWO_EFERMI case), 
-    // nelec_spin will be fixed and weights will be constrained 
+    // if nupdown is not 0(TWO_EFERMI case),
+    // nelec_spin will be fixed and weights will be constrained
     void init_nelec_spin();
-    //used to record number of electrons per spin index
-    //for NSPIN=2, it will record number of spin up and number of spin down
-    //for NSPIN=4, it will record total number, magnetization for x, y, z direction  
+    // used to record number of electrons per spin index
+    // for NSPIN=2, it will record number of spin up and number of spin down
+    // for NSPIN=4, it will record total number, magnetization for x, y, z direction
     std::vector<double> nelec_spin;
 
     virtual void print_psi(const psi::Psi<double>& psi_in, const int istep = -1)
@@ -102,7 +98,7 @@ 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
@@ -142,28 +138,24 @@ class ElecState
     void set_exx(const std::complex<double>& Eexx);
 #endif //__LCAO
 #endif //__EXX
- 
+
     double get_hartree_energy();
     double get_etot_efield();
     double get_etot_gatefield();
 
     double get_solvent_model_Ael();
     double get_solvent_model_Acav();
 
-    virtual double get_spin_constrain_energy() {
+    virtual double get_spin_constrain_energy()
+    {
         return 0.0;
     }
 
     double get_dftu_energy();
     double get_local_pp_energy();
 
-#ifdef __DEEPKS
-    double get_deepks_E_delta();
-    double get_deepks_E_delta_band();
-#endif
-
-    fenergy f_en;                                  ///< energies contribute to the total free energy
-    efermi eferm;                                  ///< fermi energies
+    fenergy f_en; ///< energies contribute to the total free energy
+    efermi eferm; ///< fermi energies
 
     // below defines the bandgap:
 

source/module_elecstate/elecstate_energy.cpp

@@ -1,10 +1,10 @@
-#include <cmath>
-
 #include "elecstate.h"
 #include "elecstate_getters.h"
 #include "module_base/global_variable.h"
 #include "module_base/parallel_reduce.h"
 #include "module_parameter/parameter.h"
+
+#include <cmath>
 #ifdef USE_PAW
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include "module_hamilt_pw/hamilt_pwdft/global.h"
@@ -103,10 +103,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
     const double* v_eff = this->pot->get_effective_v(0);
     const double* v_fixed = this->pot->get_fixed_v();
     const double* v_ofk = nullptr;
-    const bool v_ofk_flag =(get_xc_func_type() == 3 
-                            || get_xc_func_type() == 5);
+    const bool v_ofk_flag = (get_xc_func_type() == 3 || get_xc_func_type() == 5);
 #ifdef USE_PAW
-    if(PARAM.inp.use_paw)
+    if (PARAM.inp.use_paw)
     {
         ModuleBase::matrix v_xc;
         const std::tuple<double, double, ModuleBase::matrix> etxc_vtxc_v
@@ -115,19 +114,19 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
 
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
-            deband_aux -= this->charge->rho[0][ir] * v_xc(0,ir);
+            deband_aux -= this->charge->rho[0][ir] * v_xc(0, ir);
         }
         if (PARAM.inp.nspin == 2)
         {
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
-                deband_aux -= this->charge->rho[1][ir] * v_xc(1,ir);
+                deband_aux -= this->charge->rho[1][ir] * v_xc(1, ir);
             }
         }
     }
 #endif
 
-    if(!PARAM.inp.use_paw)
+    if (!PARAM.inp.use_paw)
     {
         for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
         {
@@ -137,16 +136,16 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
         {
             v_ofk = this->pot->get_effective_vofk(0);
             // cause in the get_effective_vofk, the func will return nullptr
-            if(v_ofk==nullptr && this->charge->rhopw->nrxx>0)
+            if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
             {
-                ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr");
+                ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
             }
             for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
             {
                 deband_aux -= this->charge->kin_r[0][ir] * v_ofk[ir];
             }
         }
-        
+
         if (PARAM.inp.nspin == 2)
         {
             v_eff = this->pot->get_effective_v(1);
@@ -157,9 +156,9 @@ double ElecState::cal_delta_eband(const UnitCell& ucell) const
             if (v_ofk_flag)
             {
                 v_ofk = this->pot->get_effective_vofk(1);
-                if(v_ofk==nullptr && this->charge->rhopw->nrxx>0)
+                if (v_ofk == nullptr && this->charge->rhopw->nrxx > 0)
                 {
-                    ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband","v_ofk is nullptr");
+                    ModuleBase::WARNING_QUIT("ElecState::cal_delta_eband", "v_ofk is nullptr");
                 }
                 for (int ir = 0; ir < this->charge->rhopw->nrxx; ir++)
                 {
@@ -219,7 +218,7 @@ double ElecState::cal_delta_escf() const
         if (get_xc_func_type() == 3 || get_xc_func_type() == 5)
         {
             // cause in the get_effective_vofk, the func will return nullptr
-            assert(v_ofk!=nullptr);
+            assert(v_ofk != nullptr);
             descf -= (this->charge->kin_r[0][ir] - this->charge->kin_r_save[0][ir]) * v_ofk[ir];
         }
     }
@@ -278,7 +277,7 @@ void ElecState::cal_converged()
 /**
  * @brief calculate energies
  *
- * @param type: 1 means Harris-Foulkes functinoal; 
+ * @param type: 1 means Harris-Foulkes functinoal;
  * @param type: 2 means Kohn-Sham functional;
  */
 void ElecState::cal_energies(const int type)
@@ -292,7 +291,7 @@ void ElecState::cal_energies(const int type)
     //! energy from gate-field
     this->f_en.gatefield = get_etot_gatefield();
 
-    //! energy from implicit solvation model 
+    //! energy from implicit solvation model
     if (PARAM.inp.imp_sol)
     {
         this->f_en.esol_el = get_solvent_model_Ael();
@@ -305,27 +304,19 @@ void ElecState::cal_energies(const int type)
         this->f_en.escon = get_spin_constrain_energy();
     }
 
-     // energy from DFT+U
+    // energy from DFT+U
     if (PARAM.inp.dft_plus_u)
     {
         this->f_en.edftu = get_dftu_energy();
     }
 
     this->f_en.e_local_pp = get_local_pp_energy();
 
-#ifdef __DEEPKS
-    // energy from deepks
-    if (PARAM.inp.deepks_scf)
-    {
-        this->f_en.edeepks_scf = get_deepks_E_delta() - get_deepks_E_delta_band();
-    }
-#endif
-
     if (type == 1) // Harris-Foulkes functional
     {
         this->f_en.calculate_harris();
     }
-    else if (type == 2)// Kohn-Sham functional
+    else if (type == 2) // Kohn-Sham functional
     {
         this->f_en.calculate_etot();
     }

source/module_elecstate/elecstate_energy_terms.cpp

@@ -3,8 +3,8 @@
 #include "module_elecstate/potentials/efield.h"
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"
-#include "module_hamilt_lcao/module_dftu/dftu.h"
 #include "module_hamilt_lcao/module_deltaspin/spin_constrain.h"
+#include "module_hamilt_lcao/module_dftu/dftu.h"
 
 namespace elecstate
 {
@@ -41,25 +41,15 @@ double ElecState::get_dftu_energy()
 
 double ElecState::get_local_pp_energy()
 {
-    double local_pseudopot_energy = 0.;                   // electron-ion interaction energy from local pseudopotential
+    double local_pseudopot_energy = 0.; // electron-ion interaction energy from local pseudopotential
     for (int is = 0; is < PARAM.inp.nspin; ++is)
     {
-        local_pseudopot_energy += BlasConnector::dot(this->charge->rhopw->nrxx, this->pot->get_fixed_v(), 1, this->charge->rho[is], 1)
-                                * this->charge->rhopw->omega / this->charge->rhopw->nxyz;
+        local_pseudopot_energy
+            += BlasConnector::dot(this->charge->rhopw->nrxx, this->pot->get_fixed_v(), 1, this->charge->rho[is], 1)
+               * this->charge->rhopw->omega / this->charge->rhopw->nxyz;
     }
     Parallel_Reduce::reduce_all(local_pseudopot_energy);
     return local_pseudopot_energy;
 }
 
-#ifdef __DEEPKS
-double ElecState::get_deepks_E_delta()
-{
-    return GlobalC::ld.E_delta;
-}
-double ElecState::get_deepks_E_delta_band()
-{
-    return GlobalC::ld.e_delta_band;
-}
-#endif
-
 } // namespace elecstate

source/module_elecstate/elecstate_print.cpp

@@ -1,13 +1,13 @@
 #include "elecstate.h"
 #include "elecstate_getters.h"
-#include "module_parameter/parameter.h"
 #include "module_base/formatter.h"
 #include "module_base/global_variable.h"
 #include "module_elecstate/potentials/H_Hartree_pw.h"
 #include "module_elecstate/potentials/efield.h"
 #include "module_elecstate/potentials/gatefield.h"
 #include "module_hamilt_general/module_xc/xc_functional.h"
 #include "module_hamilt_lcao/module_deepks/LCAO_deepks.h"
+#include "module_parameter/parameter.h"
 #include "occupy.h"
 namespace elecstate
 {
@@ -311,9 +311,10 @@ void ElecState::print_etot(const Magnetism& magnet,
 
     GlobalV::ofs_running << "\n Density error is " << scf_thr << std::endl;
 
-    if (PARAM.inp.basis_type == "pw") {
+    if (PARAM.inp.basis_type == "pw")
+    {
         ModuleBase::GlobalFunc::OUT(GlobalV::ofs_running, "Error Threshold", pw_diag_thr); // xiaohui add 2013-09-02
-}
+    }
 
     std::vector<std::string> titles;
     std::vector<double> energies_Ry;
@@ -378,7 +379,7 @@ void ElecState::print_etot(const Magnetism& magnet,
         if (PARAM.inp.deepks_scf) // caoyu add 2021-08-10
         {
             titles.push_back("E_DeePKS");
-            energies_Ry.push_back(GlobalC::ld.E_delta);
+            energies_Ry.push_back(this->f_en.edeepks_delta);
         }
 #endif
     }

source/module_elecstate/fp_energy.h

@@ -39,12 +39,13 @@ struct fenergy
     double esol_el = 0.0;  ///< the implicit solvation energy Ael
     double esol_cav = 0.0; ///< the implicit solvation energy Acav
 
-    double edftu = 0.0;       ///< DFT+U energy
-    double edeepks_scf = 0.0; /// DeePKS energy
+    double edftu = 0.0;         ///< DFT+U energy
+    double edeepks_scf = 0.0;   /// DeePKS energy difference
+    double edeepks_delta = 0.0; /// DeePKS energy
 
     double escon = 0.0; ///< spin constraint energy
 
-    double ekinetic = 0.0;  /// kinetic energy, used in OFDFT
+    double ekinetic = 0.0;   /// kinetic energy, used in OFDFT
     double e_local_pp = 0.0; /// ion-electron interaction energy contributed by local pp, used in OFDFT
 
     double calculate_etot();